7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Instantaneous particle acceleration model for gassolid suspensions at moderate
Reynolds numbers
S. Tenneti*, R. O. Foxtand S. Subramaniam*
Department of Mechanical Engineering, Iowa State University, Ames, IA 50011, USA
t Department of Chemical and Biological Engineering, Iowa State University, Ames, IA 50011, USA
sudheert@iastate.edu, rofox@iastate.edu, shankar@iastate.edu
Keywords: Immersed boundary method, DNS, drag laws, granular temperature
Abstract
Gassolid flows are encountered in many industrial applications such as fluidized beds and coal gasification. The
design and scaleup of such industrial devices requires a better understanding of the characteristics of gassolid
suspensions. Devicescale computational fluid dynamics (CFD) simulations that solve for average quantities such as
solid volume fraction and phasic mean velocity fields are being extensively used in the industrial design process. The
capability of these simulations to accurately predict the characteristics of gassolid flow depends upon the accuracy
of the models for unclosed terms that appear in the equations for mass, momentum and energy conservation. Hrenya
and Sinclair (1997) show that the particle granular temperature (particle velocity variance) plays an important role
in the prediction of the core annular structure in riser flows. In statistically homogeneous suspensions undergoing
elastic collisions, the particle accelerationvelocity covariance alone governs the evolution of granular temperature.
This accelerationvelocity covariance can be decomposed into a source and dissipation of granular temperature due
to hydrodynamic forces. Koch and coworkers (Koch 1990; Koch and Sangani 1999) quantified the hydrodynamic
source and dissipation terms in the granular temperature evolution using a combination of kinetic theory closure
and multiple expansion simulations at very low Reynolds numbers (Stokes flow regime). At moderate Reynolds
numbers, particle resolved direct numerical simulation (DNS) is as a viable tool to quantify the hydrodynamic source
and dissipation. In this study, DNS of freely evolving gassolid suspensions are performed using the Particleresolved
Uncontaminatedfluid Reconcilable Immersed Boundary Method (PUReIBM) that has been developed at Iowa State
University to simulate flow past fixed particle assemblies (Garg et al. 2009) and freely evolving suspensions (Tenneti
et al.). Analysis of DNS results shows that the fluctuations in the particle acceleration that are aligned with the
fluctuations in the particle velocity give rise to source in the granular temperature. It is found that simple extension of
a class of mean particle acceleration models to their corresponding instantaneous versions does not predict the correct
joint accelerationvelocity statistics that are obtained from DNS. Also such models do not give rise to any source
in the granular temperature due to hydrodynamic effects. This motivates the development of better instantaneous
particle acceleration models. It is found that a Langevin equation for the increment in the particle velocity reproduces
the DNS results for particle velocity autocorrelation in freely evolving suspensions.
Introduction
Particleladen flows are very common in many indus
trial applications like fluidized bed combustors, coal
gasification, pneumatic transport lines etc. In addition to
the industrial significance of gassolid flows, the com
plexity of the physical mechanisms that govern these
flows attracts the attention of many researchers. Further
more, a fundamental understanding of gassolids flow is
relevant due to increasing interest in carbonneutral en
ergy generation technologies such as chemical looping
combustion (Li and Fan 2008).
Computational fluid dynamics (CFD) simulations that
solve for the averaged equations of multiphase flow are
a costeffective solution for rapid evaluation of design
and scale up of industrial devices like fluidized beds.
Devicescale CFD simulations are usually based on the
EulerianEulerian twofluid approach in which averaged
equations for conservation of mass, momentum and en
ergy are written for each phase, with coupling terms rep
resenting the interphase interactions. These equations
contain unclosed terms that need to be modeled. For
example, the mean momentum conservation equation in
the particle phase requires closure of the average fluid
particle interaction force (mean drag force) and the av
erage stress in the solid particle phase. Therefore, the
predictive capability of multiphase CFD simulations de
pends on the models for interphase exchange of species,
momentum, and heat. In this work we only consider clo
sure models for the interphase momentum transfer.
In any statistical closure problem, an important mod
eling question is the adequacy of the mathematical rep
resentation to capture physical phenomena. Hrenya and
Sinclair (1997) show that it is necessary to solve the
transport equation for the particle granular temperature
to predict the coreanular structure observed in riser
flows. This shows that closure only at the level of the
mean momentum is not adequate, but a closure at the
level of second moment of particle velocities is neces
sary. But it is not clear if closure at the level of second
moments is adequate for predictive CFD simulations.
The objective of this work is to use direct numeri
cal simulations (DNS) to propose a particle acceleration
model that provides a closure for all moment equations.
In order to achieve this objective we consider the evo
lution of the oneparticle distribution function. The ad
vantage of this approach is that closure at the level of
the oneparticle distribution function implies a closure
for all the moment equations. Tenneti et al. discussed
in detail the evolution equation of the oneparticle dis
tribution function for gassolids flow and analyzed ex
isting models used to close the conditional expectation
of the particle acceleration term that appears in the evo
lution equation. They showed that simple extension of
mean particle acceleration models does not predict the
the correct particle accelerationvelocity covariance ob
tained from DNS. A consequence of the use of such
models in CFD codes is that the granular temperature
will decay to zero resulting in physically incorrect be
havior and triggering numerical instabilities.
For statistically homogeneous gassolids flow with
perfectly elastic solid particles, the particle granular
temperature evolves only due to the correlation be
tween the fluctuating particle acceleration and the fluc
tuating particle velocity. For dilute suspensions in the
Stokes flow regime Koch (1990) decomposed the par
ticle accelerationvelocity covariance into source and
sink contributions arising due to hydrodynamic interac
tions and derived the corresponding analytical expres
sions for these terms. Later Koch and Sangani (1999)
used an approximate multiple expansion method to
specify the source and sink terms for higher solid vol
ume fractions in the Stokes flow regime. At moder
ate mean flow Reynolds numbers, Tenneti et al. re
ported DNS results that confirm Koch's observation (in
the Stokes flow regime) that particle velocity fluctua
tions correlate with particle acceleration fluctuations to
generate a source in the granular temperature.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
In this study, DNS of freely evolving gassolid
suspensions are performed using the Particleresolved
Uncontaminatedfluid Reconcilable Immersed Bound
ary Method (PUReIBM) that has been developed at Iowa
State University to simulate flow past fixed particle as
semblies (Garg et al. 2009) and freely evolving suspen
sions (Tenneti et al.). PUReIBM is a particleresolved
direct numerical simulation approach for gassolid flow
with the following features that distinguish it from other
immersed boundary method approaches:
1. Uncontaminated fluid: In PUReIBM the immersed
boundary (IB) forcing is solely restricted to those
grid points that lie in the solid phase, and therefore
the flow solution in the fluid phase is uncontami
nated by the IB forcing. Consequently the velocity
and pressure in the fluid phase is a solution to the
unmodified NavierStokes equations (in contrast to
IB implementations that smear the IB forcing on to
grid points in the fluid phase adjoining solid bound
aries, resulting in solution fields that do not corre
spond to unmodified NavierStokes equations).
2. Reconcilable: In PUReIBM the hydrodynamic
force experienced by a particle is computed directly
from the stress tensor at the particle surface that is
obtained from this uncontaminated fluid flow so
lution (in contrast to IB implementations that cal
culate the hydrodynamic force from the IB forc
ing field). This feature of PUReIBM enables us
to directly compare the DNS solution with any
randomfield theory of multiphase flow. In par
ticular, for statistically homogeneous suspensions
it is shown that (Garg et al. 2009) if the volume
averaged hydrodynamic force exerted on the par
ticles by the fluid is computed from a PUReIBM
simulation, it is a consistent numerical calcula
tion of the average interphase momentum transfer
term 7' )6 (x x('))) in the twofluid the
ory (Drew 1983). This reconciles DNS results with
multiphase flow theory.
Owing to these specific advantages, it is shown else
where (Garg et al. 2009; Garg 2009) that PUReIBM is
a numerically convergent and accurate particleresolved
DNS method for gassolids flow. Its performance has
been validated in a comprehensive suite of tests:(i)
Stokes flow past simple cubic (SC) and face centered
cubic (FCC) arrangements (ranging from dilute to close
packed limit) with the boundaryintegral method of Zick
and Homsy (1982), (ii) Stokes flow past random ar
rays of monodisperse spheres with LBM simulations
of van der Hoef et al. (2005) (iii) moderate to high
Reynolds numbers (Re, < 300) in SC and FCC ar
rangements with LBM simulations of Hill et al. (2001)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and (iv) high Reynolds number flow past random arrays
of monodisperse spheres with ANSYSFLUENT CFD
package. It has also been extended to study passive
scalar transport, and validated for heat transfer from a
single isolated sphere (Garg 2009; Garg et al.). In this
work, the DNS methodology based on PUReIBM devel
oped by Tenneti et al. to simulate freely evolving sus
pensions is used to propose an instantaneous particle ac
celeration model that incorporates the effect of particle
velocity fluctuations and hydrodynamic effects of neigh
boring particles.
The rest of the paper is organized as follows. First the
instantaneous particle acceleration model is described.
Results from the simulations of freely evolving suspen
sions and verification of the Langevin model are dis
cussed next. Finally a procedure for determining the co
efficients of the acceleration model is described.
Instantaneous particle acceleration model
We propose the following stochastic model for the incre
ment in the particle velocity:
dvo = 3 (W,) dt yv}'dt + BdW,. (1)
The above equation is an isotropic form of the general
Langevin model developed by Tenneti and Subrama
niam. In the above equation dvi is the increment in the
particle velocity, v' is the fluctuation in the particle ve
locity and i \V is a Wiener process increment. Fluctua
tions in the particle velocity are defined about the mean
particle velocity, i.e. v = vj (vj). The first term on
the right hand side of Eq. 1 accounts for the effect of
the mean slip velocity. The mean slip velocity, defined
as (W) (v) (u(f)), is the relative velocity between
the solid phase mean velocity and the fluid phase mean
velocity. The second term accounts for the fluctuation in
the particle velocity and the last term models the effect
of the hydrodynamic interaction of the neighboring par
ticles. The coefficient is the inverse of the Lagrangian
particle velocity autocorrelation time. It quantifies how
long a particle retains memory of its initial velocity.
These coefficients are functions of volume fraction (y),
mean flow Reynolds number (Re,) and particle to fluid
density ratio (pp/pf). To extract a functional form for
the Langevin model coefficients, simulations of freely
evolving suspensions where, the motion of the particles
is affected by the surrounding fluid, are performed using
the DNS methodology developed by Tenneti et al.. Re
sults from the simulations freely evolving suspensions
are presented in the following section.
101
100
101
I
0
102
103
30 40
tllid
Figure 1: Plot showing the evolution of the Reynolds
number based on granular temperature (ReT)
for a freely evolving suspension of elastic par
ticles at a volume fraction of 0.1 and solid to
fluid density ratio of 100.
Results
IBM simulations of freely evolving suspensions are per
formed for solid volume fractions between 0.1 and 0.4
and for mean flow Reynolds numbers between 10 and
100. The Reynolds number based on the mean slip ve
locity between the fluid and particulate phase is defined
as
pf I (v) (fu')) d
Re, (1 ) (, (2)
P[f
where p is the solid volume fraction, pf and pf are the
density and dynamic viscosity of the fluid phase respec
tively, and d is the particle diameter. When characteriz
ing the effect of particle velocity fluctuations it is useful
to define a Reynolds number based on the granular tem
perature ReT as:
pfdT1/2
where T is the granular temperature which is given by
1
T ( vi"7). Estimation of granular temperature from
3
DNS of freely evolving suspensions is discussed in de
tail by Tenneti et al.. The evolution of ReT for a freely
evolving suspension at a volume fraction of 0.1 and solid
to fluid density ratio of 100 is shown in Fig. 1. Firstly,
we observe that for all mean flow Reynolds numbers,
the granular temperature attains a statistical steady state.
From the figure we can also see that at a given solid to
fluid density ratio, the steady state granular temperature
increases with increasing mean flow Reynolds number.
This result can be explained based on the fact that the
particles pick up energy from the fluid and the energy in
the system increases with increasing Re,. To the au
thors' knowledge, this is the first report of the effect of
mean flow Reynolds number on granular temperature.
Similar behavior of the steady granular temperature with
Rem is observed for all the volume fractions studied (0
= 0.1, 0.2, 0.3 and 0.4). The Langevin model for the
particle acceleration is verified by computing the parti
cle velocity autocorrelation function after the granular
temperature reaches a steady state.
Verification of Langevin model. In this section, verifi
cation of the Langevin model for the instantaneous par
ticle acceleration is presented. For this purpose we con
sider the increment in the particle velocity fluctuations:
dv' = yv'dt + B. 1V (4)
As described earlier, in the above equation the model co
efficient 7 is the inverse of the integral time scale of the
particle velocity autocorrelation. The particle velocity
autocorrelation function p (s) is defined as follows:
p ((s (to) v' (to + s))
( (to). '(to))
where s is the separation in time. The integral time
scale for the autocorrelation function is defined as TL
fo p (s) ds. Using this definition, we computed the in
tegral time scale from DNS after the granular tempera
ture reached a steady state. If a stochastic process obeys
the Langevin equation with an integral time scale of TL,
then its autocorrelation function should decay exponen
tially, i.e., p (s) = e/TL. We extracted the autocorre
lation function from the DNS and compared it with the
exponential function predicted by the Langevin model.
This comparison is shown in Fig. 2. We can see that for
both density ratios considered, the evolution of the au
tocorrelation function obtained from DNS matches with
the exponential decay predicted by the Langevin model.
We established that Langevin model predicts the dynam
ics of a freely evolving suspension very well after the
granular temperature reaches steady state. However, at a
given volume fraction, mean flow Reynolds number and
solid to fluid density ratio, we need to specify the coef
ficients as a function of the granular temperature so that
we can predict the evolution of the suspension using the
Langevin model. In order to do this, we have to identify
the source and dissipation of granular temperature from
the DNS data. In the next section, a method to iden
tify source and dissipation of granular temperature from
DNS is presented.
Identification of source and dissipation from DNS.
Using the Langevin model for the increment in the par
ticle velocity fluctuations (cf. Eq. 2), we can derive the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
evolution equation for the modeled granular temperature
T* as:
dT*
S 27T* +B2. (5)
dt
In the above equation, we can clearly identify that the
source for the granular temperature is B2 and that the
dissipation is 27T*. For a statistically homogeneous
suspension, the evolution equation for the granular tem
perature can be written as:
dT 2
d (A"iv"i) .
dt 3
In the above equation, the fluctuations in the accel
eration are defined about the mean particle accelera
tion, i.e. A" Aj (Aj). The instantaneous par
ticle acceleration model should model evolution of the
granular temperature correctly. In order to do this, we
have to match the source and dissipation implied by the
Langevin model to the source and dissipation obtained
from DNS. However, given the correlation (A 'v'), it
is nontrivial to uniquely decompose it into source and
dissipation. Koch (1990) derived analytical expressions
for the source and dissipation in the limit of low volume
fractions and low Reynolds numbers. Later Koch and
Sangani (1999) used approximate multiple expansions
approach to determine the source and dissipation for
dense suspensions but limited to the Stokes flow regime.
Here we present a method to extract the same from the
DNS at moderate Reynolds numbers.
The fluctuation in the acceleration experienced by the
1' "' particle is denoted A"(") and similarly, the fluctu
ation in the velocity is denoted v"( ). The fluctuating
acceleration can be written as:
A" = )v + A" i (6)
In the above equation, we decomposed the fluctuating
acceleration vector along a direction parallel to the par
ticle fluctuating velocity and along a direction perpen
dicular to it. The component of the acceleration per
pendicular to the fluctuating velocity is denoted A" (
to represent the remainder term. It is important to note
that this is not a model but an exact expression for the
fluctuating particle acceleration. We can now form the
estimate for the correlation (A' 'v') as follows:
1 1p
( A 1v\ I: A"."v
(A E V1 NP n=l
After substituting Eq. 4 into the above equation and per
forming some algebraic manipulations, we can write the
evolution equation for the estimate of granular tempera
ture as:
dT S (7)
dt
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1
0.8
.0.6
CP
0.4
0.2
0 0.5 1 1.5
s/T,
2 2.5
Figure 2: Plots showing the verification of the Langevin model. Figure 2a compares the particle velocity autocorre
lation function extracted from the DNS of freely evolving suspension (volume fraction of 0.2, mean flow
Reynolds number 20 and solid to fluid density ratio of 100) with the exponential decay predicted by the
Langevin model. Figure 2b is the same as Fig. 2a for a suspension with solid to fluid density of 100.
where the source is,
Np
S _2 1 ( ( ,) !!( ) (8)
Sn=l
3N, (8)
and the dissipation is
r C v v" (9)
Sn= 1
In the above expressions for source and dissipation,
) 1 (() + () ) and (((r (r ).
From these equations, it is clear that the particles
whose fluctuating acceleration is aligned with the fluc
tuating velocity contribute to the source in granular
temperature. Particles with the fluctuating acceleration
aligned in a direction opposite to that of the fluctuating
velocity contribute to the dissipation in granular temper
ature. This can be easily visualized from the scatter plot
shown in Fig. 3. For illustration, in this figure we show
the scatter plot of fluctuating particle acceleration and
fluctuating particle velocity obtained from the DNS of
flow past a fixed particle assembly (0 = 0.2, Rem = 20
and ReT = 16). The symbols that lie in the first and
the third quadrants denote the particles whose fluctuat
ing acceleration is aligned with the fluctuating veloc
ity. Hence, these particles contribute to the source in
granular temperature. Similarly, the symbols in the sec
ond and fourth quadrants contribute to the dissipation in
granular temperature. Tenneti et al. show that a simple
extension of mean particle acceleration model applied
to any particle velocity distribution does not produce any
scatter in the first and third quadrants. They demonstrate
that such instantaneous acceleration models do not pro
duce any source in the granular temperature evolution
equation. On the other hand, a stochastic acceleration
model (cf. Eq. 1) always results in a finite source term
for the particle granular temperature for nonzero B.
We wish to know if the steady state attained by the
suspension is a stable attractor. We can verify this by
studying suspensions with different initial conditions.
Using the formulae given by Eq. 6 and Eq. 7, we ex
tracted the source and dissipation from the DNS data and
plotted them versus the granular temperature in Fig. 4 as
a phase space plot. The state of the system is character
ized by the granular temperature T and in phase space
dT
plots, the rate of change of the state variable dt is plot
ted versus the state variable T. But the rate of change
of granular temperature is simply the difference between
the source (S) and dissipation (F). To understand the be
havior of source and dissipation separately, we plot both
S and F versus T in the plots shown in Fig. 4. Source
and dissipation are plotted for two sets of simulations
differing only in their initial conditions. For both cases,
the solid volume fraction is 0.2, mean flow Reynolds
number is 20 and solid to fluid density ratio is 100. Tri
angular symbols are the data for a suspension where the
particles initialized with zero granular temperature. In
this case, the particles pick up energy from the fluid and
hence we note that the source term is greater than dissi
pation at initial time. Square symbols show the source
DNS
Langevin Model
. ..
C.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 3: Scatter plot showing the the fluctuating parti
cle acceleration versus the fluctuating particle
velocity obtained from the DNS of flow past
a fixed particle assembly at a volume fraction
of 0.2, mean flow Reynolds number of 20 and
Reynolds number based on granular tempera
ture of 16. From the analysis on the extraction
of source and dissipation from the DNS, we
can see that the symbols that lie in the first and
the third quarter contribute to source and the
symbols that lie in the second and the fourth
quadrant contribute to the dissipation.
and dissipation data for a suspension initialized with a
granular temperature corresponding to ReT = 4. In this
case, the particles have higher energy and thus they lose
their energy to fluid. So for this case the dissipation is
initially greater than the source term. A key observation
from the figure is that both suspensions attain the same
steady granular temperature. This shows that the steady
state behavior of a statistically homogeneous gassolids
suspension is independent of the initial temperature. Fi
nally we note that for both suspensions, the source and
dissipation are equal at steady state. These observations
verify the expressions for the source and dissipation pre
sented earlier and show that this is a viable approach to
propose the acceleration model.
Conclusions
Evolution of the particle granular temperature in gas
solid suspensions at moderate Reynolds numbers is
studied by DNS of freely evolving suspensions using
PUReIBM. DNS results show that the steady granu
lar temperature increases with the mean flow Reynolds
number. This steady granular temperature is indepen
Figure 4: Phase space plot showing the variation of
source and dissipation with granular temper
ature. DNS data for two different suspensions
is shown here. The volume fraction, mean
flow Reynolds number and solid to fluid den
sity ratio are 0.2, 20 and 100 respectively. Tri
angles denote the source (filled symbols) and
dissipation (hollow symbols) and dissipation
for a suspension initialized with zero granu
lar temperature. Squares denote the data for
a suspension initialized with a finite granular
temperature corresponding to ReT = 4.
dent of the initial granular temperature of the suspen
sion. A key finding that emerges from this work is that
at steady state the particle velocity autocorrelation func
tion decays exponentially thus motivating the use of a
Langevinlike model to describe the increment in parti
cle velocity. A Langevin equation is proposed to model
the instantaneous particle acceleration in gassolid sus
pensions at moderate Reynolds numbers. Expressions
to compute the source and sink of granular temperature
from DNS are derived. Previous work by Tenneti et al.
shows that simple extension of mean particle acceler
ation models does not recover the joint acceleration
velocity statistics that are obtained from DNS. Using
the expressions for source and dissipation derived in this
work it is also shown that such models do not produce
any source for the particle granular temperature which
justifies the use of Langevinlike models to model the
instantaneous particle acceleration.
Acknowledgments
This work was supported by a Department of Energy
grant DEFC2607NT43098 through the National En
* S,ReT=4
S T, ReT = 4
* S,ReT= 0
F, ReT =0 AA
10C
 10
107
V'x)cv
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ergy Technology Laboratory. J. Mohd. Yusof. Interaction of Massive Particles with
Turbulence. PhD thesis, Cornell University, 1996.
References A. A. Zick and G. M. Homsy. Stokes flow through peri
odic arrays of spheres. J. FluidMech., 115:1326, 1982.
D. A. Drew. Mathematical modeling of twophase flow.
Annu. Rev. Fluid Mech., 15:261291, 1983.
R. Garg. Modeling and simulation of twophase flows.
PhD thesis, Iowa State University, 2009.
R. Garg, S. Tenneti, M. G. Pai, and S. Subramaniam.
Heat transfer in stokes flow past ordered and random as
semblies of monodisperse spheres. In review.
R. Garg, S. Tenneti, J. MohdYusof, and S. Subrama
niam. Direct numerical simulation of gassolids flow
based on the immersed boundary method. In S Pannala,
M. Syamlal, and T. J. O'Brien, editors, Computational
GasSolids Flows and Reacting Systems: Theory, Meth
ods and Practice. 2009. (in review).
R. J. Hill, D. L. Koch, and A. J. C. Ladd. Moderate
Reynoldsnumber flows in ordered and random arrays
of spheres. J. Fluid Mech., 448(243278), 2001.
C. M. Hrenya and J. L. Sinclair. Effects of particlephase
turbulence in gassolid flows. AIChE Journal, 43(4):
853869, 1997.
D. L. Koch. Kinetic theory for a monodisperse gassolid
suspension. Physics of Fluids, A 2(10)(17111723),
1990.
D. L. Koch and A. S. Sangani. Particle pressure and
marginal stability limits for a homogeneous monodis
perse gasfluidized bed: kinetic theory and numerical
simulations. J. Fluid Mech., 400(229263), 1999.
F. X. Li and L. S. Fan. Clean coal conversion process
progress and challenges. Fo i. i Environ. Sci., 1:248
267, 2008.
S. Tenneti and S. Subramaniam. Instantaneous particle
acceleration model for gassolid suspensions. In prepa
ration.
S. Tenneti, R. Garg, C. M. Hrenya, R. O. Fox, and
S. Subramaniam. Direct numerical simulation of gas
solid suspensions at moderate reynolds number: quan
tifying the coupling between hydrodynamic forces and
particle velocity fluctuations. In review.
M. A. van der Hoef, R. Beetstra, and J. A. M.
Kuipers. LatticeBoltzmann simulations of low
Reynoldsnumber flow past mono and bidisperse arrays
of sphere: results for the permeability and drag force. J.
Fluid Mech., 528:233254, 2005.
