Analysis of The Brinkman Equation As A Model For Flow in Porous Media

Analysis of the Brinkman equation as a model for flow in porous media
L. Durlofsky and J. F. Brady

Department of Chemical Engineering, California Institute of Technology, Pasadena, California 91125
(Received 21 May 1987; accepted 17 July 1987)
The fundamental solution or Green's function for flow in porous media is determined using
Stokesian dynamics, a molecular-dynamics-like simulation method capable of describing the
motions and forces of hydrodynamically interacting particles in Stokes flow. By evaluating the
velocity disturbance caused by a source particle on field particles located throughout a
monodisperse porous medium at a given value of volume fraction of solids t/J, and by
considering many such realizations of the (random) porous medium, the fundamental solution
is determined. Comparison of this fundamental solution with the Green's function of the
Brinkman equation shows that the Brinkman equation accurately describes the flow in porous
media for volume fractions below 0.05. For larger volume fractions significant differences
between the two exist, indicating that the Brinkman equation has lost detailed predictive value,
although it still describes qualitatively the behavior in moderately concentrated porous media.
At low t/J where the Brinkman equation is known to be valid, the agreement between the
simulation results and the Brinkman equation demonstrates that the Stokesian dynamics
method correctly captures the screening characteristic of porous media. The simulation results
for ;;;.0.05 may be useful as a basis of comparison for future theoretical work.
I. INTRODUCTION
Averaged equations describing viscous flow through porous media are of great theoretical and practical interest. At
the fundamental microscale the Stokes equations apply and
provide a complete description of the entire flow field. However, as a result of the complex and often only statistically
known geometry of the solid surfaces in the medium, solution of the Stokes equations is generally very difficult. On the
macroscopic level, Darcy's law, first established empirically
but more recently derived formally by performing appropriate volume averages of the Stokes equations, is applicable. I-4
The qualitative difference between these two descriptions of
the flow motivated Brinkman5 to suggest a general equation
that interpolates between the Stokes equation and Darcy's
law. His equation,
J.LV 2u- Vp- J.ta 2u = 0,
Vu = 0,
( 1)
2
whereJ.t is the Newtonian fluid viscosity, a- is the permeability, and u andp are the average velocity and pressure, is,
like the Stokes equation but unlike Darcy's law, second order in velocity. This is significant since it allows for the solution of flow around a particle or flow caused by motion of a
particle with no-slip boundary conditions on the surface.
The averages implicit in ( 1) should be viewed as averages
over an ensemble of different realizations of the porous medium.
On small length scales in the Brinkman equation, the
pressure gradient balances the Laplacian of the velocity and
the flow is essentially viscous. Over larger length scales,
where the velocity is slowly varying, the pressure gradient
balances the average velocity as it does in Darcy's law. The
characteristic length that distinguishes between these two
regions of scaling is the Brinkman screening length given by
the square root of the permeability a- 1. In the dilute limit
a- 1 = ({i./3) a- 112, where a is the characteristic particle
3329
Phys. Fluids 30 (11 ), November 1987
size and t/J is the volume fraction of solids. In a viscous fluid,

the velocity disturbance resulting from a point force decays
as 1/r, where r is the distance from the source point to a point
in the fluid. At large distances (r)>a- 1) from a point force
disturbance in a Brinkman medium, however, the response
is very different, with the velocity disturbance decaying faster, as l/(a 2r).
Though Brinkman's derivation ofEq. ( 1) was heuristic,
subsequent investigators have rigorously established the validity of this equation at low volume fraction of solids. 6- 12
The question of the applicability of the Brinkman equation,
i.e., how well it describes flow behavior in a porous medium,
at higher values of t/J remains open, however, although theoretical predictions of permeability based on the Brinkman
equation agree well with experimentally measured values.513 This in itself does not establish its validity (even empirically) for other than dilute systems, as the permeability is
only a single scalar quantity and is not necessarily representative of the general flow field.
In this paper we employ Stokesian dynamics to approximate the fundamental solution, or Green's function, for flow
in random porous media. Stokesian dynamics is a general,
molecular-dynamics-like simulation method capable of describing the motions of, or forces resulting from, hydrodynamically interacting particles immersed in a viscous fluid under conditions of vanishing particle Reynolds number. The
general method applied to finite systems ofhydrodynamically interacting particles has been presented recently by Durlofsky, Brady, and Bossis 14; extensions to infinite systems,
required for the present problem, are considered by Brady et
a/. 15 Once the fundamental solution is determined via Stokesian dynamics it can be compared with the solution of Brinkman's equation, allowing an assessment of the applicability
of the Brinkman equation to porous media of arbitrary volume fraction.
0031-9171 /87/113329-13$01.90
@ 1987 American Institute of Physics
3329
Downloaded 13 Jan 2006 to 131.215.225.172. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
We discuss the application of the general Stokesian dynamics methodology to the present problem in Sec. II. There
it will be seen that the problem is essentially to form the Nparticle mobility matrix, which relates the difference
between the velocity of each particle and the suspension
average velocity to the forces exerted by each particle on the
fluid. Due to the slowly decaying nature of particle interactions in Stokes flow, the effects of particles a great distance
from a given particle must be included. This is accomplished
efficiently using the Ewald summation technique, recently
presented for Stokes flow by Beenakker, 16 and shown to
yield a convergent result by applying the method developed
by O'Brien, 17 as in Brady eta/. 15 Appropriate manipulations
and averages of this mobility matrix result in the fundamental solution as well as the permeability for the system in question.
In Sec. III we compare our simulation results for the
fundamental solution with the Brinkman propagator or
Green's function. Our results are for systems of both point
forces and identical finite-sized spheres, computed both with
and without the application of the Ewald summation technique. For very dilute systems (t,b.;;;;0.01 ), the simulation results for the fundamental solution computed with Ewald
sums agree very well with the Green's function for the Brinkman equation. Because the Brinkman equation is valid as
t,b-+ 0 this is to be expected and serves as a verification of our
method. Discrepancies between the Brinkman and simulation results are evident in simulations performed without
Ewald sums, illustrating the importance, even at low t,b, of
effects from distant particles. In moderately dilute systems,
t,b = 0.05, the fundamental solution is still well described by
the Brinkman propagator, though differences are clearly apparent. The Brinkman result is seen to be no longer quantitatively applicable in moderately concentrated systems t,b
= 0.2, though it still provides a qualitative picture of particle interactions in a porous medium. Finally, we present results for permeability that agree well with the results of
Brinkman5 and Kim and Russel. 13
II. DETERMINATION OF THE FUNDAMENTAL
SOLUTION FOFI FLOW IN POROUS MEDIA
Our intent is to formulate a method that will allow us to

determine the velocity field due to a point force disturbance
in an unbounded porous medium. The disturbance (source)
is applied at a given position and the suspension velocity
(response) is measured at all points (field points) within the
porous medium. By repeating this procedure for many different realizations of the porous medium and performing
appropriate averages, the form of the fundamental sol~tion
for flow in porous media can be deduced. The particles that
comprise the porous medium are prescribed to remain fixed
in space though the forces that they must exert on the fluid to
remain so are unknown. For a given velocity disturbance at
some point in the fluid, the force that any particular particle
must exert to remain fixed depends in part on the positions of
the other particles in the system. Thus the system is fully
coupled; the response of each individual particle is affected
by every other particle.
In all that follows, we suppose that the flow on the parti3330
Phys. Fluids, Vol. 30, No. 11, November 1987
cle scale is in the Stokes regime; i.e., that the particle Reynolds number Ua/v, where Uis a characteristic velocity, a the
characteristic particle size, and v the kinematic viscosity of
the fluid, is much less than unity.
A. Formation of the N-particle mobility matrix

In Stokes flow problems involving systems ofhydrodynamically interacting particles, the velocities of all the particles can be related to the forces exerted on the fluid by each of
the particles via the mobility matrix M:
U = MF,
(2)
where U is the translationaVrotational velocity vector and F

is the force/torque vector, for all N particles in the system.
The inverse problem involves the resistance matrix R, which
relates force/torques to particle velocities;
(3)
F=RU.
The resistance matrix is the inverse of the mobility matrix,

R = M- 1
(4)
Here, M and R depend only on the instantaneous particle

configuration; they are not affected by the velocities or forces
imposed on the particles. In addition, both M and R are
symmetric and positive definite.
General methods have recently been presented for accurately approximating R (and therefore M) for both finite
and infinite systems of hydrodynamically interacting particles.1415 Because these methods are applied to dynamically
evolving configurations of particles, substantial effort is required to assure that short-range lubrication forces, which
act to prevent particles from overlapping during the course
of a dynamic simulation, are correctly included in the resistance matrix. In the present problem, however, only static
(i.e., instantaneous) configurations need be considered, and
our interest is in long-range effects, which are not influenced
by lubrication forces. Thus considerable simplification from
the more general approaches is possible in the present case.
We begin the development of the method by considering
the formulation for a finite system of particles and then pass
to the thermodyna~ic limit; i.e., let the number of particles
Nand the volume of the system Vapproach infinity keeping
the ratioN IV constant. The subsequent development is limited to systems of either point forces or identical spheres of
radius a, although it is straightforward to generalize to more
complex systems. Performing a moment expansion of the
integral representation for the Stokes flow velocity field and
applying Faxen's law, we have the following relationship
between the translational velocity of a given sphere, with
center at xa, and the other N- 1 spheres 14:
Ua-U""(Xa)
61TI-a
81TI-
=--+-X
L (1 + ;
V2) J(xa- Xp)FP,
P=l
P#a
(5)
where Ua is the velocity Of sphere a, U"" (Xa) is the imposed
flow at infinity evaluated at the sphere center, Pis the force
3330
exerted by sphere a on the fluid, and J is the free-space

Green's function or fundamental solution for Stokes flow,
J(r)
= 1/r + rr/~,
(6)
xa -
2
where r =
Xp and r = lrl. For point forces, the (a /
2
3)V J term is identically zero. Writing (5) for each of theN
particles in the system, an approximation to the mobility
matrix ofEq. (2) can be constructed. For systems of finitesized spheres, each sphere-sphere interaction is simply the
well-known Rotne-Prager tensor, the long-range part ofthe
complete two-sphere interaction, evaluated as though the
two spheres were alone in the fluid. Effects from third bodies
do not affect the two-sphere interactions until O(l/r4 ),
where r is a characteristic particle spacing. This is consistent
with the Rotne-Prager approximation, which also neglects
terms of 0( 1/~) in the two-sphere interactions. For systems
of point forces, Eq. ( 5), without the finite size (a 2 /3) V2J
term, is exact; no higher-body effects at all enter in the mobility matrix.
Inversion of the N-particle mobility matrix formed as
described above gives a far-field approximation to the resistance matrix of (3 ). As discussed by Durlofsky eta/., 14 inverting the mobility matrix actually performs all the manybody reflections among all particles. Thus, although the
mobility matrix is formed in a pairwise-additive manner, its
invert, the resistance matrix, contains many-body interactions. In fact, it is these many-body reflections, summed
upon the inversion of the mobility matrix, that give rise to
the screening characteristic of porous media. In the resistance matrix, two-body interactions are via a medium of
fixed (nonzero force) particles, in contrast to the mobility
matrix, where two-body interactions are via a medium of
force-free (nonzero velocity) particles. Therefore, the resistance interactions provide precisely the type of information
required to extract the form of the fundamental solution in
porous media.
The discussion up to this point has been limited to systems of finite numbers of particles. We now consider the
extensions required for infinite systems. To pass to the thermodynamic limit, the number of particles Nand the volume
of the system V approach infinity with the ratioN IV constant. Thus, the volume fraction of particles
(7)
for spherical particles of radius a, can be defined. To simulate an infinite system we could, theoretically, form the mobility matrix as described above for a system of N spheres
and focus only on a subsystem of N 1 spheres immersed within the larger system. For sufficiently large N (N> N 1 >1),
such a system should be representative of an unbounded suspension; in fact, a rigorously convergent expression can be
constructed asN and N 1 - oo. This type of approach is, however, very inefficient computationally; N particles must be
included in the simulation but only N 1 particles provide any
information. Instead of proceeding as described above, we
impose periodic boundary conditions, a technique widely
used in molecular dynamics and Monte Carlo simulations.
This means that we focus on a system of N particles contained within a cell that is periodically replicated throughout
all space. Taking the periodic cell to be a cube of side H, the
3331
volume fraction of particles rp is given by Eq. ( 7) with H 3

replacing V. Although periodic boundary conditions are imposed out of computational necessity, our intent is to model
random systems. Thus the periodic cell must be large enough
that positions of neighboring particles appear random to a
particle placed at the center of the periodic cell. This point
will be discussed in more detail below. In the simulation
results discussed in Sec. III we shall refer to systems of point
forces at nonzero volume fractions. What this means is simply that the particles interact as point forces; i.e., without the
finite-size term in Eq. ( 5). They are still physically spheres
of radius a, and occupy a nonzero volume, with a well-defined volume fraction as given by Eq. (7).
Because the sphere mobility interactions decay slowly,
as 1/r, the effects of spheres far from the test sphere must be
included in the simulation. Therefore, spheres in a given periodic cell must interact not only with spheres in the same
periodic cell but also with the images in other periodic cells.
This entails performing so-called lattice sums of the RotnePrager tensor of Eq. ( 5). Thus for a system of N particles
replicated periodically in space, Eq. ( 5) becomes, with no
impressed flow at infinity,
ua = __!:___ + _1_
61rp,a
81rp,
2
XL fJ=L (1 + a3 V)J<xa -Xp)FP,

y
(8)
where r numbers the unit cells and the double summation is

not performed for r = 1, {3 =a.
Beenakker 16 applied the Ewald summation technique,
previously used to handle slowly decaying Coulombic interactions, to these Rotne-Prager lattice sums, casting them
into a rapidly converging form. In his development, Beenakker specified that the total force on the particles in a unit
cell be zero, and thus the infinite sum in ( 8) is well defined
and convergent. It is essential that the convergence of these
sums be accelerated, as the slow decay of the interaction
would result in very slowly converging sums at great computational cost and questionable accuracy. Writing the Ewaldsummed version of Eq. ( 8) for each of theN particles in the
periodic cell (see Beenakker16 for details), an approximation to the mobility matrix can now be constructed that is
valid for N particles immersed in an unbounded system. We
designate the mobility matrix to which this Ewald summation technique is applied M* and write
U=M*F.
(9)
Note that the mobility matrix in Eq. ( 9), as well as all subsequent mobility and resistance matrices, relate translational
velocity and force, in contrast to the more general matrices
in Eqs. ( 2 )-( 4), which relate translational/rotational velocities to force/torque. Similarly, U now designates the particle translational velocity vector and F the force vector.
When the average force the particles exert on the fluid is
not zero, the expression ( 8) for the velocity of a particle
must be modified. This can be accomplished in a rigorous
fashion by applying a technique first proposed by O'Brien 17
to an infinite suspension of forced particles. The details of
such an approach are in Brady et a/. 15 ; here we shall only
3331
sketch the derivation. We start from an integral representation for the solution to Stokes equations for the velocity field
u(x) at a point x in the fluid in terms of integrals of the force
distribution on the particle surfaces and an integral over a
mathematical surface r of large radius that cuts through
both fluid and particles:
a= I
Sa
1
N
u(x) = - -
81rp,
Jsr
LN
a= I
Sa
1
= --
81rp,
n
-87rp,
( 10)
Ja-n dS
iR J(F)dV.
( 11)
Here, (F) is the average force the particles exert on the fluid,
n = N IV is the number density of particles, and R is the
radius of the volume enclosed by the surface r. The above

reduction is valid for point forces only, but it is straightforward to generalize to the complete problem. 15 It is now permissible to let R ..... oo , because at large distances from x the
sum becomes equivalent to a volume integral of n times the
average force, and ( 11) results in a finite convergent expression for u- (u). Physically, the integral represents a "back
flow" of fluid, relative to zero volume flux axes (u) = 0,
caused by the macroscopic pressure gradient that balances
the excess weight, (F) :1:0, of the particles. It is the velocity
relative to this average back flow that is the physically significant quantity, not its absolute value.
Using the above procedure and Faxen laws for particle
velocities, absolutely convergent expressions for the particle
velocities in Eq. ( 8) can be obtained. With the obvious
changes in notation on the particle sums and expanding the
surface integrals in ( 10) moments, we have in lieu of ( 8)
3332
81rp, r
13 = 1
(1 + a23 v2) J(F)dV.

(12)
Here J is the Green's function (6), K = - 6p,rrr/r,

r = x - y, y being a point on the surface, a is the fluid stress
tensor, and n is the outer normal to the surfaces.
Equation ( 10) is exact for rigid particles. No divergences occur because we have a finite region bounded by the
surface r. This surface is an arbitrary surface immersed in
an unbounded statistically homogeneous suspension, i.e.,
the suspension continues outside of r. If the radius of this
surface is taken to be very large, the variation in J and K will
be small over a surface element dSr that cuts many particles
and fluid. Thus in the integrand of the second integral we
may replace a and u by their suspension averages-fluid and
particle phase averages-( a) and (u). In a statistically homogeneous suspension, (a) and (u) are either constants or
linear functions of position, arising from the average pressure in (a) and a linear shear flow (if one exists) in ( u).
Using the divergence theorem, and introducing the suspension average quantities, ( 10) can be manipulated to give
u(x) - (u(x))
61rp,a
- _n_ ("'
81rp, Jo
Ja-n dS
1
- -- ( (Ja + Ku)n dS.
81rp,
a2)2
ua_ (u(xa)) = -~- + -1. L LN( 1 +-V
The constant term t/J(F) results from the finite particle size
contributions ( a 2/3) V2J. [The reduction of ( 10) for finitesized particles requires some care and is discussed in Brady
et a/. 15 ] Equation ( 12) is an absolutely convergent expression for the velocities of theN particles that have been periodically replicated throughout all space.
It is not necessary for the particles to be periodically
replicated. Equation ( 11) applies for any distribution of particles, but convergence of the difference between the sum and
the integral in ( 11 ) or ( 12) can be accelerated by using periodic replication and the Ewald summation technique. Applying the Ewald summation procedure to (12) results in
exactly the same Ewald-summed mobility matrix M* as in
Eq. (9). The only change the average force makes is to add
the suspension average velocity to the left-hand side. Thus,
in place of ( 9) when the particles are not force-free, we have
U-(u)=M*F.
(13)
Mathematically, the stipulation by Beenakker that

(F) = 0 removes a singular term in the reciprocal space lattice sum at k = 0, where k is a reciprocal wave vector. The
constant term t/J(F) and the back-flow integral in ( 12) precisely cancel this k = 0 term when the average force is not
zero. It may appear surprising that the same mobility matrix
M* is obtained whether or not the average force is zero, but
there is a simple intuitive argument that shows it must be
true. The mobility matrix is a purely geometric quantity that
describes particle interactions. It cannot depend on the velocity, forces, etc., that the particles ultimately have. Therefore, it must be the same whether or not the average force is
zero. Said differently, in writing the mobility matrix in ( 13)
the particles do not know whether the forces sum to zero or
to a finite value; the particle interactions must be the same in
the two cases.
The Ewald summed mobility matrix M* contains all the
information required to compute the fundamental solution
in porous media. Before describing the manipulations that
must be performed to extract the desired results, we shall
briefly consider how M* differs from the far-field mobility
approximations of Durlofsky et al. and Brady et al. In the
present problem, the mobility matrix relates translational
velocities to forces; no angular velocities, imposed rate of
strain, torques, or stresslets are included. In the formulations of Durlofsky et al. and Brady et al., these additional
interactions are included. In the problem at hand, because
no nonuniform bulk flow or external torques are imposed
and because the particles are distributed randomly, the angular velocity, rate of strain, torque, and stresslet unknowns
need not be included. This results in a significant computational savings; M* is reduced from 11N X 11N to 3N X 3N.
3332
We mention in passing that the procedure described

above is fully applicable to ordered as well as random systems. The bulk properties of regular systems are studied in
detail by Brady eta/. 15 through an approach much like that
described here.
B. Fundamental solution for flow In porous media and
comparison with the Brinkman equation
To determine the Green's function for flow in porous

media, we should apply an infinitesimal point force at a point
in the medium and measure the velocity response at all field
points (both fluid and solid phase points since Brinkman's
equation applies to suspension averages). This could be accomplished by using two "test" particles of very small radii
in the mobility formulation ( 13), and then forming the mobility invert to approximate theN-particle resistance matrix
a:
(14)
F = a(U- (u) ).
Recall that interactions in the resistance matrix are via a bed

of fixed particles. The Green's function can then be extracted
from ( 14) as described below.
While correct, the above approach is computationally
expensive because the procedure would have to be repeated
for many different bed particle configurations and for many
different test particle locations in order to form averages.
Instead, we shall simply use two of the bed particles as the
source and field points, applying a force to one and measuring the velocity of the other. For each configuration of bed
particles all possible pairs of particles can be used as source
and field points, and a large amount of statistical information can be obtained from a single realization. Thus, for two
particles a and/3, we form the two-particle resistance matrix
corresponding to their interactions:
~] = [a:a [ F.a
a;,. -
(1/N)(a!a + a:.a)
( 1/N) (a;,. + a~.a)
where a:a and a~.a are the 3 X 3 self-term component matrices and a:.a and a;,. are a-{3 interaction matrices of the Nparticle resistance matrix.
Ifa force is applied to one particle in a porous medium of
infinite extent, the average suspension velocity (u) as well as
the total force exerted by the particles on the fluid, .l: ~
= N (F), are identically zero. Owing to the imposition of
periodic boundary conditions in our model systems, however, the velocity disturbance is periodically replicated
throughout all space, yielding, in general, nonzero values for
the total force and the suspension average velocity. In simulation (u) or (F) must be specified to obtain a well-defined
problem, and thusweseteither (u) =Oor (F) :=0. Note, it is
not possible to prescribe both (u) and (F) because this overdetermines the system of equations. The simplification to
( 12) is obvious in the case (u) =0; for (F) =0 the following
condition on (u) can be derived from Eq. (14):
(16)
for translation of particle 8 with all other particles fixed.
Averaging (16) over all particles 8, (u) can be simply expressed as
(u) = ( 1/N)U,
(17)
where U is the translational velocity of any particle with all

others fixed. Thus, for the case (F) =0, Eq. ( 15) becomes
a:.a- (1/N)(a!a + a:.a>][ua].

a~.a- (1/N)(a~a +a~)
u.a
Equation (18) above, valid for systems in which (F)

:=0, differs from the result for systems with ( u) = 0 by terms
proportional to 1/N. Therefore, the results from two specifications will differ by an amount that tends to zero as N-+ oo,
as would be expected physically. The convergence of the
results from the two specifications with increasing N will be
clearly illustrated in Sec. III. Note that the expressions given
above allow for nonzero (u) even when (F) =0 and nonzero
(F) when (u) =0. These results, though nonintuitive, are
consistent with a macroscopic momentum balance. Specifically, in the case (u) :=0, (F) :;60, a macroscopic pressure
gradient appears which generates an average fluid velocity
that exactly balances the velocity generated by the periodically replicated translating particle. For the case (u) :;60,
(F) :=0, no macroscopic pressure gradient exists, even
though a bulk flow is present, because none is required to
balance the forces exerted by the particles on the fluid.
For a direct comparison of the two-particle simulation
3333
(15)
(18)
result with the Brinkman propagator, the two-particle resistance expressions, Eq. (15) for the case (u) :=Oand Eq. ( 18)
for the case (F) =0, must be recast into mobility expressions, accomplished through inversions. The resulting mobility matrix for either case is designated as M ~.a, where the
superscript Pindicates that the two-particle interaction is via
a porous medium.
The nondimensional velocity field created by a point
force located at xa in a Brinkman medium is
(19)
where u(x) is the velocity at a point x and~ is the force
nondimensionalized by 6rrpa U. The Brinkman propagator
/ is given by Howells8 as
f = (2/a2r)[ (1 + ar + a 2r)e-a'- 1]1

+ (6/a 2r) [1- (1 + ar + j a 2 r)e-a']rr,
(20a)
3333
expressed more concisely as
=~fa (r}l
+ ~ ga (r) (rr/r),
(20b)
where fa (r) and ga (r) are apparent from (20a), with the
factor of~ introduced to cancel the~ in Eq. ( 19). In the limit
tfJ-+ 0, a, nondimensionalized by the sphere radius a, is given
by
(20c)
Note that the velocity disturbance due to a point force in a
Brinkman medium decays far away as l/(a 2r), in contrast
to the l!r decay of the Stokes propagator, which can be recovered from ( 20) in the limit a -+ 0. Also of interest is the
observation that integration of f over a spherical surface
results in an expression that decays exponentially with r
rather than algebraically.
For any isotropic, homogeneous medium, the propagator or fundamental solution can be expressed in a form analogous to Eq. [20(b) ]; therefore the propagator determined
by our simulations for flow in porous media, designated P,
can be written as
P = ~fp (r)l
+ ~ gp (r) (rr/r),
(21)
where the velocity field due to a point force is given by replacing f with Pin Eq. (19). By determiningjp(r) and
gp(r), the porous medium propagator is fully specified.
From our simulations we determine M ~. whose off-diagonal component matrix relates the velocity of a field particle
located at Xp to the force exerted by a source particle at xa ,
and is therefore the porous medium propagator P. If all the
particles in the simulation are point forces, the off-diagonal
component matrix of M ~ gives the Green's function for
flow in such a porous medium; on the other hand, if the
particles are all finite-sized spheres, M ~P relates the velocity
of sphere {3 to the force exerted by sphere a.
The two functions fp (r) and gp ( r) are determined by
evaluating the velocity disturbance caused by a source particle as measured by the motion of the field particles. In the
actual simulation, the following two summations, evaluated
at discrete values of r, where r is the distance from the source
particle to the field particle, are performed:
m,
m1
p= I
0-t/J)-LU~--
41rr
=
1
m2
(1-t/J)-
Lm,
p
=I
=!
L L
a- 2
= ~ {1/t/J )K.
(24)
In the actual simulations, particles are placed randomly,

via random sequential addition, in the periodic cell. Particles
are not permitted to overlap. As mentioned above, it is important that the periodic cell be large enough that our results
are not influenced by the imposition of periodic boundary
conditions. In a Brinkman medium, the screening length,
beyond which sphere interactions can be expected to differ
significantly from those in a Stokes fluid, is the square root of
the permeability, a- 1 To ensure that the model system is
indeed sufficiently large, we require that the ratio of the periodic cell size, H [cf. Eq. (7) ], to the Brinkman screening
length be large compared to unity; e.g.,
aH;::;:,3.4N 113tfJ 116 > 1.
(25)
In most of our simulations, N = 125; t/J varies from 0.002

( aH;:::;:, 6) to 0.2 ( aH;:::;:, 13). Thus, as we shall see in Sec. III,
our system is indeed large enough to display screening.
All of the simulations discussed in the next section were
performed on a Cray X-MP supercomputer. In simulations
involving the application of the Ewald summation technique, lattice sums are computed over 125 periodic cells, and
approximately six minutes of CPU time are required. Over
90% of this time is spent computing the lattice sums; thus
simulations of 125 particles without the application of the
Ewald summation technique require less than 30 sec of CPU
time.
u 1 dS
(22a)
undS
hemisphere
[fp(r) +gp(r)]Ff,
(22b)
where m 1 is the total number of particles contained within a

spherical shell of thickness !J..r situated a distance r from the
source particle, and m 2 is the total number of particles contained within a hemispherical shell oriented such that its
base is in the 2-3 plane. The subscript 1 refers to the 1 direction, n is the unit normal, and UP is the velocity of particle /3.
The factors of ( 1 - tfJ) correct for the fact that all field points
sampled are considered fluid elements, i.e., the particle {3 is
assumed force-free in ( 15) in computing the summations in
(22), while in reality a fraction tfJ are fixed in space and
would have given zero to the sum. By evaluating each of
3334
1 (1
N
N
)
(23)
K-'=-trR!p ,
3
N a=IP=!
where R!p, nondimensionalized by 61T'p,a, is as defined
above. The dimensionless permeability a- 2 is related to the
drag coefficient by
sphere
[fp(r) +!gp(r)] Ff,
Upn-21T'r
these sums for imposed forces in the 1, 2, and 3 directions

and by considering each particle in tum as the source particle, a large body of statistical information is generated, yielding accurate estimates for fp ( r) and g P ( r).
The average force exerted by a fixed particle immersed
in a uniform flow ( nondimensionalized by 611'p,a U), referred
to as the drag coefficient K - 1, can be determined directly
from Eq. ( 14 ). Specifically, the drag coefficient is given by
Ill. RESULTS
In this section we present simulation results for particle
interactions in porous media, computed both with and without the application of the Ewald summation technique, and
compare them with the Brinkman solution at four different
values of t/J: O.Dl, 0.002, 0.05, and 0.2. The comparison is in
terms ofthe functionsf(r) andg(r) and the drag coefficient
or resistivity K - 1 All distances are nondimensionalized by
the sphere radius a. The results presented at a given value of
t/J are averaged over three independent, random realizations.
In each realization every particle in tum is considered as the
source particle, yielding a large body of statistical information. At given values of t/J and N, several sets of simulations
are generally performed to gauge the differences between
results computed with different system specifications; e.g.,
finite-sized spheres versus point forces or (u) =0 vs (F) =0.
3334
In these cases, the same three realizations are used for each
set of simulations. The variations between the results of the
three realizations are generally quite small. In thefp (r) and
gp (r) results, the standard deviations of the three-realization average are usually smaller than the size ofthe symbols
on the figures by the third or fourth data point after r = 2
(typically rz 3.5 to rz 5) and remain so for all larger values
of r. For the first few data points, the standard deviation is
generally about 10% of the average for both fp (r) and
gp ( r). More specific information concerning the variations
between the three realizations will be cited only in those
cases where the variations differ significantly from the general observations cited above.
We first consider results for= 0.01. Figures 1(a) and
1(b) display the simulation results, shown as X 'sand + 's,
forfp(r) andgp(r) for a system of 125 point forces. In this
and all subsequent figures the X 's correspond to simulations
of point forces for which (u) =0 and the + 's to simulations
of point forces for which (F) =0. The results in Figs. 1(a)
and 1(b) are for simulations performed with the application
of the Ewald summation technique. The solid lines are the
Brinkman medium results, given by Eqs. (20), with a as
given by the infinite dilution result Eq. (20c), and the
dashed lines are pure fluid (Stokes propagator) results;
0.15
f 5 (r)
=g5 (r) = (i)r- 1 Thehalf-boxwidth (H /2 = 18.7)
is indicated by a vertical line on the abscissa. Both sets of

simulation results agree quite well with the Brinkman medium results and deviate from one another only slightly and by
a nearly uniform offset. As discussed in Sec. II, the two sets
of simulation results deviate from one another by an 0( 1/N)
amount; a larger variation between them will be evident in
theN= 27 results shown below.
Note that the simulation results agree with the Brinkman propagator over the entire range of r, up to r-;:::;32. At
the parameter values for these simulations ( = 0.01, N
= 125) the box size H is 37.4. The maximum distance
between two particles is ,jjH /2 ( = 32.4), the distance from
the center of the periodic cell to a comer. Because portions of
spherical surfaces of radius r extend beyond the periodic cell
(a cube of side H), for values of r beyond H /2, the summations in Eq. ( 19) are not actually over the surfaces of spheres
of radius r but over the portion of the surface of a sphere of
radius r that falls within a cube of side H. This geometrical
limitation results in a diminution in the amount of statistical
information as r ..... ,jjH /2, but does not affect the validity of
Eq. ( 19) or the agreement between the simulation results
and the Brinkman equation in the region H /2 < r < ,J3H /2.
Because the Brinkman equation is rigorously valid for a
random, infinite system as J ..... O, we expect our simulation
I
I
<I>
=0.01
0.15~~~--~---.-~-,...---,-----,
I
I
0.12
I
I
(a)
0.09
f (r)
''
''
0.06
' ' ........
........... ....
0.03
(a)
\
0.09
-------------
f(r)
' ' ',

''....................
0.06
-- ------ ---
0.03
<jo=O.OI
\
0.12
---------
xXxxxxxxxxxxxxxxxxxxxxx)(
-0.032
12
17
22
27
32
7
17
12
22
27
32
<j>=O.OI
0.26
(b)
I
I
I
I
I
I
I
g(r)
0.14
0.08
(b)
I
I
I
g(r)
\
''
' ......
... ... _
0.02
----
I
\
0.)4
0.08
',
' '' .............. ... _
-----------)(
0.02
-0.04
12
17
22
27
FIG. I. Comparison of 1 = 0.01 theoretical and simulation results for (a)

/(r)and (b) g(r). Simulations performed for 125 point forces with the application of the Ewald summation technique. The X 's are (u) =0 results;
+ 's are (F) =0 results, the solid curves are the Brinkman propagator
fs (r) and g 8 (r) functions, and the dashed curves are the Stokes propagator. The vertical line on the abscissa in these and all subsequent figures indicates the half-box width.
)(
)(
;-;-)(')(l<'>c'-Jt-f!C..-,r
32
7
3335
12
17
22
27
32
FIG. 2. Comparison of 1 = 0.01 theoretical and simulation results for (a)

/(r) and (b) g(r). Simulations performed for 125 point forces without the
application of the Ewald summation technique. The X 's are (u) sO simulation results, the solid curves are the Brinkman propagator, and the dashed
curves are the Stokes propagator. Note the deviation of the simulation and
Brinkman results for r > H /2.
3335
results to recover the Brinkman propagator at low 4> for sufficiently large N. Figures 1(a) and 1(b) indicate that a periodically replicated random system of 125 point forces does
indeed behave as a Brinkman medium when the particles are
specified to be fixed in space. This can be taken as verification that the inversion of the mobility matrix M* results in
all many-body scatterings (at least at the level of point
forces), yielding a medium that behaves fundamentally differently than either a pure fluid or a suspension afforce-free
particles.
We next present results for simulations of point forces at
the same parameter values as above (t/> = 0.01, N = 125),
but without the application of the Ewald summation technique. In these simulations, periodic boundary conditions
are imposed but no lattice sums are performed; i.e., particles
interact only with their nearest neighbors and not with an
infinite replication of images. The results forfp (r) andgp (r)
areshowninFigs. 2(a) and 2(b). Again the X'scorrespond
to point-force simulation results with (u) =0, the solid lines
to the Brinkman propagator and the dashed line to the
Stokes propagator. Though agreement with the Brinkman
solution is good for both f~ctions in the range 2<r< 17,
agreement beyond r-;:::; 20 becomes poor, particularly for
gp(r), which tends to the Stokes solution for r>27. Comparison of Figs. 2(a) and 2(b) with Figs. 1(a) and 1(b)
clearly illustrates the importance of the effect of distant particles, even for 4> as low as 0.01. As will be seen below, the
effect of distant particles is similarly important for 4> = 0.05.
Having established that theN= 125 point-force simulations at 4> = 0.01 with the application of Ewald sums do indeed reproduce Brinkman's result, the effect of the size of the
periodic cell will now be considered. Shown in Figs. 3 (a)
and 3(b) arefp (r) andgp (r) for systems of27 point forces
for 4> = 0.01 (H /2 = 11.2) for simulations performed with
Ewald sums. Again, the X 's correspond to (u) =0 and the
+ 's to (F) = 0. The standard deviations in the fp ( r) and
gp (r) data points for the N = 27 simulation results are
slightly higher than those observed with N = 125. lngp (r),
standard deviations of 0(0.01 to 0.02), for both the (u) =0
and (F) =0 simulations, exist in the second, third, and sixth
data points (r = 3.45, 4.42, 7.33 ). At larger values of r, the
standard deviations for the N = 27 simulation fp (r) and
gp(r) functions are, though small [-0(0.001), smaller
than the size of the symbols] , about a factor of 3 larger than
those for theN= 125 simulation functions.
It is apparent from Figs. 3 (a) and 3 (b) that, though the
trends agree with the Brinkman result, the offsets between
the simulation and Brinkman results and between the two
simulation results are noticeably greater than for the N
= 125 simulations [compare Figs. l(a) and l(b)].
Further, these variations are statistically significant. As
shown in Sec. II the two sets of simulation results deviate
from one another by an 0( liN) amount. To determine the
variation of the offset between the Brinkman and simulation
results with system size, additional simulations with (u) =0
were performed at 4> = 0.01 for systems of 64 and 90 point
forces. The offset is quantified by computing the average of
the difference between the simulation and Brinkman results
over the range H /2 < r < .,f3H /2. Least squares fits for log3336
oj>=O.OI
',',
f (r)
(a)
' ......... .....
........ .............
------------
++++++
..
++++
14
17
-0.03
oj>= 0.01
+I
\
0.2
20
(b)
\
\
g (r)
0.14
' ', ....

~
0.08
........
........
-- -- ----- ----
0.02
-0.04
14
17
20
FIG. 3. Comparison of <,6 = 0.01 theoretical and simulation results for (a)
f(r) and (b) g(r). Simulations performed for 27 point forces with Ewald
sums. The X 's are (u) =0 simulation results, the + 's are (F) =0 simulation results, the solid curves are the Brinkman propagator, and the dashed
curves are the Stokes propagator. The deviations between the two simulation results for /(r) and between the simulation results and the Brinkman
function are noticeably larger than for theN= 125 simulations [see Fig.
l(a)].
log plots of offset in both fp ( r) and g P ( r) vs N for

27 ..;;N..;; 125 give slopes of - 0.94 0.05 and - 0. 77
0.03, respectively. Thus over this range of N the average
offset scales approximately as N - 0 86 Only considering the
range 64.-;N.-;125, the slopes are - 1.03 0.13 and
- 0.87 0.01 for the offsets infp (r) andgp (r), respectively, giving an average offset that scales approximately as
N - 0 95 Therefore, it appears that the difference between the
Brinkman and simulation results scales approximately as
N - I , as does the offset between the two sets of simulation
results. Further, it is evident from Figs. 1(a) and 3(a) that
for simulations with (F) =O,fp(r) approaches/B (r) from
above and for simulations with (u) =0 the approach is from
below. Such a trend does not exist for the gp ( r) functions.
We next briefly consider simulation results for very dilute systems: 4> = 0.002. These simulations, performed for
125 point forces with Ewald sums and (u) =0, are of interest
because they allow comparison over a large range of r (HI
2 = 32.0). Figures 4(a) and 4(b) display the results for
fp(r) andgp(r); agreement with the Brinkman solution is
excellent over the entire range of the results, r < 54. Note the
only slight difference between the Stokes and Brinkman g ( r)
functions in Fig. 4(b) and the agreement of the simulation
3336
0.2 r r r - - - , - - - - - - - , - , - - - - - , - - - - - ,
I
I
I
(a)
I
I
f ( r)
\
\
''
0.04
',, ... ,
"x
------ --------
xxxxxxxxxxxxxxxxxx
-0.0 4 2L _ _ ___..LI5----2L8---'---l.41------='54
0.34 . . - - - - - . - - - - - - - , . . . - - - , - - - - - . - - - - - - - ,
0.28
(b)
0.22
where/ is given by Eq. (20a). The (aB2/B0 ) 2 V2 f term

in (26a) arises because V2 f satisfies the Helmholtz equation, giving V4 / =a 2V2/ . 18 Note that V2/ decays exponentially with r, so its effect is only noticeable at small r.
Further, the dilute limit approximation for a, Eq. ( 20c), will
not be used but rather the value for a obtained from the
simulation [ cf. Eq. ( 24) ] .
The results for fp(r) and gp(r) for simulations performed with Ewald sums are displayed in Figs. 5(a) and
5(b). The X 's correspond to point force, (u) =0 simulation
results; the + 's to point force, (F) =0 simulation results;
O's to finite sized sphere, ( u) = 0 simulation results; the solid
curves to the Brinkman propagator / 8 (r) and g8 (r) functions computed using the point force simulat;i!-:>n result for a
(aPf = 0.6775), and the broken curves to the / 8 (r) and
g8 (r) functions, corresponding to f 0 > in Eq. (26), computed using the finite-sized sphere simulation result for a
(a.r. = 0.6348). By rz6, the solid and broken curves for
both / 8 (r) and g 8 (r) are in good agreement despite the
slightly different a values.
0.12 r - - - - - , - - - - - - - - - , - - , - , - - , - - - - - - - - ,
oj>0.05
0
(a)
FIG. 4. Comparison of tfi = 0.002 theoretical and simulation results for {a)
/(r) and {b) g{r). Simulations performed for 125 point forces with Ewald
sums. The X 's are (u) =0 simulation results, the solid curves are the Brinkman propagator, and the dashed curves are the Stokes propagator.
results with the latter. Despite the only slight quantitative

differences between the simulation (or Brinkman) results
and the pure fluid (Stokes) results, it is important to emphasize that the simulation and Brinkman/(r) and g(r) functions decay as 1/(a 2 ~) at larger while the pure fluid functions decay as 1/r. This fundamental difference can be
readily appreciated if one computes the spherical average of
the Brinkman and Stokes propagators: the Brinkman propagator behaves as e-ar/r, while the Stokes propagator still
only decays as 1/r.
We now tum from very dilute systems to a consideration of moderately dilute porous media: <fJ = 0.05. The simulation results presented at <fJ = 0.05 are both for systems of
125 point forces and 125 identical finite-sized spheres. When
considering systems of finite-sized spheres, comparison of
the simulation results with the Green's function for the
Brinkman equation may no longer be appropriate. Rather,
the simulation results should be compared to the Brinkman
equation equivalent of the Rotne-Prager tensor for Stokes
flow; i.e., the Brinkman Green's function plus twice the Faxen law contribution, 8
/(1)
{1 +[~2 +(~2rl v2} /.
B 0 (a) = 1 +a
B 2 (a)
3337
(ea-
+!a2,
B0 )/a
(26a)
(26b)
(26c)
_l
-0.04 2:-----::6-----:,70---'--------:'",4:------:,8
0.36r-----,------,----,----,----
(b)
0.28
g (r)
0
*
0
*~
~ ~
1/J
-0.04 2~----;!-6--------:,~0--'------;';,4:----~,8
FIG. 5. Comparison of tfi = 0.05 theoretical and simulation results for {a)
/(r) and {b) g{ r). Simulations performed for 125 particles with the application of Ewald sums. The X 's are ( u) = 0 point force simulation results, the
+ 's are (F) =0 point force simulation results, the O's are (u) =0 simulation results for finite-sized spheres, the solid curves are Brinkman propagator functions (a= 0.6775), and the broken curves are Brinkman RotnePrager functions {a = 0.6348). The deviations between the simulation and
Brinkman results indicates the loss of accuracy of the Brinkman equation at
tfi=0.05.
3337
Though quantitative agreement between the Brinkman

functions and the simulation results does not exist except at
large values of r, qualitatively the results are in fair agreement even at moderate values of r. Both the simulation and
Brinkman results, for both point forces and finite-sized
spheres, display negative minimums in/(r) at moderate values of r. The fp ( r) results for systems of finite-sized spheres
lie above those for systems of point forces, just as the/B (r)
functions of Eq. ( 26) lies above the Brinkman propagator
fB (r) function. In the g(r) functions both the simulation
results for systems of finite-sized spheres and the Brinkman
Rotne-Prager result level off near r = 2 while the simulation
results for systems of point forces and the Brinkman propagator decrease monotonically. In all cases, however, quantitative differences are apparent between the simulation and
Brinkman results.
It is our belief that the discrepancies between the simulation and Brinkman results at t/J = 0.05 are real and will persist as N-+ oo. The small deviations between the (u) =0 and
(F) =0 point force results indicate that the simulation results have essentially converged for increasing N, and considerable quantitative differences still exist between the simulation and Brinkman results, particularly ing(r). Further,
simulations performed with systems of27 point forces do not
indicate that the simulation results tend toward the Brinkman results with increasing N. For simulations in which
(F) :::0, theN= 27 fp (r) results lie consistently above the
N = 125 fp (r) results. At r = 2.28, where the largest absolute deviation for r<; 10 between the two occurs ( -0.01), the
discrepancy in theN= 27 results relative to theN= 125
results is 16%. The fact that the minimum in fp (r) for
N = 27 simulations lies above that for N = 125 simulations
suggests that the minimum becomes more negative, and thus
deviates more from the Brinkman result, as N increases. The
gp (r) results for the two simulations deviate only slightly; at
r = 2.28 the N = 27 simulation result exceeds that of the
N = 125 simulation result by 1.5%, while at r = 3.82 the
N = 27 result falls below theN= 125 result by 6%. Thus,
there is no evidence that the t/J = 0.05 point force simulation
results tend to the Brinkman solution with increasing N.
The Brinkman equation provides an exact description of
flow in random porous media as t/J-+0. Rigorous derivation
of the Brinkman equation requires only point force interactions among the fixed particles-no higher moments are
necessary. Therefore, the quantitative differences between
the point force simulation results and the Brinkman propagator at t/J = 0.05 are significant; they indicate that the
Brinkman equation has begun to lose accuracy for systems of
point forces at t/J as low as 0.05.
The next simulation results presented are for t/J = 0.05,
N = 125 with (u) :::0 but without the application of the
Ewald summation technique. Simulations of systems of
point forces under these conditions resulted in nonpositive
definite mobility matrices for some realizations, while simulations of systems of finite-sized spheres yielded well-behaved mobility matrices in all cases. Results for fp (r) and
gp (r) for finite sized spheres are shown in Figs. 6(a) and
6(b). The open circles are the simulation results and the
solid and broken curves are the Brinkman propagator and
3338
0.12
(a)
0.08~
0.04
f(r)
0
\
K~
0
0
0
-0.04
0.05
0
0
0
0
0
0
-0.08
I
2
10
18
14
0.36
~0.05
0.28
(b)
g(rl
oooooooooooooooo
-0.0 4 2.'----~6------,I'::-O--'------:'cl4:-------:'18
FIG. 6. Comparison of~= 0.05 theoretical and simulation results for (a)
f(r) and (b) g(r). Simulations performed for 125 finite-sized spheres without Ewald sums. The O's are ( u) = 0 simulation results, the solid curves are
Brinkman propagator functions (a= 0.6775), and the broken curves are
Brinkman Rotne-Prager functions (a= 0.6348). Thef(r) simulation results are aphysical for r > H /2.
Brinkman Rotne-Prager results, respectively, with the same

a's as in Figs. 5(a) and 5(b). (As discussed below, simulations performed without the application of the Ewald sum
technique do not yield reasonable a values.) The results for
fp (r) agree qualitatively with the Ewald summed t/J = 0.05
simulation results [Fig. 5 (a) 1for 2 < r < 10, but beyond this
the trend is entirely different. For r>H /2, thefp (r) simulation results are rather aphysical; they increase in magnitude
with increasing separation. The gp (r) results appear more
reasonable, though they deviate considerably from the results of Fig. 5(b) and actually tend more toward the Stokes
solution [not shown in Fig. 6(b) 1than the Brinkman solution at larger. This was also the case for the t/J = 0.01 gp (r)
results for simulations performed without Ewald sums [see
Fig. 2(b) 1
The last set of simulation results to be considered is for
moderately concentrated suspensions, t/J = 0.2. In simulations with N = 125, the mobility matrix lost positive definiteness in all cases when the Ewald summation technique
was not applied (for systems of either point forces or finitesized spheres) and also when the Ewald summation technique was applied to systems of point forces. Thus the only
results presented are for systems of 125 finite-sized spheres
simulated with Ewald sums. Thefp ( r) and gp ( r) results for
3338
these simulations are shown in Figs. 7(a) and 7(b). The

open circles correspond to the ( u) = 0 simulations, the asterisks to (F) =0 simulation results, the solid curves to the
Brinkman propagator using the value for a obtained from
the simulation (a= 1.989), and the broken curves to the
Brinkman Rotne-Prager functions using the same value for
a. Thestandarddeviationsinthefp(r) andgp(r) simulation
results are very small for all values of r-in all cases smaller
than the size of the symbols. The simulation and Brinkman
results deviate significantly, though the trends are still in
qualitative agreement. Note that the deviations between the
two are similar to those observed at t/J = 0.05; thefp (r) simulation results display a more negative minimum than do the
/ 8 (r) curves, while thegp (r) simulation results lie considerably abovetheg8 (r) curves at small to moderate values ofr.
The t/J = 0.2 simulation results should be interpreted
more as an indication of the loss of the validity of the Brinkman equation than as an accurate estimate of the actual porous media/(r) andg(r) functions. At t/J = 0.2, characteristic nearest-neighbor interparticle spacings r 12 [r12 = (~1T/
t/J) 113 ] are 2.76 sphere radii, to be contrasted with 7.48, 12.8,
and 4.38 at t/J values of 0.01, 0.002, and 0.05, respectively.
The relatively near spacing of the spheres at t/J = 0.2 indicates that higher moments in the expansion of the integral
0.01
4> 0.2
. . .
(a)
*
0
0
0
f (r)
0
-0.03
-0.04
10
12
0.16
4> 0.2
0.13
0.1
g(r)
!!
!!
(b)
expression for the Stokes flow velocity field may be required;

in other words, the Rotne-Prager tensor may no longer suffice as a description of two-sphere mobility interactions.
Further, lubrication may need to be introduced into theresistance matrix to accurately account for very short-ranged
interactions. These extensions are available and straightforward, 1415 but require considerably more computation time
than the present method and were therefore not implemented. Also, the use of bed particles rather than infinitesimal
test particles may affect the behavior near r = 2 at this high
value of t/J. The t/J = 0.2 simulation results may find additional use as a "known" result to serve as a comparison with
a theoretical treatment of nondilute porous media where
sphere-sphere interactions are approximated by the RotnePrager tensor.
Finally, we consider the results for the drag coefficient
K - 1, for the simulations discussed above. Reasonable results
for the drag coefficient are obtained only in simulations performed with the application of Ewald sums. Without Ewald
sums, for 125 particles, K - 1 = 0.114 at tfJ = 0.01 for point
forcesandK - I = 0.0712at(J = 0.05 forfinite-sizedspheres.
These results are clearly unrealistic; K - 1 must exceed unity
and increase with increasing t/J.
Displayed in Table I are drag coefficient results with
simulations performed with Ewald sums. The standard deviations presented correspond to variations between the
averages for each of the three distinct realizations. In the last
column are the results of Kim and Russel 13; for t/J<0.05 their
low t/J expansion result [their Eq. (2.14)] is presented and
for t/J = 0.2 their numerical value. Several trends are apparent from the simulation results. Most importantly, K - I increases with increasing t/J; this is not the case when the Ewald
summation technique is not applied. Values forK - 1 are consistently higher for simulations of systems of point forces
than for simulations of systems of finite-sized spheres, presumably a result of the slightly different nature of the manybody reflections that occur upon inversion of the mobility
matrix for a system of finite-sized spheres compared to those
for a system of point forces.
The K - I values presented in Table I increase with decreasing system size, although this trend does not appear to
be statistically significant; i.e., the standard deviation also
increases with decreasing system size, as would be expected.
!!
!t
I!
TABLE I. Results for the drag coefficient K - 1 defined in Eq. ( 24). All
simulations were performed with the application of the Ewald summation
technique.
0.04
0.02 z'----....J4'-----'6-.....!....---'-B---..LI0-----'12
FIG. 7. Comparison of 1 = 0.2 theoretical and simulation results for (a)

/(r) and (b) g(r). Simulations performed for 125 finite-sized spheres with
Ewald sums. The O's are (u) =0 simulation results, the s are (F) =0 simulation results, the solid curves are Brinkman propagator functions
(a= 1.989), and the broken curves are Brinkman Rotne-Prager functions
(a= 1.989). The Brinkman equation is not valid at 1 = 0.2.
3339
Point force
or
finite size
0.002
0.01
0.01
0.01
0.05
0.05
0.05
0.2
125
125
125
27
125
125
27
125
PF
PF
FS
PF
PF
FS
PF
FS
K-1
(K-I)a
1.098 0.024
1.259 0.044
1.237 0.043
1.320 0.082
2.040 0.032
1.791 O.D35
2.178 0.126
4.396 0.056
1.102
1.280
1.280
1.280
1.981
1.981
1.981
4.61
Results from Kim and Russel. 13

3339
For very small systems, however, a statistically significant

increase in the drag coefficient would be expected to occur.
Due to the imposition of periodic boundary conditions, a
small system is not truly random and would behave somewhat as an ordered, periodic system. Reference to Kim and
Russel 13 (their Table III) shows that the drag coefficient is
consistently greater (at least for <0.45) for ordered as
compared to random systems.
Finally, we note that the agreement between the simulation results for the drag coefficient and the results of Kim
and Russel is quite good; within 10% in all cases. This would
be expected for low c,6, but it is somewhat surprising that the
agreement persists for c,6 as large as 0.2. At c,6 = 0.2, the
Brinkman equation does not provide a quantitative basis for
computing sphere interactions in a porous medium, as is
evident from Figs. 7(a) and 7(b). Although the= 0.2simulation results for fp(r) and gp(r) are only approximate
(see discussion above), the drag coefficient result is expected
to be more reliable, due to the relatively small effect oflubrication in the calculation of the drag coefficient (refer to
Brady eta/. 15 for a discussion of this point). In their permeability calculations, however, Kim and Russel approximate
sphere-sphere interactions as via a Brinkman medium at all
values of c,6. Though this approximation loses accuracy at
moderate and high c,6, there are apparently compensating
effects in their calculations that render the final result for the
permeability accurate.
IV. CONCLUSIONS
In this paper, we have applied the general methodology
of Stokesian dynamics to determine the form of the fundamental solution for flow in porous media. In simulations of
dilute porous media, the results clearly show that the system
behaves as a Brinkman medium, with long-ranged interactions screened by intervening fixed particles, rather than as a
viscous fluid. This "effective medium" behavior of the simulated porous medium comes about upon the inversion of a
properly constructed N-particle mobility matrix, which itself derives from a moment expansion of the integral representation of the Stokes velocity field. Thus the effective properties of the medium arise naturally out of the Stokesian
dynamics methodology; they need not be postulated a priori.
Indeed, we have presented our analysis of the Brinkman propagator as if it were obvious that Stokesian dynamics necessarily gives the correct answer. In retrospect, it is obvious
that Stokesian dynamics is correct, and we hope the present
study provides a convincing proof. In a subsequent publication, we shall show that the effective interactions among particles in sedimenting suspensions can be determined in much
the same way as those of porous media using the Stokesian
dynamics method. 19
The results presented in Sec. III demonstrate the agreement between the simulation and the Brinkman equation at
low c,6, but quantitative differences between the two, indicating the loss of validity of the Brinkman equation, are evident
for ;;;.0.05. The results also show the importance of including the effects of distant particles via the Ewald summation
technique; at c,6 = 0.01 and c,6 = 0.05 qualitative inaccuracies
3340
appear at large distances from the source particle in simulations performed without Ewald sums. Over more restricted
regions (r<H /2), however, the simulations performed
without the application of Ewald sums do provide reasonable results. It may seem surprising that the Brinkman equation starts to break down at what appears to be a rather low
volume fraction, but it should be realized that a characteristic interparticle spacing at c,6 = 0.05 is only slightly larger
than four particle radii; particles are actually rather close
together. Permeability or drag coefficient calculations only
yield realistic values when the Ewald summation technique
is applied; in these cases agreement with the self-consistent
results of Brinkman5 and the results of Kim and Russel 13 is
consistently good up to c,6 = 0.2, the highest value of c,6 considered.
In addition to presenting the form of the fundamental
solution for flow in porous media at various values of c,6, the
simulation results at c,6 = 0.05 and c,6 = 0.2 may find additional use as a basis of comparison for future theoretical
work. Rubinstein 12 rigorously derived the Brinkman equation by considering a dilute system of fixed spheres approximated as point forces, but his diluteness criterion is highly
restrictive. In another paper, Rubenstein20 suggests that subsequent theoretical approaches, aimed at rigorously describing porous media at higher values of c,6, may involve higher
multipoles. In this case, the theoretical results may be directly comparable to our simulation results for systems of finitesized spheres.
The procedure presented in this paper is appropriate for
studying at most moderately concentrated suspensions. Extensions to more concentrated systems are developed and
applied to the study of bulk properties (the nature of the
fundamental solution is not considered) of ordered systems
by Brady et a/. 15 ; disordered systems are presently under
study. For disordered systems, a Monte Carlo method is applied to assure that hard-sphere distributions are obtained;
the simple random sequential addition used in the present
study reaches a percolation threshold at moderate values of
volume fraction and subsequently fails. Further, at high volume fractions the full Stokesian dynamics method, which
includes lubrication interactions, must be used.
ACKNOWLEDGMENTS
We wish to thank J. Lester and R. Phillips for their

assistance with some of the computer programming and G.
Bossis for many useful discussions.
Computer time was provided on the San Diego Supercomputer Center Cray X-MP through a National Science
Foundation allocation. Partial support for this work was
through National Science Foundation Grant No. CBT8696067.
H. P. G. Darcy, Lesfontanes publiques de Ia ville de Dijon (Dalmont,

Paris, 1856).
2
J. B. Keller, in Statistical Mechanics and Statistical Methods in Theory and
Applications, edited by U. Landman (Plenum, New York, 1977), pp. 631644.
3340
S. Whitaker, Transport in Porous Media 1, 3 ( 1986).

H. Brenner, Philos. Trans. R. Soc. London Ser. A 297, 81 ( 1980); P.M.
Adler and H. Brenner, Phys. Chern. Hydrodyn. 5, 245 (1984).
~H. C. Brinkman, Appl. Sci. Res. A 1, 27 (1947).
6
C. K. W. Tam, J. Fluid Mech. 38, 537 (1969).
7
S. Childress, J. Chern. Phys. 56, 2527 (1972).
8
1. D. Howells, J. Fluid Mech. 64,449 ( 1974).
9
E. J. Hinch, J. Fluid Mech. 83, 695 ( 1977).
1
K. F. Freed and M. Muthukumar, J. Chern. Phys. 68, 2088 ( 1978).
"M. Muthukumar and K. F. Freed, J. Chern. Phys. 70, 5875 (1979).
4
3341
12
J. Rubenstein, J. Fluid Mech. 170, 379 (1986).

S. Kim and W. B. Russel, J. Fluid Mech. 154, 269 ( 1985).
14
L. Durlofsky, J. F. Brady, and G. Bossis, J. Fluid Mech. 180, 21 ( 1987).
'~J. F. Brady, R. Phillips, J. Lester, and G. Bossis, submitted to J. Fluid
Mech.
16
C. W. J. Beenakker, J. Chern. Phys. 85, 1581 (1986).
17
R. W. O'Brien, J. Fluid Mech. 91, 17 ( 1979).
18
S. Kim and W. B. Russel, J. Fluid Mech. 154,253 (1985).
19
1. F. Brady and L. Durlofsky, submitted to Phys. Fluids.
20
1. Rubenstein, J. Stat. Phys. 44, 849 (1986).
13
3341

Analysis of The Brinkman Equation As A Model For Flow in Porous Media

Uploaded by

Copyright:

Available Formats

Analysis of The Brinkman Equation As A Model For Flow in Porous Media

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Analysis of The Brinkman Equation As A Model For Flow in Porous Media

Uploaded by

Copyright:

Available Formats

Analysis of the Brinkman equation as a model for flow in porous media

L. Durlofsky and J. F. Brady

J.LV 2u- Vp- J.ta 2u = 0,

Phys. Fluids 30 (11 ), November 1987

size and t/J is the volume fraction of solids. In a viscous fluid,

@ 1987 American Institute of Physics

Our intent is to formulate a method that will allow us to

Phys. Fluids, Vol. 30, No. 11, November 1987

A. Formation of the N-particle mobility matrix

where U is the translationaVrotational velocity vector and F

The resistance matrix is the inverse of the mobility matrix,

Here, M and R depend only on the instantaneous particle

V2) J(xa- Xp)FP,

exerted by sphere a on the fluid, and J is the free-space

Phys. Fluids, Vol. 30, No. 11, November 1987

volume fraction of particles rp is given by Eq. ( 7) with H 3

XL fJ=L (1 + a3 V)J<xa -Xp)FP,

where r numbers the unit cells and the double summation is

n = N IV is the number density of particles, and R is the

radius of the volume enclosed by the surface r. The above

(1 + a23 v2) J(F)dV.

Here J is the Green's function (6), K = - 6p,rrr/r,

Phys. Fluids, Vol. 30, No. 11, November 1987

Mathematically, the stipulation by Beenakker that

We mention in passing that the procedure described

To determine the Green's function for flow in porous

Recall that interactions in the resistance matrix are via a bed

where U is the translational velocity of any particle with all

a:.a- (1/N)(a!a + a:.a>][ua].

Equation (18) above, valid for systems in which (F)

Phys. Fluids, Vol. 30, No. 11, November 1987

f = (2/a2r)[ (1 + ar + a 2r)e-a'- 1]1

expressed more concisely as

In the actual simulations, particles are placed randomly,

aH;::;:,3.4N 113tfJ 116 > 1.

In most of our simulations, N = 125; t/J varies from 0.002

where m 1 is the total number of particles contained within a

[fp(r) +!gp(r)] Ff,

these sums for imposed forces in the 1, 2, and 3 directions

Phys. Fluids, Vol. 30, No. 11, November 1987

=g5 (r) = (i)r- 1 Thehalf-boxwidth (H /2 = 18.7)

is indicated by a vertical line on the abscissa. Both sets of

' ' ........

' ' ',

' '' .............. ... _

FIG. I. Comparison of 1 = 0.01 theoretical and simulation results for (a)

FIG. 2. Comparison of 1 = 0.01 theoretical and simulation results for (a)

Phys. Fluids, Vol. 30, No. 11, November 1987

' ......... .....

' ', ....

log plots of offset in both fp ( r) and g P ( r) vs N for

where/ is given by Eq. (20a). The (aB2/B0 ) 2 V2 f term

results with the latter. Despite the only slight quantitative

{1 +[~2 +(~2rl v2} /.

Phys. Fluids, Vol. 30, No. 11, November 1987

Though quantitative agreement between the Brinkman

Phys. Fluids, Vol. 30, No. 11, November 1987

Brinkman Rotne-Prager results, respectively, with the same