J. Fluid Mech.
(2008),
vol.
597,
pp.
305–341.
c
2008 Cambridge University Press
doi:10.1017/S0022112007009834 Printed in the United Kingdom
305
Collective diffusion in sheared colloidal
suspensions
ALEXANDER M. LESHANSKY
1
,
JEFFREY F. MORRIS
2
AND
JOHN F. BRADY
3
1
Department of Chemical Engineering, Technion–IIT, Haifa, 32000, Israel
2
Benjamin Levich Institute and Department of Chemical Engineering, City College of New York,
New York, NY 10031, USA
3
Division of Chemistry & Chemical Engineering, California Institute of Technology,
Pasadena, CA 91125, USA
(Received
3 June 2007 and in revised form 18 October 2007)
Collective diffusivity in a suspension of rigid particles in steady linear viscous flows
is evaluated by investigating the dynamics of the time correlation of long-wavelength
density fluctuations. In the absence of hydrodynamic interactions between suspended
particles in a dilute suspension of identical hard spheres, closed-form asymptotic
expressions for the collective diffusivity are derived in the limits of low and high
P
́
eclet numbers, where the P
́
eclet number
Pe
=
̇
γa
2
/D
0
with
̇
γ
being the shear rate
and
D
0
=
k
B
T/
6
π
ηa
is the Stokes–Einstein diffusion coefficient of an isolated sphere
of radius
a
in a fluid of viscosity
η
. The effect of hydrodynamic interactions is studied
in the analytically tractable case of weakly sheared (
Pe
1) suspensions.
For strongly sheared suspensions, i.e. at high
Pe
, in the absence of hydrodynamics
the collective diffusivity
D
c
=6
D
s
∞
,where
D
s
∞
is the long-time self-diffusivity and both
scale as
φ
̇
γa
2
, where
φ
is the particle volume fraction. For weakly sheared suspensions
it is shown that the leading dependence of collective diffusivity on the imposed flow is
proportional to
D
0
φ
Pe
ˆ
E
,where
ˆ
E
is the rate-of-strain tensor scaled by
̇
γ
, regardless
of whether particles interact hydrodynamically. When hydrodynamic interactions
are considered, however, correlations of hydrodynamic velocity fluctuations yield
a weakly singular logarithmic dependence of the cross-gradient-diffusivity on
k
at
leading order as
ak
→
0 with
k
being the wavenumber of the density fluctuation.
The diagonal components of the collective diffusivity tensor, both with and without
hydrodynamic interactions, are of
O
(
φ
Pe
2
), quadratic in the imposed flow, and finite
at
k
=0.
At moderate particle volume fractions, 0
.
10
φ
0
.
35, Brownian Dynamics (BD)
numerical simulations in which there are no hydrodynamic interactions are performed
and the transverse collective diffusivity in simple shear flow is determined via time
evolution of the dynamic structure factor. The BD simulation results compare well
with the derived asymptotic estimates. A comparison of the high-
Pe
BD simulation
results with available experimental data on collective diffusivity in non-Brownian
sheared suspensions shows a good qualitative agreement, though hydrodynamic
interactions prove to be important at moderate concentrations.
1. Introduction
Over a century after Einstein (1905) published his seminal paper on Brownian
diffusion, the theory of particle diffusion in colloidal dispersions is still under
development. Despite the tremendous progress in understanding of diffusion inspired
306
A. M. Leshansky, J. F. Morris and J. F. Brady
by the Einstein work over the last 100 years, the theory of diffusion in non-equilibrium
flowing suspension has developed only in the last 30 years. The major obstacle
appears to be the non-trivial coupling between Brownian and hydrodynamic forces
due to shearing motion. The influence of these factors is not simply additive, and
shear-induced enhancement of particle diffusivity or, more generally, dispersion is
anticipated as in the analogous problem of passive scalar transport (Taylor 1953;
Novikov 1958; Elrick 1962; Batchelor 1979).
As is the case at equilibrium, in suspensions undergoing steady shear two distinct
diffusive processes can be identified: the particle
self-diffusivity
, related to the rate
of change of the variance in the position of a tagged particle, and the
collective
diffusivity
(or gradient diffusivity) corresponding to a rate at which macroscopic
spatial inhomogeneities in particle distribution relax. Also, as in a macroscopically
quiescent suspension, the particles are expected to behave diffusively on two separate
time scales characterized by short- and long-time diffusion coefficients (see Pusey
1991). The latter results from the hindering effect of the direct interactions (e.g.
hydrodynamic or steric) between suspended particles as they move through the
suspension. In sheared suspensions the spatial anisotropy of the imposed flow results
in diffusion coefficients that are anisotropic and the variance in position of a particle
may grow faster than linearly in time owing to the Taylor dispersion driven by the
position-dependent imposed flow.
The early studies of diffusion in sheared systems (e.g. san Miguel & Sancho 1979)
used a transformation to a coordinate system co-moving with the shearing motion
to remove the linear flow from the governing equations. In Frankel & Brenner
(1991), a generalized Taylor dispersion theory was developed to study transport of an
isolated particle in linear flows. An appropriate transformation of the time coordinate
was applied to remove the position-dependent advection of the bulk flow, but the
complexity of the analysis may not allow it to be applied to study of the transport of
interacting particles or collective diffusion.
At the extreme of strong shearing, there have been a number of studies of shear-
induced or
hydrodynamic
diffusion (Eckstein, Bailey & Shapiro 1977; Leighton &
Acrivos 1987; Breedveld
et al.
1998). Trajectory analysis was applied in Acrivos
et al.
(1992) to calculate the self-diffusivity, scaling as
φ
̇
γa
2
in the flow direction,
with
̇
γ
being the shear rate and
a
the sphere radius. The
O
(
φ
2
) self- and collective-
diffusion coefficients in the velocity-gradient direction were determined by Wang,
Mauri & Acrivos (1996 and 1998, respectively) from trajectory analysis including
three-particle interactions in a monolayer of spheres.
†
In da Cunha & Hinch (1996),
the diffusion coefficients were derived from the pair-trajectory calculations of spheres
with slight surface roughness. A similar analysis was applied in Wang & Mauri (1999)
to calculate a collective-diffusion coefficient from analysis of pairwise interactions at
O
(
φ
). Nevertheless, applicability of trajectory analysis will probably remain limited
to studies of diffusion in the absence of Brownian transport, because it is not clear
how to introduce stochastic forces into the formulation.
An alternative approach was recently proposed by Mauri (2003) where homogeni-
zation techniques were applied to derive a constitutive relation for the coarse-grained
†
The hydrodynamic encounter between two perfectly smooth spheres in Stokes flow does not
lead to a net lateral displacement of either sphere due to the linearity and time-reversibility of
the equations of motion. Thus, at the pair level a non-hydrodynamic microscopic mechanism (e.g.
particle roughness, repulsive forces, Brownian forces, etc.) must be included to provide microscopic
irreversibility that would lead to net displacement and diffusion upon many successive encounters.
Collective diffusion in sheared colloidal suspensions
307
averaged volumetric flux of non-Brownian particles in an inhomogeneous and non-
stationary viscous flow. Though the general constitutive relation for the collective
diffusivity in terms of the velocity cross-correlation was derived, no comparison with
existing experimental results or theoretical predictions have been performed. As in the
trajectory technique discussed above, generalization of the method of Mauri (2003) to
allow consideration of combined hydrodynamic and Brownian transport is missing.
Morris & Brady (1996) presented a methodology stemming from the experimental
technique of dynamic light scattering (Berne & Pecora 1976) to calculate particle
self-diffusivity in sheared suspensions, and applied it to the case of weak shearing
(small P
́
eclet number,
Pe
). Although successful theories referring to dynamic light
scattering had been proposed to study diffusion in quiescent colloidal dispersions
of hydrodynamically interacting particles (Russel & Glendinning 1981; Rallison &
Hinch 1986; Brady 1994), the generalization of the method by Morris & Brady
allowed investigation of particle transport in non-equilibrium sheared suspensions.
The central idea of the equilibrium theories is based upon the relation between
diffusivity (self- and collective) and the rate of temporal decay of the autocorrelation
in number density fluctuations (decorrelation in scattered light is due to uncorrelated
particle motions). This idea has been shown to have the same role in the theory of self-
diffusivity in sheared suspensions (Morris & Brady 1996). Under shear the variance
in particle position is not readily related to diffusion as in a quiescent suspensions
owing to the Taylor dispersion driven by the position-dependent imposed flow, while
the Fourier-transform method was shown to be useful in identifying particle self-
diffusivity by analogy with the treatment of passive scalar transport (see Batchelor
1979). The short-time self-diffusivity (the product of the average mobility of a particle
and
k
B
T
) and the long-time self-diffusivity (due to the correlation of the flux of
the tagged particles with the microstructural perturbation caused by their motion)
were determined at
O
(
φ
); note that the microstructure associated with the short-time
self-diffusivity is that of the steadily sheared suspension at the shear rate of interest.
In a later paper, Brady & Morris (1997) demonstrated that the same methodology
is applicable to study the effect of strong shear on self-diffusivity and to predict
the
O
(
φ
) coefficient of the
̇
γa
2
self-diffusivity in a general linear flow. The obvious
advantage of the Fourier-transform method is that it places the analysis of quiescent
and flowing suspensions on the same footing. It was suggested in Morris & Brady
(1996) that the approach can be applied to calculation of the collective diffusivity in
sheared colloidal suspensions, as done for the quiescent suspension by Rallison &
Hinch (1986), but no such study has been attempted prior to the present work.
Although self-diffusion coefficients can be calculated using appropriate kinematic
descriptions, such as the mean-square displacement or velocity autocorrelation
function of a particle, and these apply even in highly non-equilibrium systems
(Marchioro & Acrivos 2001; Drazer
et al.
2002; Sierou & Brady 2004), evaluation
of the gradient diffusivity is less obvious. For macroscopically quiescent suspensions
there is a clear physical description of the short-time collective diffusivity through the
wave-vector-dependent diffusion coefficient
D
(
k
), which depends on both structural
and hydrodynamic effects (Russel & Glendinning 1981; Pusey 1991; Segr
́
e, Behrend &
Pusey 1995),
D
(
k
)
D
0
=
H
(
k
)
S
(
k
)
,
(1.1)
where
k
is the magnitude of the wave-vector of the density fluctuation,
H
(
k
) represents
purely hydrodynamic effects, and structural effects are contained in the static structure
308
A. M. Leshansky, J. F. Morris and J. F. Brady
factor
S
(
k
). When
ak
1, we are in the long-wavelength limit where the scale of
the density variation is large relative to the particle size. The relaxation of long-
wavelength density fluctuations yields a particle flux which is the same as if a small
constant density gradient persisted everywhere (Rallison & Hinch 1986). Therefore,
as
ka
→
0, one expects that
D
(
k
) will asymptote to the constant collective or down-
gradient diffusivity,
D
c
. Also in the limit of vanishing
k
, the magnitude of the density
fluctuations of a system, as characterized by the static structure factor
S
(0), is related to
its osmotic or isothermal compressibility
S
(0) =
k
B
T
(
∂n/∂Π
)
T
(Berne & Pecora 1976;
Pusey 1991). Therefore,
S
(0)
−
1
can be interpreted as a thermodynamic force associated
with the density gradient which drives the diffusion. The hydrodynamic factor,
H
(
k
),
can be expressed through an integral over the velocity cross-correlation (Pusey 1991;
Segr
́
e
et al.
1995). For a dilute colloidal suspension of hard spheres,
S
(0)
∼
1
−
8
φ
due to excluded volume effects, while the collective mobility is given by the mean
sedimentation velocity, relative to zero volume-flux axes
H
(0)
∼
(6
π
ηa
)
−
1
(1
−
6
.
55
φ
),
leading to the classical result
D
c
=
k
B
T
H
(0)
/S
(0)
∼
D
0
(1 + 1
.
45
φ
) (Batchelor 1976);
more accurate hydrodynamic data yield a slightly different value of 1
.
47
φ
for the last
term in parentheses (Rallison & Hinch 1986).
In Marchioro & Acrivos (2001), a novel numerical scheme for calculating the
gradient diffusivity in non-Brownian suspensions, which should also be applicable
to the case with Brownian motion, using Stokesian Dynamics simulations has
been proposed. The essence of the method is that the suspension microstructure
relaxes at a rate proportional to the gradient-diffusion coefficient. Thus, if an initial
microstructure is sufficiently different from that of a suspension achieved after
extended steady shearing, probing the rate of relaxation from the initial toward
the steady microstructure provides a means of determining the gradient diffusivity.
A different approach was recently proposed in Leshansky & Brady (2005), where
it was demonstrated that representation (1.1) holds for the tensorial wave-vector-
dependent coefficient
D
of suspensions in arbitrary steady linear flow. Again the
method was applied to describe the relaxation of the long-wavelength density
fluctuations,
ak
→
0. The work showed that the tensorial hydrodynamic factor
H
can be expressed through the averaged integral of the velocity cross-correlation
function and
S
(0) can be interpreted as that corresponding to the non-equilibrium
shear-induced suspension microstructure. This approach is conceptually different from
that of Marchioro & Acrivos (2001) as it concerns the decay of the fluctuations in
number density rather than relaxation of the macroscopic density itself. For the self-
diffusivity this approach yields a correct description of the self-diffusivity in terms of
the integral over the velocity autocorrelation derived earlier from purely kinematic
arguments by Marchioro & Acrivos (2001) and Sierou & Brady (2004). The results
for the shear-induced collective-diffusion coefficients (in the absence of Brownian
transport) in the velocity-gradient and vorticity directions obtained in Leshansky &
Brady from accelerated Stokesian Dynamics simulations (Sierou & Brady 2001)
show very good agreement with available experimental data. The derivation is
independent of the details of the microscale dynamics and thus should hold quite
generally.
The present paper is a natural continuation of the work in Morris & Brady (1996),
Brady & Morris (1997), and Leshansky & Brady (2005). We address the problem of
calculating the collective diffusivity in a suspension of hard spheres in an arbitrary
steady linear flow at small Reynolds number. The light-scattering-based method of
Morris & Brady is applied to investigation of the rate of decay of the density
fluctuations.
Collective diffusion in sheared colloidal suspensions
309
The paper is organized as follows. In
§
2 we show how the collective diffusivity can
be identified via the dynamic structure factor approach and derive the closed-form
asymptotic expression for the tensorial diffusion coefficient in the dilute limit for
weakly and strongly sheared suspensions. In
§
3 we briefly describe the details of the
Brownian Dynamics simulations and present results for the transverse diffusivities
along both principal directions of a simple shear flow (velocity gradient and vorticity).
In
§
4 a brief summary and concluding remarks are provided. The Appendix concerns
the long-wavelength microstructural deformation of a suspension in simple shear flow
in terms of the non-equilibrium static structure factor
S
(
k
=0).
2. Theoretical development
2.1.
Dynamic structure factor approach
To understand how we can define the particle diffusivity in linear flow, consider the
following conservation equation for the local particle number density,
n
(
x
,t
):
∂n
∂t
+
̇
Γ
·
x
·
∇
n
+
U
·
∇
n
=
−∇
·
j
,
(2.1)
where
̇
Γ
is the constant velocity-gradient tensor,
U
is the bulk average velocity
measured at an arbitrary field point,
x
0
, from which the bulk shear velocity is
referenced and
j
is the diffusive flux of particles. Neglecting memory effects for time
scales over which the diffusion is stationary (either short or long time scale), the flux
should be expressible as a generalized Fick’s law:
j
=
−
D
(
x
−
x
)
·
∇
n
(
x
,t
)d
x
,
(2.2)
where the non-local kernel is identified as the
effective diffusivity
with dimensions
(
LT
)
−
1
. With the definition of the Fourier transform (
F
)pair
h
(
k
)=
h
(
x
)e
−
i
k
·
x
d
x
,
h
(
x
)= (2
π
)
−
3
h
(
k
)e
i
k
·
x
d
k
,
where
k
is the wave-vector in Fourier space, and using the identity
G
(
x
−
y
)
·
W
(
y
)d
y
=(2
π
)
−
3
W
(
k
)
·
G
(
k
)e
i
k
·
x
d
k
,
the expression for flux (2.2) can be rewritten as
j
=
−
(2
π
)
−
3
D
(
k
)
·
i
k
n
(
k
,t
)e
i
k
·
x
d
k
.
(2.3)
Applying the Fourier transform to (2.3) and using the identity
F{
e
i
x
·
(
k
−
k
)
}
=
(2
π
)
3
δ
(
k
−
k
), we arrive at the Fourier-space expression for
j
:
j
=
−
D
(
k
)
·
i
k
n
(
k
,t
)
,
(2.4)
where the
k
-dependent diffusivity
D
(
k
) has the dimensions
L
2
/T
. The Fourier
transform applied to the left-hand side of (2.1) combined with (2.4) converts (2.1) to
∂
n
∂t
−
k
·
̇
Γ
·
∇
k
n
−
i
k
·
U
n
=
−
k
·
D
(
k
)
·
k
n.
(2.5)
This equation is analogous to that which describes the evolution of a passive scalar
field in a linear flow, in which case there would be a constant isotropic diffusivity
310
A. M. Leshansky, J. F. Morris and J. F. Brady
D
=
D
I
. Note that the spreading or dispersion of the passive scalar due to the shear
flow, characterized by the second moment
xx
n
(
x
,t
)d
x
, may grow faster than
linearly in time. This is the well-known Taylor dispersion (Taylor 1953; Novikov
1958) and is accounted for by the
k
·
̇
Γ
·
∇
k
n
term on the left-hand side of (2.5). The
general solution of (2.5) for an arbitrary linear flow can be found for a constant
k
-independent tensorial coefficient
D
(Batchelor 1979).
In a steadily sheared suspension the average density fluctuation is zero and
therefore we consider the two-point time autocorrelation of the local particle
density,
F
(
k
,t
)=
n
(
k
,t
)
n
∗
(
k
,
0)
,where
∗
indicates a complex conjugate and
angle brackets denote an appropriate statistical averaging (either ensemble or time
averaging). Multiplying both sides of (2.5) by
n
∗
(
k
,
0) the evolution equation for
F
becomes
∂F
∂t
−
k
·
̇
Γ
·
∇
k
F
−
i
k
·
U
F
+
n
k
·
̇
Γ
·
∇
k
n
∗
=
−
k
·D
(
k
)
·
k
F,
(2.6)
where we have removed the hat over
D
for simplicity. The above equations for
n
and
F
have the expected form for a diffusive process in a linear flow, and can be used
to identify the proper expressions for the particle diffusivities either from statistical
mechanical theory or from dynamic simulations of multiparticle systems. Note that
the last term on the left-hand side of (2.6) may lead to a modification of the particle
dispersion in comparison with the dispersion of a passive scalar as described by (2.5),
where this term is missing. In Morris & Brady (1996), the working hypothesis was
that the relaxation of the scalar field
F
is governed by (2.5) simply by analogy and,
therefore, this term was omitted. (Retaining the last term on the left-hand side of
(2.6) does not alter the outcome of the analogous analysis for the self-diffusivity in
Morris & Brady (1996) as the diffusive variation of
F
is
O
(
k
2
), while linear flow
causes variation which is independent of
k
.)
As in macroscopically quiescent suspensions (e.g. Russel & Glendinning 1981) the
limiting value of
D
as
k
→
0 is identified as the collective-diffusion coefficient
D
c
.
To make analytical progress, the particle number density at any point in the
suspension
x
is defined in terms of distributions as
n
(
x
,t
)=
N
α
=1
δ
(
x
−
x
α
(
t
))
,
(2.7)
where
δ
is the Dirac delta function, and
x
α
(
t
) is the position of particle
α
at time
t
.
The spatial Fourier transform of
n
(
x
)is
ˆ
n
(
k
,t
)=
e
i
k
·
x
N
α
=1
δ
(
x
−
x
α
(
t
)) d
x
=
N
α
=1
e
i
k
·
x
α
.
(2.8)
The density autocorrelation function, known as the
dynamic structure factor
or as
the intermediate scattering function (primarily in colloidal and polymer literature), is
given by
F
(
k
,t
)=
1
N
ˆ
n
(
k
,t
)
ˆ
n
∗
(
k
,
0)
=
1
N
N
α,β
e
i
k
·
(
x
α
(
t
)
−
x
β
(0))
,
(2.9)
where the angular brackets
denote an ensemble average over the joint probability
P
N
(
x
N
(
t
)
,t
;
x
N
(0)
,
0) of particles being at
x
N
(0) at time 0 and at
x
N
(
t
) at time
t
.The
constant background density level of
n
in (2.9) is irrelevant and thus
F
describes the
Collective diffusion in sheared colloidal suspensions
311
autocorrelation of density
fluctuations
(Berne & Pecora 1976). Integrating over the
joint probability density (2.9) yields, for statistically indistinguishable particles,
F
(
k
,t
)=
1
N
1
N
!
N
α,β
e
i
k
·
(
x
α
(
t
)
−
x
β
(0))
P
N
(
x
N
(
t
)
,
x
N
(0)) d
x
N
d
x
N
(0)
=
1
N
!
(e
i
k
·
(
x
1
(
t
)
−
x
1
(0))
+(
N
−
1)e
i
k
·
(
x
2
(
t
)
−
x
1
(0))
×
P
N
d
x
N
(
t
)d
x
N
(0)
=
F
s
(
k
,t
)+
N
−
1
N
!
e
i
k
·
(
x
2
−
x
1
(0))
P
N
d
x
N
(
t
)d
x
N
(0)
,
(2.10)
where we have introduced the self-dynamic structure factor (or self-intermediate
scattering function)
F
s
(
k
,t
).
Following Rallison & Hinch (1986), we further define
ˆ
P
N
(
x
N
,t
;
k
)=
P
N
(
x
N
|
x
N
(0))
P
0
N
(
x
N
(0))e
−
i
k
·
x
1
(0)
d
x
N
(0)
,
(2.11)
where the transitional probability
P
N
(
x
N
|
x
N
(0)) is governed by a Smoluchowski
equation (Morris & Brady 1996)
∂P
N
∂t
+
N
α
=1
∇
α
·
j
α
=0
,
(2.12)
with initial value
δ
(
x
N
−
x
N
(0)) at
t
= 0. In (2.12),
j
α
is the probability flux associated
with particle
α
, given by
j
α
=
U
α
P
N
−
N
β
=1
D
αβ
P
N
·
∇
β
(
ln
P
N
+
V
)
,
(2.13)
where
D
αβ
=
k
B
T
M
αβ
with
M
αβ
the hydrodynamic mobility of particle
α
due to a
force on particle
β
,
V
(
r
N
) is the interparticle potential (independent of absolute
position) and
U
α
is the hydrodynamic velocity that particle
α
would have due to
external forces and bulk imposed flow, in the absence of Brownian motion and an
interparticle potential.
Integrating in (2.10) over the initial coordinates
x
N
(0) using (2.11) yields the
expressions for
F
and
F
s
, respectively,
F
=
1
N
!
[1 + (
N
−
1)e
i
k
·
r
2
]
ˆ
P
N
e
i
k
·
x
1
d
x
1
d
r
N
,F
s
=
1
N
!
ˆ
P
N
e
i
k
·
x
1
d
x
1
d
r
N
,
(2.14)
where the new coordinates
x
1
and
r
N
≡
(
r
2
,...,
r
N
)
,
have been introduced where
x
1
=
x
1
and
r
α
=
x
α
−
x
1
,α
2, with the original variables
denoted by primes. The formulation for the collective diffusivity can be expressed in
terms of the configuration-dependent perturbation function
f
N
(
r
N
,t
;
k
),
ˆ
P
N
≡
ˆ
P
1
P
0
N
−
1
|
1
[1 +
f
N
]
≡
ˆ
P
1
Q,
(2.15)
where
ˆ
P
1
(
x
1
,t
;
k
)=
1
(
N
−
1)!
ˆ
P
N
d
r
N
312
A. M. Leshansky, J. F. Morris and J. F. Brady
is the probability density for particle 1 chosen at random (Morris & Brady 1996) and
the function
P
0
N
−
1
|
1
(
r
N
)isthe
steady non-equilibrium conditional
probability of finding
N
−
1 particles given that there is particle 1 at
x
1
.
Using the Smoluchowski equation (2.12), governing the conservation of
ˆ
P
N
(
x
1
,
r
N
;
k
), and integrating by parts, we arrive at
̇
F
=
1
N
!
[1 + (
N
−
1)e
i
k
·
r
2
]e
i
k
·
x
1
∂
ˆ
P
N
∂t
d
x
1
d
r
N
=i
k
·
1
N
!
[
ˆ
j
1
+(
N
−
1)
ˆ
j
2
e
i
k
·
r
2
]e
i
k
·
x
1
d
x
1
d
r
N
,
(2.16)
where
ˆ
j
α
is the probability flux associated with particle
α
given by (2.13), but now in
the new coordinates (
x
1
,
r
N
), instantaneously centred at
x
1
,
ˆ
j
α
=
U
α
ˆ
P
N
−
D
α
1
·
∇
1
ˆ
P
N
−
N
β
=2
(
D
αβ
−
D
α
1
)
ˆ
P
N
·
∇
β
(ln
ˆ
P
N
+
V
)
,
1
α
N.
(2.17)
For an arbitrary linear flow
U
α
may be written as
U
α
=
̇
Γ
·
(
r
α
+
x
1
)+
U
∗
+
U
α
(
r
N
)
,α
2
,
(2.18)
where
̇
Γ
is the constant traceless velocity-gradient tensor,
U
∗
=
U
∞
(
x
0
)
−
̇
Γ
·
x
0
with
U
∞
being the bulk average velocity measured at an arbitrary field point,
x
0
,and
U
α
is the configuration-dependent velocity fluctuation. Equation (2.18) holds for particle
1 upon substitution
r
1
= 0. Substituting (2.17) into (2.16) and taking into account that
F
s
=
N
−
1
ˆ
P
1
e
i
k
·
x
1
d
x
1
we integrate out
x
1
to obtain
̇
F
=
k
·
̇
Γ
·
∇
k
F
−
1
N
n
k
·
̇
Γ
·∇
k
n
∗
+i
k
·
U
∗
F
−
F
s
kk
:
1
(
N
−
1)!
[
D
11
+(
N
−
1)
D
21
e
i
k
·
r
2
]
Q
d
r
N
+
F
s
i
k
·
1
(
N
−
1)!
[
U
1
+(
N
−
1)
U
2
e
i
k
·
r
2
]
Q
d
r
N
−
F
s
i
k
·
1
(
N
−
2)!
[
D
12
−
D
11
+(
D
22
−
D
21
)e
i
k
·
r
2
]
Q
·
∇
2
[ln
Q
+
V
]d
r
N
−
F
s
i
k
·
1
(
N
−
3)!
(
D
23
−
D
21
)e
i
k
·
r
2
Q
·
∇
3
[ln
Q
+
V
]d
r
N
.
(2.19)
The first two terms on the right-hand side of (2.19) result from the undisturbed linear
velocity in (2.18):
i
k
·
1
N
!
[
̇
Γ
·
x
1
+(
N
−
1)
̇
Γ
·
(
r
2
+
x
1
)e
i
k
·
r
2
]e
i
k
·
x
1
ˆ
P
N
d
x
1
d
r
N
=
1
N
n
∗
(
k
,
0)
k
·
̇
Γ
·∇
k
n
(
k
,t
)
=
k
·
̇
Γ
·
∇
k
F
−
1
N
n
k
·
̇
Γ
·∇
k
n
∗
.
(2.20)
We can simplify (2.19) further by noting that the total velocity disturbance must
vanish, i.e.
[
U
1
+(
N
−
1)
U
2
]
Q
d
r
N
=0
,
Collective diffusion in sheared colloidal suspensions
313
so that the second integral on the right-hand side in (2.19) can be rewritten
F
s
i
k
·
1
(
N
−
2)!
U
2
(e
i
k
·
r
2
−
1)
Q
d
r
N
.
(2.21)
Next we note that the particle mobility
D
11
is given by the isolated-particle value
D
0
I
plus the contribution from
N
−
1 other particles
†
,i.e.
D
11
=
D
0
I
+(
N
−
1)
D
11
.
Further, for identical particles
D
12
−
D
11
+
D
22
−
D
21
≡
0
.
Thus we have
̇
F
=
k
·
̇
Γ
·
∇
k
F
−
1
N
n
k
·
̇
Γ
·∇
k
n
∗
+i
k
·
U
∗
F
−
F
s
Dk
2
−
F
s
kk
:
1
(
N
−
2)!
(
D
11
+
D
21
e
i
k
·
r
2
)
Q
d
r
N
+
F
s
i
k
·
1
(
N
−
2)!
U
2
(e
i
k
·
r
2
−
1)
Q
d
r
N
−
F
s
i
k
·
1
(
N
−
2)!
[(
D
22
−
D
21
)(e
i
k
·
r
2
−
1)]
Q
·
∇
2
[ln
Q
+
V
]d
r
N
−
F
s
i
k
·
1
(
N
−
3)!
(
D
23
−
D
21
)e
i
k
·
r
2
Q
·
∇
3
[ln
Q
+
V
]d
r
N
.
(2.22)
Substituting (2.15) into
F
as given by (2.14) and integrating over
x
1
, it appears that
F
is linear in
F
s
:
F
(
k
,t
)=
F
s
1
(
N
−
1)!
[1 + (
N
−
1)e
i
k
·
r
2
]
Q
d
r
N
.
(2.23)
Here we have made use of the fact that
F
s
(0
,t
) = 1, so that
F
s
(
k
,t
)
δ
(
k
)=
δ
(
k
). Since
the average perturbation must be zero,
P
0
(
N
−
1)
|
1
f
N
d
x
N
=0
,
it follows from (2.23) that in the low-
k
limit
F
(
k
,t
) asymptotes to
S
(
k
),
F
(
k
,t
)=1+
n
(
g
(
r
)
−
1)e
i
k
·
r
d
r
+
O
(
k
2
)=
S
(
k
)+
O
(
k
2
)
,
(2.24)
where
g
(
r
) is the steady pair-distribution function defined by
P
0
1
|
1
=
ng
(
r
)and
S
(
k
)is
the
non-equilibrium static structure factor
. In the derivation of (2.24) we have ignored
the delta-function term
Nδ
(
k
); this Fourier transform of the large scattering volume
V
is irrelevant to our subsequent analysis (see Rallison & Hinch 1986).
Quantities are scaled as
r
∼
a,
U
∼
̇
γa,
k
∼
a
−
1
,
D
∼
D
0
and
t
∼
a
2
/D
0
, and the
P
́
eclet number is defined as
Pe
=
̇
γa
2
/D
0
,where
̇
γ
≡
(
̇
Γ
:
̇
Γ
)
1
/
2
.Wedonot
employ alternative symbols to denote dimensionless quantities and all quantities
are dimensionless unless stated otherwise in the following. By comparing with (2.6)
the short-time effective transport properties may be determined by considering their
†
It will be shown later that velocity fluctuations cannot be added in the same pairwise manner,
i.e.
U
i
=
N
j
=
i
U
ij
, since it would lead to divergence of the resulting integral.
314
A. M. Leshansky, J. F. Morris and J. F. Brady
contribution to the initial value of the logarithmic derivative
∂
t
log
F
,
∂
t
log
F
(
k
,
0) =
Pe
k
·
̇
Γ
·
∇
k
ln
F
−
Pe
(
NS
0
)
−
1
n
k
·
̇
Γ
·∇
k
n
∗
+
Pe
i
k
·
U
∗
−
1
S
0
k
2
−
1
S
0
kk
:
1
(
N
−
2)!
(
D
11
+
D
21
e
i
k
·
r
2
)
P
0
d
r
N
+
Pe
S
0
i
k
·
1
(
N
−
2)!
U
2
(e
i
k
·
r
2
−
1)
P
0
d
r
N
−
1
S
0
i
k
·
1
(
N
−
2)!
[(
D
22
−
D
21
)(e
i
k
·
r
2
−
1)]
P
0
·
∇
2
[ln
P
0
+
V
]d
r
N
−
1
S
0
i
k
·
1
(
N
−
3)!
(
D
23
−
D
21
)e
i
k
·
r
2
P
0
·
∇
3
[ln
P
0
+
V
]d
r
N
,
(2.25)
where
P
0
≡
P
0
N
−
1
|
1
to simplify notation and
S
0
≡
S
(0) = 1 +
n
(
g
(
r
)
−
1) d
r
. The short-
time or initial collective velocity, giving the
O
(
k
) term in (2.25), is
U
c
0
≡
U
∗
−
1
S
0
1
(
N
−
3)!
(
D
23
−
D
21
)
P
0
·
∇
3
[ln
P
0
+
V
]d
r
N
.
(2.26)
If the initial distribution is the equilibrium one, then the integral in (2.26) is zero.
This integral is likely to be zero for most flows, based on symmetry arguments.
Expanding the exponential e
i
k
·
r
2
and identifying the terms quadratic in
k
, the initial
collective diffusivity is found from (2.25) at
k
→
0as
D
c
0
≈
1
S
0
I
+
1
S
0
1
(
N
−
2)!
[
D
11
+
D
21
]
P
0
d
r
N
−
F
kk
−
1
S
0
1
(
N
−
2)!
r
2
(
D
22
−
D
21
)
P
0
·
∇
2
[ln
P
0
+
V
]d
r
N
−
1
S
0
1
(
N
−
3)!
r
2
(
D
23
−
D
21
)
P
0
·
∇
3
[ln
P
0
+
V
]d
r
N
,
(2.27)
where
F
kk
is the
k
2
tensorial coefficient of the small-
k
expansion of the velocity
fluctuation integral – the term on the right-hand side of (2.25) containing
U
2
.In
deriving (2.27) we have expanded the factor e
i
k
·
r
2
for small
k
(except for the velocity
fluctuation integral) assuming that all integrals over
r
N
are convergent.
At long times, the perturbation
f
N
is
O
(
k
) due to the long-wave nature of the
response and can be written as
f
N
=i
k
·
b
N
,where
b
N
(
r
N
,t
) is a vector function
independent of
k
(Morris & Brady 1996; Brady & Morris 1997). Substituting this
ansatz into (2.22) we find that the long-time collective velocity is unchanged,
U
c
∞
=
U
c
0
,
(2.28)
while the long-time collective diffusivity differs from its initial value only in the terms
involving third-particle effects,
D
c
∞
=
D
c
0
−
1
S
0
1
(
N
−
3)!
b
N
(
D
23
−
D
21
)
·
∇
3
[ln
P
0
+
V
]d
r
N
−
1
S
0
1
(
N
−
3)!
(
D
23
−
D
21
)
P
0
·
∇
3
b
N
d
r
N
.
(2.29)
We conclude that when hydrodynamic interactions are neglected, there is no change
in the collective diffusivity with time for an arbitrary P
́
eclet number and particle
Collective diffusion in sheared colloidal suspensions
315
concentration. In the following, we therefore omit the subscripts ‘0’ and ‘
∞
’of
D
c
for
simplicity.
2.2.
No hydrodynamic interactions
In this section, we develop results given in general form above for the case in which
the particles have no hydrodynamic interactions. Results are developed completely
(to numerical values with appropriate tensor form) for the low- and high-
Pe
limits.
Following from the point just made at the end of
§
2.1, we see that in the absence
of hydrodynamic interactions between particles, the long-time collective diffusivity
equals its initial value, and (2.27) reduces to
D
c
=
1
S
0
I
−
1
(
N
−
2)!
r
2
P
0
∇
2
[ln
P
0
+
V
]d
r
N
.
(2.30)
The influence of flow is found in the
Pe
-dependent structure, i.e. via both
P
0
and
S
0
.
Referring to (2.27), we see that third-particle effects enter through the final integral
which vanishes when hydrodynamic interactions are neglected, and thus (2.30) is not
simply a pair-level result. We can further write
1
(
N
−
2)!
r
2
P
0
∇
2
[ln
P
0
+
V
]d
r
N
=
n
r
g
(
r
)
∇
[ln
g
(
r
)+
V
0
2
]d
r
,
where
V
0
2
=
1
(
N
−
2)!
VP
0
N
−
2
|
2
(
r
3
,...,
r
N
|
r
2
)d
r
3
···
d
r
N
is the potential of mean force averaged over the conditional probability
P
0
N
−
2
|
2
of
finding
N
−
2 particles at
r
3
,...,
r
N
, given that there is a particle at
r
2
. Hence, the
collective diffusivity becomes
D
c
=
1
S
0
I
−
φ
3
4
π
r
g
(
r
)
∇
[ln
g
(
r
)+
V
0
2
]d
r
,
(2.31)
which is valid for arbitrary
Pe
and particle volume fraction
φ
=
4
3
π
a
3
n
.At
thermodynamic equilibrium the particles are Boltzmann distributed,
g
eq
= exp(
−
V
0
2
),
and the equilibrium collective diffusivity is isotropic and is (in the dimensional form)
D
c
eq
=
D
0
S
eq
(0)
=(6
π
ηa
)
−
1
∂Π
∂n
,
(2.32)
as it should be for non-hydrodynamically interacting particles. To make progress in
the general case, i.e. out of equilibrium, (2.31) can be simplified in the dilute limit,
where for hard spheres
g
eq
(
r
)=1 at
|
r
|
>
2and
V
0
2
= 0. After integration by parts,
(2.31) becomes
D
c
=
1
S
0
I
+
I
φ
3
4
π
r>
2
(
g
(
r
)
−
1
)
d
r
+
φ
6
π
r
=2
nn
(
g
(2)
−
1
)
d
Ω
,
(2.33)
where
n
=
−
r
/
2 is the outer unit normal vector to the surface of contact,
r
=2, and
g
(2) denotes the contact value of the pair-distribution function (which has angular
dependence on the solid angle
Ω
). We have used
∇
g
(
r
)=
∇
(
g
(
r
)
−
1) in order to
neglect the integral at infinity. The first term plus first integral on the right of (2.33)
can be written in terms of the static structure factor deviation from its equilibrium
316
A. M. Leshansky, J. F. Morris and J. F. Brady
value at
k
=0 as
D
c
=
1
S
0
I
(
1+
S
0
)
+
φ
6
π
nn
(
g
(2)
−
1)d
Ω
,
(2.34)
where
S
0
=
S
0
−
S
eq
0
and
S
eq
0
∼
1
−
8
φ
. Note that at equilibrium,
g
(2) = 1,
S
0
=0
and thus (2.32) is recovered.
High P
́
eclet number
For linear flows at large P
́
eclet number, the asymptotic solution for the pair-
distribution function
g
(
r
) obtained in the radial balance approximation (Brady &
Morris 1997) gives contact value of
g
(2)
∼−
2
3
Pe
γ
r
+
O
(1)
,
for
γ
r
=
n
·
ˆ
E
·
n
<
0
,
(2.35)
where
Pe
1and
ˆ
E
ij
=(
̇
Γ
ij
+
̇
Γ
ji
)
/
2
̇
γ
is the dimensionless rate-of-strain tensor.
Length is scaled with the sphere radius
a
. It follows from (2.34) that for a dilute
suspension
D
c
=
1
S
0
(1 +
S
0
)
I
−
φ
Pe
4
π
γ
r
<
0
γ
r
nn
d
Ω
+
O
(1)
.
(2.36)
Although the shear-induced diffusion enhancement given by the integral in (2.36)
can be calculated for an arbitrary linear flow, the asymptotic evaluation of
S
(0)
in the large-
Pe
limit is a formidable task even without hydrodynamic interactions.
While the radial balance approximation can be exploited to determine the leading-
order solution for
g
(
r
) in the compressional quadrants, it is not applicable in the
extensional quadrants (
γ
r
>
0), where low-
g
wake-like regions extend from the surface
of contact to some large distance
r
∼
L
(
Pe
) downstream where they are smeared by
a weak transverse diffusion. In the compressional quadrants
g
(
r
)
∼
Pe
in the vicinity
of the surface of contact (2.35) while the thickness of the radial boundary layer is
∼
Pe
−
1
. Therefore, the boundary layer at small
φ
contributes the value of
O
(1) to the
volume integral of (
g
−
1) and its positive contribution to the
S
0
is estimated as
O
(
φ
) at most. On the other hand, the long low-
g
wakes in the extensional quadrants
may diminish the value of
S
0
considerably due to the flow-induced ‘excluded volume’
effect. The microstructural deformation at large
Pe
was determined numerically by
expanding
g
(
r
) into spherical harmonics and integrating (
g
−
1) in the entire domain.
The negative (i.e. ‘excluded volume’) contribution of wake-like regions to
S
0
in the
simple shear flow, was estimated to be of
O
(
φ
√
Pe
) (see the Appendix for details).
Since for the dilute suspension
S
0
=1
−
8
φ
+
S
0
, one can re-write the first term in
(2.36) as
1
S
0
(1 +
S
0
)=1+
8
φ
1+
S
0
;
the effect of the suspension microstructure deformation on the collective diffusivity is
always small and the second term in (2.36), which scales linearly with
Pe
, yields the
leading-order contribution to
D
c
.
Previous calculations by Brady & Morris (1997) showed that the dimensional
long-time self-diffusivity for the same problem is
D
s
∞
=
D
0
I
−
̇
γa
2
φ
2
3
π
γ
r
<
0
γ
r
nn
d
Ω.
(2.37)
Collective diffusion in sheared colloidal suspensions
317
Comparing only the leading high-P
́
eclet-number terms gives
D
c
=6
D
s
∞
=
−
a
2
η
2
9
∂
Σ
∂φ
,
where
Σ
is the contribution to the bulk stress from contact integrals in the com-
pressional quadrants,
γ
r
<
0:
Σ
=
η
̇
γφ
2
9
π
γ
r
<
0
γ
r
nn
d
Ω.
For example, in the simple shear flow,
̇
Γ
ij
=
δ
i
1
δ
2
j
, we have in the dimensional form
D
c
11
=
D
c
22
=2
D
c
33
=
32
15
π
̇
γa
2
φ.
(2.38)
Low P
́
eclet number
At low P
́
eclet number the pair-distribution function can be written
g
(
r
)=
g
eq
(
r
)[1+
p
(
r
)], where
p
is the entire
Pe
-dependent perturbation to the equilibrium micro-
structure. For a dilute suspension,
g
eq
= 1, and we have, in place of (2.34),
D
c
=
1
S
0
(1 +
S
0
)
I
+
φ
6
π
nn
p
(2)d
Ω
.
(2.39)
Here the shear-induced deformation of the static structure factor from its equilibrium
value can be written as
S
0
−
S
eq
0
=
S
0
=
φ
3
4
π
r>
2
p
(
r
)d
r
+
O
(
φ
2
)
.
(2.40)
In Brady & Vicic (1995), the asymptotic expansion of
p
to
O
(
Pe
2
) has been
determined with and without hydrodynamic interactions between particles. Without
hydrodynamic interactions, the contact value of
p
is given by
p
(2) =
−
Pe
2
3
n
·
ˆ
E
·
n
+
Pe
2
8
27
n
·
ˆ
E
·
ˆ
Ω
·
n
+
8
63
n
·
ˆ
E
·
ˆ
E
·
n
+
128
945
ˆ
E
:
ˆ
E
+
O
(
Pe
5
/
2
)
,
where
ˆ
Ω
ij
=(
̇
Γ
ij
−
̇
Γ
ji
)
/
2
̇
γ
is the dimensionless bulk vorticity tensor of the imposed
linear flow and the notation
n
i
ˆ
E
ik
ˆ
Ω
kj
n
j
for the first term in brackets is employed.
Since the static structure factor can be written as
S
0
=1
−
8
φ
+
S
0
+
O
(
φ
2
), and in
a weakly sheared suspension
S
0
=
o
(
φ
) (see the Appendix), the first term in (2.39)
can be re-written as
1
S
0
(1 +
S
0
)
I
=1+8
φ
+
O
(
φ
2
)
.
Therefore
D
c
becomes
D
c
=(1+8
φ
)
I
−
φ
Pe
32
15
ˆ
E
+
φ
Pe
2
128
135
ˆ
E
·
ˆ
Ω
+
128
315
ˆ
E
·
ˆ
E
+
1024
945
(
ˆ
E
:
ˆ
E
)
I
+
O
(
φ
Pe
5
/
2
,φ
2
)
.
(2.41)
While a uniformly valid asymptotic expansion of
p
(
r
)beyond
O
(
Pe
2
) may include
terms with non-integer powers of
Pe
due to non-negligible convection effects, it was
found that the ‘outer’ solution is slaved to the ‘inner’ expansion, and thus the contact
value
p
(2) possesses an expansion in integer powers of
Pe
up to
O
(
Pe
5
/
2
).
318
A. M. Leshansky, J. F. Morris and J. F. Brady
2.3.
Hydrodynamic interactions
We now turn to the more difficult problem of evaluating the collective diffusivity with
hydrodynamic interactions. We shall only consider the analytically tractable case of
a weakly sheared suspension,
i.e.
the low-P
́
eclet-number limit. In the work we make
use of certain prior results on the structural distortion induced by weak shear flow.
The added complexity of hydrodynamics in the dilute theory is largely the result
of long-range hydrodynamic interactions, with expected effects on the integrals of
velocity fluctuations. Nonetheless, results will be found to be qualitatively similar to
the results in the absence of hydrodynamics presented in the previous section, with
one important exception being a weakly divergent (as log
k
for
k
→
0) diffusivity with
tensor form given by the rate of strain, and hence having an impact only on the
off-diagonal diffusivity terms in simple shear flow.
From (2.27) after neglecting the direct third-particle effects we have
D
c
=
1
S
0
I
+
φ
3
4
π
(
D
11
+
D
21
)d
r
+
φ
3
4
π
1
S
0
(
D
11
+
D
21
)
p
d
r
−
r
(
D
22
−
D
21
)
·
∇
p
d
r
−
F
kk
+
O
(
φ
2
)
,
(2.42)
where
p
(
r
) denotes the entire
Pe
-dependent uniformly valid perturbation to the
equilibrium microstructure, and
F
kk
represents the
O
(
k
2
) tensorial coefficient of the
small-
k
expansion of the velocity fluctuation integral:
φ
Pe
3
4
π
1
S
0
i
k
·
r>
2
U
2
(e
i
k
·
r
−
1) d
r
.
(2.43)
The integral in the first square brackets in (2.42) is not convergent as it represents the
sedimentation velocity and must be properly renormalized, as shown by Batchelor
(1976). Together with the excluded volume
O
(
φ
)partof
S
0
, the first square brackets
give the equilibrium collective diffusivity,
D
c
eq
(
Pe
=0)=
1
1
−
8
φ
I
+
φ
3
4
π
(
D
11
+
D
21
)d
r
=
I
(1 + 1
.
47
φ
)+
O
(
φ
2
)
.
(2.44)
For a weakly sheared dilute suspension,
S
0
=1
−
8
φ
+
S
0
+
O
(
φ
2
), and the first
square brackets in (2.42) can be rewritten as
1
S
0
I
+
φ
3
4
π
(
D
11
+
D
21
)d
r
=
I
(1 + 1
.
47
φ
)
−
I
S
0
+
O
(
φ
2
)
.
(2.45)
The third integral on the right-hand side of (2.42) can be integrated by parts to give
r
(
D
22
−
D
21
)
·
∇
p
d
r
=
−
8
r
=2
n
·
(
D
22
−
D
21
)
n
p
(2)d
Ω
+
S
∞
n
·
(
D
22
−
D
21
)
n
r
3
p
(
r
)d
Ω
−
p
r
∇
·
(
D
22
−
D
21
)d
r
−
(
D
22
−
D
21
)
p
d
r
−
I
p
d
r
,
(2.46)
where we have expanded the scaled isolated-particle diffusivity as
D
22
=
I
+
D
22
.The
first integral is zero since
D
22
−
D
21
=
D
r
/
2 and the relative diffusivity,
D
r
, vanishes
Collective diffusion in sheared colloidal suspensions
319
at contact. The surface integral at infinity is also zero, as
p
decays exponentially at
large
r
Pe
−
1
/
2
almost everywhere.
Finally, the integral containing
D
22
in (2.46) can be combined with that of
D
11
in
(2.42) to give
(
D
11
+
D
21
)
p
d
r
−
r
(
D
22
−
D
21
)
·
∇
p
d
r
=2
D
11
p
d
r
+
1
2
p
r
∇
·
D
r
d
r
+
I
p
d
r
.
(2.47)
All integrals in (2.47) are absolutely convergent. Moreover, the regular asymptotic
expansion of
p
up to
O
(
Pe
2
) (Brady & Vicic 1995) can be used in the first two integrals
since
D
11
∼
r
−
4
and
∇
·
D
r
∼
r
−
5
, while
p
1
∼
r
−
3
and
p
2
∼
r
−
1
. The last term in (2.47),
if multiplied by (3
/
4
π
)
φ/S
0
, and using (2.40) gives
I
S
0
+
o
(
φ
). This term requires
the uniformly valid expansion of
p
(see Appendix A), but fortunately, it cancels out
exactly with the last term in (2.45) and
D
c
0
up to
O
(
Pe
2
) can be calculated through
the regular asymptotic expansion of
p
. Hence, the resulting gradient diffusivity (2.42)
will have the form
D
c
=
I
(1 + 1
.
47
φ
)+
φ
3
4
π
2
D
11
p
d
r
+
1
2
p
r
∇
·
D
r
d
r
−
F
kk
+
O
(
φ
2
)
.
(2.48)
Now we turn to the velocity fluctuation integral (2.43) which is divergent if
hydrodynamic interactions are considered in a pairwise additive manner. Expanding
the phase factor e
i
k
·
r
for small
k
leads to
F
kk
in terms of the integral of
U
2
r
over the whole space, and for neutrally buoyant spheres in linear flows
U
2
=
(
U
2
−
U
1
)
−
̇
Γ
·
r
∼
r
−
2
. Therefore, a small-
k
expansion under the integral sign
in (2.43) leads to a divergent integral. The divergence of ensemble averages due
to long-range hydrodynamic interactions in unbounded suspensions is well-studied
(Batchelor 1972; Hinch 1977). Thus, the velocity fluctuation integral needs to be
re-normalized in the standard fashion for non-convergent hydrodynamic interactions.
Note that this term would be non-convergent even if we had not expanded e
i
k
·
r
(meaning it is not the weighting with
r
that causes difficulty), and we cannot simply
add velocity fluctuations in a pairwise manner. Following Brady
et al.
(1988), we write
the velocity of particle 2 due to the motion of particle 1 averaged over all positions
of particle 2 to obtain
i
k
·
r>
2
(
U
2
−
U
,
∞
2
)(e
i
k
·
r
−
1) d
r
+
kk
:
r<
2
rU
,
∞
2
d
r
,
where we have expanded the exponential in the last integral. This integral arises from
the ‘backflow’ being defined over the whole space, not just that outside the particle
at
r
= 2. Here,
U
,
∞
2
is the far-field part of the velocity disturbance that decays as
r
−
2
.
The difference
U
2
−
U
,
∞
2
decays as
r
−
4
and the integral is absolutely convergent. In
this case we cannot simply expand the exponential in the first term as the
r
weighting
results in a conditionally convergent integral. We can, however, extract terms that
decay faster than
r
−
4
by defining
U
2
−
U
,
∞
2
=
U
,
0
2
−
U
,
∞
2
+
U
2
,
(2.49)
320
A. M. Leshansky, J. F. Morris and J. F. Brady
where
U
,
0
2
contains terms up to
r
−
4
and
U
2
is the remainder. Thus, the velocity
fluctuation integral becomes finally
i
k
·
r>
2
U
,
0
2
−
U
,
∞
2
(e
i
k
·
r
−
1)d
r
+
kk
:
r<
2
rU
,
∞
2
d
r
−
kk
:
r>
2
rU
2
d
r
,
(2.50)
where all integrals are absolutely convergent. Other than the integral involving the
r
−
4
dependence of
U
2
, the integrations to evaluate
D
c
are straightforward. This integral
may be written
i
k
·
r>
2
U
,
0
2
−
U
,
∞
2
(e
i
k
·
r
−
1) d
r
=i
k
·
r>
2
r
·
E
·
A
(
−
5)
(
r
)
−
B
(
−
5)
(
r
)
nn
+
B
(
−
5)
(
r
)
I
e
i
k
·
r
d
r
,
(2.51)
noting that, due to symmetry, the discarded term (involving the
−
1) vanishes. Here
A
(
−
5)
and
B
(
−
5)
are the
r
−
5
far-field dependences of the hydrodynamic functions
A
(
r
)and
B
(
r
) specifying the radial dependence of the velocity disturbance in a pair
interaction (Batchelor & Green 1972)
(
U
2
−
U
1
)
−
̇
Γ
·
r
=
−
r
·
E
·
[
A
(
r
)
nn
+
B
(
r
)(
I
−
nn
)
]
,
and for identical spheres of unit radius, these are
A
(
−
5)
=
−
8
r
−
5
,
and
B
(
−
5)
=
16
3
r
−
5
.
Substituting these values into (2.51) we arrive at
i
k
j
ˆ
E
lm
40
3
r>
2
n
j
n
l
n
m
r
4
e
i
k
·
r
d
r
−
16
3
r>
2
δ
mj
n
l
r
4
e
i
k
·
r
d
r
.
(2.52)
Recognizing that the integrals over the angular dependence are proportional to
kkk
and
k
I
, the entire result must be proportional to
k
·
ˆ
E
·
k
. To perform the integrations
we write the integrals as
i
k
j
ˆ
E
lm
r>
2
δ
mj
n
l
r
4
e
i
k
·
r
d
r
=
k
l
ˆ
E
lj
k
j
K
(
k
)
k
,
(2.53)
i
k
j
ˆ
E
lm
r>
2
n
j
n
l
n
m
r
4
e
i
k
·
r
d
r
=
−
k
j
ˆ
E
lm
∂
3
∂k
j
∂k
l
∂k
m
M
(
k
)
=
−
k
l
ˆ
E
lm
k
m
M
k
3
−
M
k
2
,
(2.54)
where
K
(
k
)=
r>
2
e
i
k
·
r
r
5
d
r
and
M
(
k
)=
r>
2
e
i
k
·
r
r
7
d
r
.
The evaluation of expressions involving both these functions in (2.53)–(2.54) is
straightforward and, for instance, for
K
yields
K
(
k
)
k
=
1
k
r>
2
i
μ
e
i
krμ
r
4
d
r
=
4
π
k
3
∞
r
=2
kr
cos
kr
−
sin
kr
r
4
d
r
=
4
π
3
α
3
[
α
cos
α
+
α
3
Ci(
α
)
−
(1 +
α
2
)sin
α
] (2.55)
Collective diffusion in sheared colloidal suspensions
321
where
μ
is the cosine of the angle between
k
and
n
,
α
=2
k
and Ci denotes the cosine
integral function. Taking the limit
k
→
0wefind
K
(
k
)
k
=
4
π
3
ln
k
+ln2+
γ
E
−
4
3
+
O
(
k
2
)
,
where
γ
E
0
.
5772 is the Euler constant. Analogously,
M
k
3
−
M
k
2
=
4
π
15
ln
k
+ln2+
γ
E
−
23
15
+
O
(
k
2
)
.
Inserting these expressions into (2.53)–(2.54) and gathering according to (2.52) we
finally obtain that for small
k
i
k
·
φ
r>
2
(
U
,
0
2
−
U
,
∞
2
)(e
i
k
·
r
−
1) d
r
=
−
32
π
3
φ
ln
k
+ln2+
γ
E
−
7
5
k
·
ˆ
E
·
k
+
O
(
k
4
)
.
(2.56)
Noting that
U
,
∞
2
=
−
5
r
−
2
n
·
ˆ
E
·
nn
=
−
5
r
−
2
γ
r
n
,
the second integral in (2.50) over the excluded volume
r<
2 yields
r<
2
rU
,
∞
2
d
r
=
−
16
π
3
ˆ
E
.
(2.57)
The last integral in (2.50), evaluated numerically, yields
r>
2
rU
2
d
r
=
−
ˆ
E
:
nnnn
d
Ω
∞
2
(
A
−
B
)
r
4
d
r
−
ˆ
E
·
nn
d
Ω
∞
2
B
r
4
d
r
=
−
8
π
15
ˆ
E
(12
.
5)
−
4
π
3
ˆ
E
(0
.
23) =
−
4
π
3
ˆ
E
(5
.
23)
,
(2.58)
where
A
(
r
)and
B
(
r
) are the
o
(
r
−
5
) truncated hydrodynamic functions according to
(2.49). Thus, collecting all the contributions in (2.50) we finally obtain the contribution
of the velocity fluctuation integral to the collective diffusivity:
φ
Pe
3
4
π
1
S
0
i
k
·
U
2
(e
i
k
·
r
−
1) d
r
=
−
8
φ
Pe
(ln
k
+ln2+
γ
E
−
1
.
55)
k
·
ˆ
E
·
k
+
O
(
k
2
,φ
2
)
.
(2.59)
Note that the abnormal scaling
∼
ln
k
is affecting the components of
D
c
which have
the geometry of the imposed flow
ˆ
E
. Therefore, in the case of a simple shear flow it
only influences the off-diagonal collective-diffusion coefficient
D
c
12
(
D
c
13
=
D
c
23
=0by
symmetry), while the diagonal components of
D
c
ii
are truly diffusive, as we shall see
below.
To complete the asymptotic evaluation of
D
c
at
O
(
Pe
) we also must evaluate
integrals on the right-hand side of (2.47):
r>
2
p
1
D
11
d
r
=
−
ˆ
E
:
nnnn
d
Ω
∞
2
q
(
r
)
x
a
11
(
r
)
−
y
a
11
(
r
)
r
2
d
r
=
8
π
15
ˆ
E
(0
.
44) (2.60)
and
r>
2
p
1
r
∇
·
D
r
d
r
=
−
ˆ
E
:
nnnn
d
Ω
∞
2
q
(
r
)
Z
(
r
)
r
3
d
r
=
8
π
15
ˆ
E
(
−
2
.
57)
,
(2.61)
322
A. M. Leshansky, J. F. Morris and J. F. Brady
where
q
(
r
) is given in (A 37). In (2.61) we exploited the fact that
∇
·
D
r
≡
n
Z
(
r
)=
n
G
(
r
)+2
G
(
r
)
−
H
(
r
)
r
.
Here
x
a
11
,
y
a
11
,
Z
(
r
),
G
(
r
)and
H
(
r
) are known mobility functions (Jeffrey & Onishi
1984; Kim & Karrila 1991).
Gathering all contributions, the collective diffusivity is given to
O
(
Pe
)by
D
c
=(1+1
.
47
φ
)
I
+ 8 (ln
k
+
C
)
φ
Pe
ˆ
E
+
O
(
φ
2
,
Pe
2
)
,
(2.62)
where
C
=ln2+
γ
E
−
1
.
64
−
0
.
37. Note that the coefficient of the term linear in
ˆ
E
is negative for small
k
, in agreement with the sign found without hydrodynamic
interaction.
To complete the evaluation of
D
c
to
O
(
Pe
2
), it follows from (2.48) that we also
must determine
φ
Pe
2
3
4
π
2
D
11
p
2
(
r
)d
r
+
1
2
r
∇
·
D
r
p
2
(
r
)d
r
+
rU
2
p
1
(
r
)d
r
,
(2.63)
where the last integral comes from small-
k
expansion of the velocity fluctuation
integral (2.43). The first two integrals are convergent as discussed just after (2.47),
while the last integral is convergent since
p
1
decays as
r
−
3
for large
r
and
U
2
∼
r
−
2
.It
is expected from the form of
p
1
,
p
2
and
U
2
that the
O
(
Pe
2
) contribution is quadratic
in the imposed flow and is a linear combination of
ˆ
E
·
ˆ
Ω
,
ˆ
E
·
ˆ
E
and (
ˆ
E
:
ˆ
E
)
I
as was the
case without hydrodynamics (2.41). The constant coefficients are to be determined by
numerical integration. The velocity fluctuation integral in component form is
r
i
U
2
j
p
1
(
r
)d
r
=
ˆ
E
kl
ˆ
E
rs
n
i
n
j
n
k
n
l
n
r
n
s
d
Ω
(
A
(
r
)
−
B
(
r
))
q
(
r
)
r
4
d
r
+
ˆ
E
jk
ˆ
E
rs
n
i
n
k
n
r
n
s
d
Ω
B
(
r
)
q
(
r
)
r
4
d
r
=4
π
2
105
ˆ
E
ks
ˆ
E
ks
δ
ij
I
1
+
8
105
I
1
+
2
15
I
2
ˆ
E
ik
ˆ
E
kj
,
(2.64)
where the numerical evaluation of the radial integrals in the two terms on the right-
hand side yields
I
1
=
∞
2
(
A
(
r
)
−
B
(
r
))
q
(
r
)
r
4
d
r
=20
.
90 and
I
2
=
∞
2
B
(
r
)
q
(
r
)
r
4
d
r
=
2
.
36, respectively, and the surface integrals are evaluated using the identity
1
4
π
n
α
1
n
α
2
...n
α
m
m even
d
Ω
=
1
(
m
+ 1)!!
p
δ
α
1
α
2
δ
α
3
α
4
...δ
α
m
−
1
α
m
,
where
p
is the number of all possible index permutations. The first two integrals in
(2.63) can be written in the form
nn
p
2
(
r
)
T
(1
,
2)
(
r
)d
r
,
(2.65)
where
T
(1)
=
x
a
11
(
r
)
−
y
a
11
(
r
)and
T
(2)
=
rZ
(
r
) for the first and second integral,
respectively. Substituting the known expansion for
p
2
(
r
)intermsof
h
2
(
r
)–
h
5
(
r
)
(Brady & Vicic 1995, Appendix) and using the above identity for integrals over the