J. Fluid Mech.
(2004),
vol.
506,
pp.
285–314.
c
2004 Cambridge University Press
DOI: 10.1017/S0022112004008651 Printed in the United Kingdom
285
Shear-induced self-diffusion in non-colloidal
suspensions
By ASIMINA SIEROU
AND
JOHN F. BRADY
Division of Chemistry and Chemical Engineering, California Institute of Technology,
Pasadena, CA 91125, USA
(Received
9 August 2001 and in revised form 17 November 2003)
Self-diffusion in a monodisperse suspension of non-Brownian particles in simple
shear flow is studied using accelerated Stokesian dynamics (ASD) simulation. The
availability of a much faster computational algorithm allows the study of large systems
(typically of 1000 particles) and the extraction of accurate results for the complete
shear-induced self-diffusivity tensor. The finite, and often large, autocorrelation time
requires the mean-square displacements to be followed for very long times, which is
now possible with ASD. The self-diffusivities compare favourably with the available
experimental measurements when allowance is made for the finite strains sampled
in the experiments. The relationship between the mean-square displacements and
the diffusivities appearing in a Fokker–Planck equation when advection couples to
diffusion is discussed.
1. Introduction
Shear-induced particle diffusion plays an important role in the behaviour of
concentrated non-Brownian suspensions and is responsible for a variety of interesting
rheological phenomena in such systems. The first direct experimental study of shear-
induced diffusion was reported by Eckstein, Bailey & Shapiro (1977); a variation
of their experimental technique was later used by Leighton & Acrivos (1987) and
Phan & Leighton (1999), followed by the introduction of an alternative experimental
method by Breedveld
et al.
(1998, 2001). Most of the aforementioned experiments
were limited to the determination of the self-diffusion coefficients in the velocity
gradient and vorticity directions only.
Theoretical work on shear-induced diffusion is complicated by the fact that with
only hydrodynamic interactions, the simplest two-body system either does not exhibit
diffusive behaviour (in the velocity-gradient and vorticity directions) or exhibits
a singular behaviour (an infinite diffusivity in the direction of the fluid velocity).
Acrivos
et al.
(1992) studied the diffusion coefficient parallel to the fluid velocity
by introducing a mechanism (an additional pair of particles) that results in a finite
diffusivity, while Wang, Mauri & Acrivos (1996) studied the diffusion perpendicular to
the fluid velocity by examining three-particle interactions. Da Cunha & Hinch (1996)
studied the two-particle problem in the presence of surface roughness – necessary to
create diffusive motion – and Brady & Morris (1997) introduced residual Brownian
motion and a hard-sphere interparticle force to break the symmetry of the two-particle
limit. Marchioro & Acrivos (2001) and Drazer
et al.
(2002) used Stokesian dynamics
simulations to study shear-induced diffusion in concentrated suspensions and ad-
dressed the question of whether self-diffusion can be present in a purely deterministic
system (e.g. in the absence of any roughness, residual Brownian motion, etc.).
286
A. Sierou and J. F. Brady
The purpose of this work is to use a recently developed algorithm for the calculation
of hydrodynamic interactions among particles (Sierou & Brady 2001) that allows the
simulation of much larger systems than were previously possible. Although calculation
of other rheological properties (e.g. viscosity) has proved to be accurate even with
quite small systems, the calculation of the self-diffusion coefficient requires significantly
larger systems. Access to systems with numbers of particles of the order of 1000 allows
accurate calculation of the self-diffusion coefficient and a study of its dependence on
the number of particles and the volume fraction.
The complete diffusion tensor can now be determined and simulation results for the
off-diagonal components and the coefficient parallel to the fluid velocity are presented
for the first time for a non-Brownian system. The presence of an external velocity
field makes the determination of diffusive displacements in the flow, or longitudinal,
direction difficult, and only very limited experimental results are available for these
components (Breedveld 2000). Foss & Brady (1999) calculated the longitudinal self-
diffusion coefficient for a system of Brownian particles at very high P
́
eclet numbers
by following the theoretical analysis of Morris & Brady (1996) which suggested that
the affine shearing motion could be subtracted directly at each instant in time and
the diffusivity could be calculated by taking into account only the non-affine displace-
ments. Unfortunately, as we shall show, this is not correct and the affine displacement
couples with the non-affine displacement to give an additional contribution to the
longitudinal long-time self-diffusivity. The calculation of the diffusion coefficient in
the flow direction is still a formidable task for dilute suspensions, however, because
of the strong influence of the periodicity on the particle displacements.
When advection couples to diffusion, the mean-square displacements at long time
may no longer grow linearly in time. For example, in simple shear flow, the mean-
square displacement in the flow direction grows as
t
3
, with a coefficient proportional
to the diffusivity in the velocity-gradient direction (Elrick 1962). Because of this strong
temporal growth, it can be very difficult to extract the underlying diffusive motion
in the flow direction from measurements of the mean-square displacement. We show
from an analysis of the Fokker–Planck equation both how to define the underlying
shear-induced diffusivities and how to extract their values from either simulation or
experiment.
In the next section, we outline the simulation method and the parameters used.
In
§
3.1, we discuss how the diffusion coefficients are determined from both the
mean-square displacements and the velocity autocorrelation function. The results for
the diffusivities in the velocity-gradient and vorticity direction are presented and
compared with experiment. Discrepancies between simulation and some experimental
results can be attributed to the fact that most of the experimental techniques attempt
to determine the self-diffusivities for very small strains where the behaviour is not yet
diffusive; such a case is presented in detail in
§
3.3. In
§
4, the diffusivity in the fluid
velocity direction is discussed along with the only non-zero off-diagonal component of
the diffusivity tensor. Finally, we conclude in
§
5 with a brief summary and suggestions
for future work.
2. The simulation method
A new implementation of the conventional Stokesian dynamics (SD) method, called
accelerated Stokesian dynamics (ASD), is used in these simulations. The methodology
of ASD has been laid out in Sierou & Brady (2001) where there is a full discussion.
(See also Sierou & Brady 2002 for an application of ASD to rheology.)
Shear-induced self-diffusion in non-colloidal suspensions
287
We consider the Stokes flow of a suspension of hard spheres in a fluid of viscosity
η
undergoing simple shear with shear rate
̇
γ
. As has been demonstrated in the past
(Melrose & Ball 1995; Dratler & Schowalter 1996), in the absence of Brownian
motion, a repulsive interparticle force is necessary to prevent the formation of infinite
clusters that ‘jam’ and result in excessive particle overlaps. Although the choice of
this force is somewhat arbitrary, its introduction captures the behaviour of a real
physical system, since the presence of residual Brownian forces or particle roughness,
would also prevent the particles from touching. This force was chosen to be of the
form (Bossis & Brady 1984):
F
p
(
αβ
)
=
F
0
τ
e
−
τ
1
−
e
−
τ
e
(
αβ
)
,
(2.1)
where 6
π
ηa
2
̇
γ
F
p
(
αβ
)
is the force exerted on sphere
α
by sphere
β
,
F
0
is a dimensionless
constant representing the magnitude of the interparticle force,
τ
is related to its range,
=
r
−
2 is the spacing between the surfaces of the two spheres and
e
(
αβ
)
is the unit
vector connecting the centres of the two spheres. In the results reported here, the
value of
τ
was set at 1000, while
τF
0
was chosen to be unity. Although the exact
magnitude and form of this interparticle force can play an important role in the
values of some rheological properties (e.g. viscosity and normal stress differences),
especially for high-volume fractions, simulations using other parameters gave very
similar results for the self-diffusion tensor. (At very low-volume fractions, however,
where two-particle interactions dominate, the range and amplitude of the repulsive
force can strongly influence the self-diffusivities (Drazer
et al.
2002)).
All lengths are non-dimensionalized by the particle radius
a
, all times by the charac-
teristic time
̇
γ
−
1
.Thetimestepusedrangedfrom
t
=5
×
10
−
3
to 5
×
10
−
4
depending
on the volume fraction (higher-volume fractions required smaller time steps). A fourth-
order Adams–Bashforth integration scheme was used to update the particle positions.
3. The self-diffusivities in the velocity gradient (
D
yy
) and vorticity
(
D
zz
) directions
3.1.
Numerical determination of the self-diffusivities
We now describe the calculation of the self-diffusion coefficients in the velocity-
gradient and vorticity directions. (The calculation of the diffusivity in the fluid velocity
direction is subject to further limitations and will be described in the following section.)
Dimensional arguments show that the self-diffusivities scale as
̇
γa
2
, as these are the
only time and length scales in the system (in the absence of interparticle forces). All
diffusivities reported are therefore normalized by
̇
γa
2
.
Let
x
,
y
and
z
denote the directions of the flow, the velocity-gradient and the
vorticity, respectively. Then, for a suspension undergoing simple shear flow in the
x
-direction with velocity gradient in the
y
-direction, the mean-square displacements
will grow with time according to
y
(
t
)
y
(
t
)
∼
2
D
yy
t,
z
(
t
)
z
(
t
)
∼
2
D
zz
t,
at long times after diffusive motion has been established. Here, and hereinafter,
the angle brackets denote an average over all particles in the system, while
D
yy
and
D
zz
denote the self-diffusivities in the velocity-gradient and vorticity directions,
respectively. The self-diffusivities are defined as the time rate of change of half the
mean-square displacements, and can be calculated in a straightforward manner by a
number of different, but equivalent, methods.
288
A. Sierou and J. F. Brady
16
12
8
4
800
600
400
200
0
yy
(
D
yy
= 8.50
×
10
–3
)
zz
(
D
zz
= 4.10
×
10
–3
)
yy
(
a
)
γ
t
.
a
2
zz
a
2
10
1
10
0
10
–1
10
–2
10
–3
10
–4
10
–1
10
0
10
1
10
2
t
t
2
(
b
)
yy
a
2
zz
a
2
Figure 1.
The mean-square displacement curves
yy
and
zz
as a function of time (or strain)
̇
γt
for a system of
N
= 1000,
φ
=0
.
20 for strains up to 800, in (
a
) a linear and (
b
) a log–log
plot. The values of the resulting diffusivities are also given. In this and subsequent figures, all
mean-square displacements are in units of
a
2
, and all diffusivities are in units of
̇
γa
2
.
First, the mean-square displacement curve can be constructed for a long simulation
run. Using a large number of particles, this approach results in a mean-square
displacement curve that is linear in time at long times, as shown in figure 1 for a
system with
N
= 1000 at volume fraction
φ
=0
.
2. It is clear from the logarithmic plot
that, at very short times, the mean-square displacement shows a quadratic temporal
behaviour, while after a strain of about 20 the final linear regime sets in and continues
until the end of the run. The quadratic behaviour corresponds to the regime where the
particle motion is deterministic and still strongly correlated – the diffusive behaviour
has not yet been established. It should be noted that, although after a strain of about
Shear-induced self-diffusion in non-colloidal suspensions
289
0.8
0.6
0.4
0.2
30
25
20
15
10
5
0
yy
(
D
yy
= 8.26
×
10
–3
)
zz
(
D
zz
= 3.97
×
10
–3
)
yy
(
a
)
a
2
zz
a
2
10
–6
10
–2
10
–1
10
0
10
1
10
2
10
–4
10
–2
10
0
γ
t
t
t
2
.
(
b
)
yy
a
2
zz
a
2
Figure 2.
The mean-square displacement curves
yy
and
zz
as a function of time (or strain)
̇
γt
for a system of
N
= 1000,
φ
=0
.
20 averaged over a total strain of (
a
)30and(
b
) 100. The
values of the resulting diffusivities are also given.
one the quadratic scaling no longer holds, a transitional region exists up to strains of
at least 10 before a clear linear behaviour is evident. The slope of the resulting curve
for
t>
20 can be readily evaluated and the resulting diffusivity is given in figure 1.
Second, this same long simulation run can be split into shorter independent runs,
and the resulting mean-square displacement curves from each shorter run averaged.
Although an infinite time limit is implied for obtaining the long-time self-diffusivity,
the particle does not necessarily need to travel a very long distance to reach this long-
time asymptote; it simply needs to encounter enough particle–particle interactions to
sample the dynamic microstructure and ‘forget’ its initial configuration and velocity.
Figure 2 is for the same simulations as in figure 1, but now over a shorter time interval.
290
A. Sierou and J. F. Brady
The long run of the previous case was split into smaller runs of 30–100 strains, the
results averaged and presented again in linear and logarithmic plots. Because of the
increased statistical data, the slope of the mean-square displacement curve can be
evaluated with greater accuracy than before. The final result for the diffusivity is
very close to the value reported in figure 1 and within the statistical accuracy of
the data. Although the linear behaviour of the curves in the short and long runs is
identical, the initial transient behaviour, as is apparent from the logarithmic plots, can
be quantitatively different (qualitatively, a
t
2
behaviour is apparent for both cases).
This is because for the shorter runs the sequence of independent configurations are
all sampled after the system has reached its steady-state configuration, while the
beginning of the overall mean-square displacement curve in figure 1 corresponds to
an initially random configuration. As a result, a different quadratic behaviour is seen.
Finally, the diffusivities can be evaluated from the integral of the velocity
autocorrelation function using the well-known expression:
D
yy
=
∫
∞
0
lim
t
→∞
u
y
(
τ
+
t
)
u
y
(
t
)
d
τ,
where
u
y
is replaced by
u
z
for
D
zz
. The time it takes for the velocity autocorrelation
to decay to zero corresponds to the time it takes the mean-square displacement curve
to reach the linear regime; for the diffusivities to be evaluated correctly, the integral
must be computed for at least that minimum time. As shown by Marchioro & Acrivos
(2001), the velocity autocorrelation has a negative contribution that must be included
to obtain the correct diffusivity. For the case shown in figure 1, the integral of the
velocity autocorrelation function gave
D
yy
=8
.
30
×
10
−
3
and
D
zz
=3
.
96
×
10
−
3
.
We could also improve the statistics by averaging over a number of realizations
or initial conditions; however, with the large number of particles used and the long
simulation times, our results appear to be sufficiently accurate.
The three approaches to determine the self-diffusivity gave consistent results for
all volume fractions studied. The time necessary for the diffusive behaviour to be
reached is a strong function of the volume fraction, and, as has been reported in
the past, increases with decreasing volume fraction. For very dilute suspensions, the
particle interactions are limited and a larger strain is necessary until the particles have
sampled a sufficiently large number of independent collisions and the linear regime
is reached.
Before we proceed to the presentation of the majority of our simulation results,
it is worth mentioning briefly how we define their statistical properties. The values
of the diffusivities obtained from the three aforementioned approaches are in all
cases in excellent agreement, so their simple average is used as the reported value
of the diffusivity. The calculation of the standard deviation, however, needs to be
discussed in more detail. The error associated with the calculation of any element of
the self-diffusion tensor can result from one (or more) of the following sources. First,
errors derived from the solution of the equations of motion – the assumptions used
in the development of ASD combined with the ‘standard’ numerical errors present
in any numerical calculation (finite time steps, machine accuracy, etc.) – always lead
to a numerical error in the calculation of any property. The statistical nature of
the particle displacements and the fact that the magnitude of the numerical error is
usually small makes the numerical inaccuracies of very little importance and therefore
we shall not consider them further. Secondly, there is the uncertainty associated with
the number of independent mean-square displacement curves available. The mean-
square displacement curves shown in figures 1 and 2 are already averaged over
Shear-induced self-diffusion in non-colloidal suspensions
291
all particles, and, for the case of figure 2, all the shorter independent runs are
also averaged. In order to obtain the correct values for the self-diffusion tensor,
it must be assured that a sufficient number of mean-square displacement curves
have been sampled, or in other words, that all possible particle interactions have
been sampled. This can be assured by examining systems with large
N
, or a large
number of configurations of systems with smaller
N
, as has been suggested in
the past (Foss & Brady 1999). The calculation of the error then corresponds to
grouping particles/configurations together and calculating the standard deviation of
the resulting measurements. This is the error routinely reported in our results and
is calculated either as a temporal average over smaller segments of the mean-square
displacement curve for the long runs of figure 1, or as an average over different groups
of configurations for the case of figure 2. (Note that figure 2 presents the average
mean-square displacement curve over all possible configurations and particles for a
given run; only the mean can be calculated from this curve.) Finally, the calculation
of the slope of the mean-square displacement curve, or the numerical integration of
the velocity autocorrelation function, introduces a numerical error. In addition to
any uncertainty associated with the numerical calculation, the most important issue
is the correct time regime (or strain) over which to determine the diffusivities. Very
long tails in the velocity autocorrelation function can be erroneously discarded and
visual observation of when the process becomes linear can be misleading and result
in large errors. Although we believe that this source of error can be potentially very
important, we have tried to make a conservative estimate in our calculation of the
minimum strain needed for diffusive behaviour to set in, and are confident that the
error from this source is not significant. Therefore, for the remainder of this paper
the error bars reported correspond to uncertainties associated with simply the volume
of the statistical data, uncertainties that for most cases are less that 10%. Finally, the
last source of uncertainty, system size and the effect of the periodic unit cell, will be
discussed in more detail below.
3.2.
The dependence on the number of particles
N
The introduction of a new faster methodology allows the calculation of the self-
diffusion coefficients for numbers of particles of the order of 10
3
, and we can therefore
present a reliable analysis of the dependence of the diffusivities on the number of
particles,
N
. Figure 3 shows the dependence of
D
yy
on
N
, for the case of
φ
=0
.
35 (the
interparticle force, defined in (2.1) with the parameters mentioned, is always present
unless otherwise noted). There is a rather sharp increase in the diffusivity for small
values of
N
(
N<
100), followed by a levelling off. The same trend is apparent for
all elements of the self-diffusion tensor. No theoretical expression is available for the
N
dependence of the long-time self-diffusivity, and so in figure 3 we simply present
a
N
−
1
best fit that seems to describe the data satisfactorily. Plotting the individual
mean-square displacements (not shown) for different numbers of particles also verifies
that for
N>
100 all curves are statistically indistinguishable. A possible explanation
of this behaviour may lie in the fact that after the size of the cell is sufficiently large
(or
N
is sufficiently large), all of the important particle interactions are included and
therefore the number of particles plays no role (other than increasing the statistical
accuracy, i.e. giving curves that appear smoother). When the number of particles is
very small, the corresponding size of the cubic unit cell is also small (for
N
= 27,
L
=6
.
86
a
for
φ
=0
.
35). This means that particle interactions for separations greater
than
L
occur through periodic images, which may explain the smaller diffusivities
observed for very small number of particles. (As discussed in
§
4.2 a similar, but much
292
A. Sierou and J. F. Brady
0.06
0.05
0.04
0.03
0.02
0.01
200
400
600
800
1000
0
N
φ
= 0.35
ASD,
N
= 64–1000
SD,
N
= 16–64
best fit for
D
∞
=
D
N
(1+
β
/
N
)
γ
a
2
.
D
yy
Figure 3.
The dependence of
D
yy
on the number of particles
N
, with
N
ranging from 16 to
1000, for
φ
=0
.
35. A very sharp increase is observed for small volume fractions, followed by
a relatively constant value.
φD
yy
D
zz
0
.
10
0
.
0017
±
0
.
0003
0
.
0011
±
0
.
0002
0
.
15
0
.
0045
±
0
.
0006
0
.
0024
±
0
.
0004
0
.
20
0
.
0084
±
0
.
0010
0
.
0040
±
0
.
0006
0
.
25
0
.
0171
±
0
.
0020
0
.
0070
±
0
.
0007
0
.
30
0
.
0310
±
0
.
0040
0
.
0117
±
0
.
0010
0
.
35
0
.
0460
±
0
.
0050
0
.
0185
±
0
.
0020
0
.
40
0
.
0620
±
0
.
0060
0
.
0290
±
0
.
0030
0
.
45
0
.
0583
±
0
.
0070
0
.
0450
±
0
.
0040
0
.
50
0
.
0580
±
0
.
0070
0
.
0520
±
0
.
0050
Table 1.
The dependence of
D
yy
and
D
zz
on the volume fraction
φ
for
N
= 1000.
more apparent, effect of the size of the unit cell is present in the calculation of
D
xx
for low volume fractions.) For the remainder of this paper all results will correspond
to systems of
N
= 1000 particles and no further corrections to the infinite system
limit will be attempted since, as is apparent in figure 3, such corrections would be
very small and within the statistical accuracy of the calculations.
3.3.
The dependence on the volume fraction
Figures 4 and 5 and table 1 show the dependence of
D
yy
and
D
zz
on the volume
fraction
φ
. Although qualitatively most sets of experimental or simulation results are
in agreement, quantitatively discrepancies are present that need to be explained and
examined in more detail. As was mentioned in
§
3.1, the notion that the long-time
self-diffusivities can be calculated for strains that are not particularly long has been
used, and most of the recent experimental and simulation results presented in figures 4
and 5 are calculated over a large number of configurations, but for relatively small
strains (Breedveld
et al.
1998, 2001; Foss & Brady 1999). Unfortunately, as shown
Shear-induced self-diffusion in non-colloidal suspensions
293
0.16
0.12
0.08
0.04
0
0.1
φ
0.2
0.3
0.4
0.5
0.6
Phan & Leighton (1999)
Breedveld
et al.
(1998)
Breedveld
et al.
(2001)
Foss & Brady (1999)
ASD
Eckstein
et al.
(1977)
Leighton & Acrivos (1987)
γ
a
2
.
D
yy
Figure 4.
The dependence of
D
yy
on the volume fraction
φ
. Accelerated Stokesian dynamics
(ASD) results are for a system of
N
= 1000 and compared with a number of experimental
and simulation results.
0.08
0.06
0.04
0.02
0
0.6
0.5
0.4
0.3
0.2
0.1
φ
Phan & Leighton (1999)
Breedveld
et al.
(1998)
Breedveld
et al.
(2001)
Foss & Brady (1999)
ASD
γ
a
2
.
D
zz
Figure 5.
The dependence of
D
zz
on the volume fraction
φ
. Accelerated Stokesian dynamics
(ASD) results are for a system of
N
= 1000 and compared with a number of experimental
and simulation results.
294
A. Sierou and J. F. Brady
0.4
0.3
0.2
0.1
yy
5
4
3
2
1
0
ASD
Breedveld
et al.
(2001)
Breedveld
et al.
(2001)
-- linear best-fit
(
a
)
a
2
yy
a
2
0.6
0.5
0.4
0.3
0.2
0.1
0
20
15
10
5
Strain,
γ
t
ASD
Breedveld
et al.
(2001)
Breedveld
et al.
(2001)
-- best-fit (
D
yy
= 0.021)
ASD best-fit (
D
yy
= 0.0083)
.
(
b
)
Figure 6.
The mean-square displacement
yy
from ASD simulations for
φ
=0
.
20,
N
= 1000
is compared with the equivalent results given by Breedveld
et al.
(2001) for strains of (
a
)5
and (
b
) 20. The resulting diffusivities and slopes of the mean-square displacement curve are
also shown.
below and as was discussed by Marchioro & Acrivos (2001), the strains used for these
studies are too small and therefore do not correspond to the linear regime, as a result
they overestimate the diffusivity.
Figure 6(
a
) is adapted from Breedveld
et al.
(2001) and shows the mean-square
displacement curve,
yy
,for
φ
=0
.
20 for strains up to 5. (We choose to examine
the case of a relatively small volume fraction in detail because, as is apparent in
figure (8) of Breedveld
et al.
(2001), the accuracy in their determination of the mean-
square displacement curve appears to be greatly reduced when the volume fraction
is increased owing to limitations of the experimental technique; this is not observed
Shear-induced self-diffusion in non-colloidal suspensions
295
0.14
0.12
0.10
0.08
0.06
0.04
0.02
0.6
0.5
0.4
0.3
0.2
0.1
0
φ
Breedveld
et al.
(1998)
Breedveld
et al.
(2001)
ASD
ASD -- small strains
γ
.
a
2
D
yy
Figure 7.
The self-diffusion coefficient
D
yy
as calculated from ASD for long and short strains
and compared with the experimental results (corresponding to the same short strains) of
Breedveld
et al.
(2001).
in the numerical simulations.) The dashed curve represents the ‘linear-fit’ performed
by Breedveld
et al.
(2001); the value of the diffusivity that they obtain is also given.
The solid line corresponds to our mean-square displacement for the same volume
fraction for a system of
N
= 1000 particles. The agreement between the two curves is
remarkably good over the entire range of experimentally accessible strains. Although
this is the case for the mean-square displacement curves, the resulting diffusivities are
not in good agreement and the value reported by Breedveld
et al.
(2001) is over twice
as large as the value we calculate. The reason for this is apparent in figure 6(
b
), where
the same mean-square displacement curve is presented but now for strains up to 20.
It is now clear that the linear regime has not yet been reached at a strain of 5, and
strains of at least 10 need to be sampled for the long-time diffusive behaviour to be
reached. In figure 6(
b
), we also show the change in slope in the
yy
curve and how,
if the strain is not sufficiently long, the slope simply corresponds to a transitional
regime between the
t
2
and
t
behaviours and overestimates the diffusivity. The same
conclusion is reached by studying the velocity autocorrelation function: the long tail
present in the velocity autocorrelation function is not taken into account when short
strains are used and, since the velocity autocorrelation function is negative at those
strains, the value of the diffusivity is overestimated (Marchioro & Acrivos 2001).
Figure 7 presents our simulation results, but now instead of calculating the slope of
the displacement curves in the linear regime, we calculate them for the same strains
reported by V. Breedveld (personal communication). We now observe a much better
agreement between the numerical and experimental data for all volume fractions
(although for larger volume fractions, some discrepancies are still present).
Although the limitation on strain in the experiments of Breedveld
et al.
(2001)
explains the difference between simulation and these experiments, this is not
necessarily the case for the other experimental results. Data by Leighton & Acrivos
296
A. Sierou and J. F. Brady
0.12
0.10
0.08
0.06
0.04
0.02
zz
5
4
3
2
1
0
ASD
Breedveld
et al.
(2001)
Strain,
γ
t
.
a
2
Figure 8.
The mean-square displacement
zz
as calculated from ASD simulations for
φ
=0
.
20,
N
= 1000 is compared with the equivalent results given by Breedveld
et al.
(2001) for strains of 5.
(1987) and Phan & Leighton (1999) are calculated for larger strains, but with a
different method that collects data on relatively few particles/realizations, and the
origin of the differences between these experiments and the simulations is not known.
We also wish to briefly comment here on the data by Foss & Brady (1999). In
their approach, very small strains were also used, strains that according to figure 6(
b
)
should overestimate the resulting mean-square displacement coefficient significantly.
As is apparent from figure 4, however, this is not the case and our results are in
general agreement with the results reported by Foss & Brady (1999). We believe that
this is the result of the very small systems used by Foss & Brady (1999);
N
=27for
all of their simulations. As was mentioned in the previous subsection, a very strong
N
dependence is present for such a small number of particles and the diffusivities
are underestimated significantly. It is the combination of these two sources of error –
the very small number of particles and the very small strains used – that cancel each
other and lead to a diffusivity of the correct magnitude.
Figure 8 presents the
zz
mean-square displacement curve for
φ
=0
.
20 as given by
Breedveld
et al.
(2001) compared with our simulation results; again the agreement
is reasonable, although not nearly as good as in figure 7. As shown in figure 9, the
agreement between the short strain simulations and the experimental data is not
as good for higher volume fractions as was the case for
D
yy
.Thisseemstobea
consequence of the different functional dependence on the volume fraction (see below
for more detail); in our data,
D
zz
is an increasing function of the volume fraction
(and as a result so are the values calculated for smaller strains), while Breedveld’s
data show a very sharp increase followed by a plateau. We have no explanation for
this qualitative difference. It seems that although the method developed by Breedveld
et al.
(1998, 2001) has some very attractive characteristics that can reproduce accurate
results for small strains, larger strains are needed to capture the long-time diffusive
behaviour and some adjustments are necessary in the method to allow for larger
observation windows.
Shear-induced self-diffusion in non-colloidal suspensions
297
0.10
0.08
0.06
0.04
0.02
0.6
0.5
0.4
0.3
0.2
0.1
0
φ
Breedveld
et al.
(1998)
Breedveld
et al.
(2001)
ASD
ASD -- small strains
γ
.
a
2
D
zz
Figure 9.
The self-diffusion coefficient
D
zz
as calculated from ASD for long strains, short
strains and compared with the experimental results (corresponding to the same short strains)
given by Breedveld
et al.
(2001).
We now offer a brief discussion of the qualitative behaviour of
D
yy
and
D
zz
as a
function of volume fraction. For small volume fractions,
φ<
0
.
30, a
φ
2
behaviour is
evident for both diffusivities. This is expected since, in the absence of any interparticle
force, the interactions between only two particles lead to zero net displacements (the
particles always return to their original streamlines), and the presence of a third
particle is necessary to generate non-zero displacements leading to a diffusivity of
O
(
φ
2
). As shown by da Cunha & Hinch (1996), however, with an interparticle force, a
net displacement in the velocity-gradient and vorticity directions can be achieved with
only two particles, leading to an
O
(
φ
) dependence of the diffusivity. For the range
of the interparticle force used here and the minimum separations that result, the size
of this
O
(
φ
) correction is extremely small and is overwhelmed by the displacements
caused by the presence of a third particle. Evidently, the three-particle effects are
significant even at the lowest volume fraction studied (
φ
=0
.
10). Much lower-volume
fractions would be necessary to capture and verify the order
φ
correction to the
diffusivity. We briefly note here that it is not clear whether the presence of an
interparticle force is necessary for diffusive motion to appear. It has been suggested
in the past (Marchioro & Acrivos 2001) that even in the absence of an interparticle
force, diffusive-like characteristics could be present owing to the chaotic nature of
the many-body hydrodynamic interactions. If that is indeed the case, the order
φ
correction to the diffusivity might not be dominant even for extremely small volume
fractions (Drazer
et al
. 2002).
As is apparent from figure 4, and in qualitative agreement with most experimental
results, the diffusivity in the velocity-gradient direction appears to plateau at volume
fractions about
φ
=0
.
35
−
0
.
40. This is not the case for the diffusivity in the vorticity
direction however, which is seen to increase roughly as
φ
2
for all volume fractions
studied. The possible exception being
φ
=0
.
50, where the increase is less apparent;
298
A. Sierou and J. F. Brady
the experimental results for this case suggest the presence of an earlier plateau in
the value of
D
zz
. The presence of such a plateau in either
D
yy
or
D
zz
has not been
explained.
4. The longitudinal (
D
xx
) and off-diagonal (
D
xy
) self-diffusivities
4.1.
Theoretical approach
The relationship between the complete diffusivity tensor (including
D
xx
and
D
xy
[=
D
yx
]) and the mean-square displacement for a Brownian particle suspended in a
fluid undergoing simple shear flow has been the subject of numerous studies. Elrick
(1962) calculated the
xx
and
xy
mean-square displacements for such a Brownian
particle and found them to grow in time according to:
xx
=2
D
xx
t
+2
D
xy
̇
γt
2
+
2
3
D
yy
̇
γ
2
t
3
,
(4.1
a
)
xy
=2
D
xy
t
+
D
yy
̇
γt
2
.
(4.1
b
)
Diffusive motion in the
y
-direction,
D
yy
, couples to the advective flow in the
x
-
direction to give a mean-square displacement in the flow direction that grows cubically
in time. In obtaining (4.1), it is assumed that all diffusivities are constant in time and
that the ‘random’ Brownian step the particle takes is diffusive at all times.
Using (4.1) for a system of non-Brownian particles poses further limitations, e.g. the
motion in the
y
-direction is not diffusive at short times and therefore (4.1) can only
be valid at very long times where the
t
3
(or
t
2
) behaviour will dominate and the linear
term
D
xx
t
(or
D
xy
t
) will be difficult to determine accurately. It is therefore desirable
to derive an alternative formulation that not only allows the accurate calculation of
the diffusivities, but also serves as a proof of the validity of (4.1) for a system of
non-Brownian particles. Following the equations of motion for our system, we show
how to calculate the corresponding mean-square displacements and correctly deduce
the diffusion coefficients.
The position of a given particle (initially at the origin) in the
x
-direction is given
by:
x
α
(
t
)=
∫
t
0
U
x,α
(
τ
)d
τ
=
∫
t
0
(
U
h
x,α
(
τ
)+
U
∞
x,α
(
τ
)
)
d
τ
=
X
h
α
+
∫
t
0
̇
γy
α
(
τ
)d
τ,
(4.2)
where,
U
h
x,α
and
U
∞
x,α
correspond to the instantaneous velocities of particle
α
owing to
the presence of all the other particles and to the external velocity field, respectively.
For the case of simple shear flow, the external velocity field at the centre of the
particle is simply
U
∞
x,α
=
̇
γy
α
. In a simulation where the position of each particle is
known exactly and the particle velocities are already decomposed into an affine (
U
∞
)
and a non-affine (
U
h
) part, each term of equation (4.2) is known exactly and can be
manipulated separately. Foss & Brady (1999) used this decomposition to claim that
the long-time self-diffusion coefficient in the
x
-direction can be calculated from only
the non-affine part of the motion (denoted
X
h
in (4.2)), specifically from :
D
xx
= lim
t
→∞
1
2
d
d
t
X
h
(
t
)
X
h
(
t
)
.
(4.3)
Shear-induced self-diffusion in non-colloidal suspensions
299
The theoretical analysis of Morris & Brady (1996) also indicated that such an
operation was possible and would result in the correct value of the self-diffusivity for
any value of the P
́
eclet number.
A fundamental difference between the shear-induced diffusivity of a non-Brownian
suspension and the motion of a Brownian particle is the finite (and large) correlation
time for the non-affine displacements. In the case of a Brownian particle, there
exists a time interval small enough so that the particle position has not changed
significantly, but also large enough so that the random component of the particle’s
velocity is completely uncorrelated with the random component of the velocity at the
previous time interval; this time interval is given by the particle momentum relaxation,
τ
B
=
m/
6
π
ηa
, which is much smaller than the time scale of configurational changes,
a
2
/D
. Such a time scale does not exist for the case of sheared non-Brownian particles.
The hydrodynamic component of the particle velocity, which is what gives rise
to diffusive-like motion for sheared suspensions, is only a function of the particle
configuration; as a result, in order for this velocity to change significantly (or the
velocity autocorrelation function to approach zero), the particle configuration also
has to change significantly, and for some cases (e.g. small volume fractions) strains up
to 10–20 may be necessary. As was demonstrated in the previous section, this can be
a serious limitation for the calculation of
D
yy
and
D
zz
if long enough intervals are not
sampled; the consequences are far more apparent, however, for the calculation of
D
xx
,
since now we also need to consider the coupling between the affine and non-affine
motions over the long correlation times.
Starting from the displacement (4.2), the mean-square displacements
xx
and
xy
can be readily calculated:
x
(
t
)
x
(
t
)
=
X
h
(
t
)
X
h
(
t
)
+
〈
(
∫
t
0
̇
γy
(
τ
)d
τ
)
2
〉
+2
〈
X
h
(
t
)
∫
t
0
̇
γy
(
τ
)d
τ
〉
,
(4.4)
and
x
(
t
)
y
(
t
)
=
X
h
(
t
)
y
(
t
)
+
〈
y
(
t
)
∫
t
0
̇
γy
(
τ
)d
τ
〉
.
(4.5)
Differentiating with respect to time and using (4.5) we have
d
d
t
x
(
t
)
x
(
t
)
=
d
d
t
X
h
(
t
)
X
h
(
t
)
+2
̇
γ
x
(
t
)
y
(
t
)
+2
〈
d
X
h
(
t
)
d
t
∫
t
0
̇
γy
(
τ
)d
τ
〉
,
(4.6)
and
d
d
t
x
(
t
)
y
(
t
)
=
d
d
t
X
h
(
t
)
y
(
t
)
+
̇
γ
y
(
t
)
y
(
t
)
+
〈
d
y
(
t
)
d
t
∫
t
0
̇
γy
(
τ
)d
τ
〉
.
(4.7)
Although derived for the case of a sheared non-Brownian suspension, the
above equations are also valid for the motion of a Brownian particle, where the
hydrodynamic displacements
X
h
and
y
are now replaced by the equivalent random
Brownian displacements. The fundamental difference between these two cases lies
in the calculation of the last terms in (4.6) and (4.7). For a Brownian particle, the
‘diffusive distance’ a particle has travelled in the
y
-direction, represented by the
integral
∫
t
0
̇
γy
(
τ
)d
τ
, is uncorrelated with the particle’s instantaneous random velocity
(d
y/
d
t
or d
X
h
/
d
t
), since the particle velocities are uncorrelated between time steps
t
(the Brownian relaxation time
τ
B
t
, at least within an accuracy
O
(
t
)). This is not
the case for hydrodynamically interacting non-Brownian particles; now a particle’s
position, or equivalently the distance the particle has travelled in the
y
-direction, is
300
A. Sierou and J. F. Brady
strongly correlated with the instantaneous velocity of the particle, since this velocity is
purely a result of a particle’s position and only changes significantly when the particle
configuration changes significantly. It can be argued that this last term will only
introduce a contribution of the order of the correlation time, but now the correlation
time is large and scales as
̇
γ
−
1
, and the last term in (4.6) or (4.7) becomes constant
for long times, contributing a linear term to the total mean-square displacement.
The calculation of
xx
as given in (4.4) or (4.6) accurately describes the particle
displacements and is all the information needed for the motion of the particles in
the flow direction. The calculation of a self-diffusion coefficient – corresponding to
a ‘diffusive’-like term in a master or Fokker–Planck equation – from this mean-
square displacement still needs to be addressed, however. It has been suggested
(Morris & Brady 1996; Foss & Brady 1999) that the self-diffusion coefficient in the
longitudinal direction can simply be calculated by only taking into account the non-
affine displacement according to (4.3), or simply the first term on the right-hand side
of (4.6). The presence of the last term on the right-hand side of (4.6) or (4.7) makes
such a suggestion suspect, and, as will be demonstrated below, the contribution of this
term is constant at long times. Note that the second term on the right-hand side of
(4.7) is proportional to
̇
γD
yy
t
at long times and therefore gives a
D
xy
(= (d
/
d
t
)
xy
/
2)
growing as
t
and an
xy
growing as
t
2
as in (4.1
b
). Similarly, the
xy
-term on the
right-hand side of (4.6) gives a
t
3
growth for
xx
. These are the well-known (Elrick
1962) convectively enhanced mean-square displacements of a (Brownian) particle and
cause no conceptual difficulties. Note also that the only non-zero cross-term is
xy
.
To address the issue of the proper definition of the ‘diffusivity’ in the direction of
the shearing motion, we write a master equation for the probability density of finding
a particle described by (4.2) (Van Kampen 1992):
∂P
∂t
=
(
1
2
d
d
t
xx
)
∂
2
P
∂x
2
+
(
1
2
d
d
t
yy
)
∂
2
P
∂y
2
+
(
1
2
d
d
t
zz
)
∂
2
P
∂z
2
+2
(
1
2
d
d
t
xy
)
∂
2
P
∂x∂y
,
(4.8)
where
x
=
y
=
z
= 0, and it has been assumed that the derivatives of the second
moments are only functions of time (i.e. they have no spatial dependence). Equation
(4.8) simply follows from a second-order expansion of the probability distribution
of finding a marked particle and is a very general form valid for a large number
of systems. Note that there is no explicit shearing motion, i.e. no advection term in
(4.8), because at the macroscopic scale the marked particle has averaged over the
microscale shearing flow and wanders ‘diffusively’.
The solution of (4.8) can be written as (Van Kampen 1992):
P
=
1
1
/
2
exp
(
−
x
2
2
σ
x
2
−
y
2
2
σ
y
2
−
z
2
2
σ
z
2
−
2
xy
2
σ
xy
)
,
(4.9)
where
=
σ
z
2
(
σ
x
2
σ
y
2
−
σ
2
xy
)
,
(4.10
a
)
σ
x
2
=
B
11
−
B
2
12
/B
22
,
(4.10
b
)
σ
y
2
=
B
22
−
B
2
12
/B
11
,
(4.10
c
)
σ
z
2
=
B
33
,
(4.10
d
)
σ
xy
=
B
12
−
(
B
11
B
22
)
/B
12
,
(4.10
e
)
Shear-induced self-diffusion in non-colloidal suspensions
301
and
B
is the matrix of mean-square displacements
B
11
=
xx
,
(4.11
a
)
B
22
=
yy
,
(4.11
b
)
B
33
=
zz
,
(4.11
c
)
B
12
=
xy
.
(4.11
d
)
It is straightforward to verify that (4.9) is a solution of (4.8). This is still a very general
result since the functional form for the second moments of (4.11) has not yet been
specified. Using the functional forms for d
xx
/
d
t
and d
xy
/
d
t
given by (4.6) and
(4.7), equation (4.8) can be rewritten as:
∂P
∂t
=
−
y
∂P
∂x
+
D
xx
(
t
)
∂
2
P
∂x
2
+
D
yy
(
t
)
∂
2
P
∂y
2
+
D
zz
(
t
)
∂
2
P
∂z
2
+2
D
xy
(
t
)
∂
2
P
∂x∂y
,
(4.12)
where
D
xx
(
t
)=
1
2
d
d
t
X
h
X
h
+
〈
d
X
h
d
t
∫
t
0
̇
γy
(
τ
)d
τ
〉
,
(4.13
a
)
D
yy
(
t
)=
1
2
d
d
t
yy
,
(4.13
b
)
D
zz
(
t
)=
1
2
d
d
t
zz
,
(4.13
c
)
D
xy
(
t
)=
1
2
d
d
t
X
h
y
+
1
2
〈
d
y
d
t
∫
t
0
̇
γy
(
τ
)d
τ
〉
.
(4.13
d
)
Note that the
x
(
t
)
y
(
t
)
and
y
(
t
)
y
(
t
)
terms of (4.6) and (4.7) (multiplying the
second-order partial derivatives of
P
) have been rewritten as the convective flux term
(
−
y∂P/∂x
) with
P
given by (4.9). This is a straightforward calculation:
1
2
∂
2
P
∂x
2
2
xy
+
∂
2
P
∂xy
yy
=
P
(
x
2
σ
x
2
σ
x
2
+
y
2
σ
xy
σ
xy
+
2
xy
σ
x
2
σ
xy
−
1
σ
x
2
)
xy
+
P
(
xy
σ
x
2
σ
y
2
+
x
2
σ
x
2
σ
xy
+
y
2
σ
y
2
σ
xy
+
xy
σ
xy
σ
xy
)
yy
=
P
(
x
2
B
2
22
B
12
(
B
22
B
11
−
B
2
12
)
2
+
y
2
B
2
12
B
12
(
B
22
B
11
−
B
2
12
)
2
−
2
xyB
22
B
2
12
(
B
22
B
11
−
B
2
12
)
2
−
B
22
B
12
(
B
22
B
11
−
B
2
12
)
+
xyB
2
22
B
11
(
B
22
B
11
−
B
2
12
)
2
−
x
2
B
2
22
B
12
(
B
22
B
11
−
B
2
12
)
2
−
y
2
B
11
B
12
B
22
(
B
22
B
11
−
B
2
12
)
2
+
xyB
2
12
B
22
(
B
22
B
11
−
B
2
12
)
2
+
B
22
B
12
(
B
22
B
11
−
B
2
12
)
)
=
P
−
y
2
B
12
+
xyB
22
(
B
22
B
11
−
B
2
12
)
,
and
−
y
∂P
∂x
=
P
(
xy
σ
x
2
+
y
2
σ
xy
)
=
P
−
y
2
B
12
+
xyB
22
(
B
22
B
11
−
B
2
12
)
.
302
A. Sierou and J. F. Brady
We have denoted the diffusivities as functions of time so that (4.12) and (4.13) are
valid at each instant in time, and not only for very long times when the motion
becomes truly diffusive. As has been mentioned, at long times, both contributions to
D
xx
and
D
xy
become constant with time – it is those values of the diffusivity that we
wish to calculate.
From the master equation (4.8) describing the ‘macroscopic’ behaviour, we have
derived the conventional convection–diffusion equation (4.12) where, by construction,
the diffusion coefficients are given by (4.13). It is these diffusion coefficients that we
wish to calculate, and it is now apparent that in order to calculate the diffusivity in the
x
-direction we need more information than simply the non-affine displacements. The
theoretical analysis of Morris & Brady (1996) also agrees with the above conclusion;
unfortunately, an error in their earlier calculation (which has been corrected) led to
the omission of the last terms of (4.13) for
D
xx
and
D
xy
.
For simplicity of notation we rewrite (4.13) as:
D
xx
(
t
)=
D
xx
(
t
)+
D
corr
xx
(
t
)
,
(4.14
a
)
D
xy
(
t
)=
D
xy
(
t
)+
D
corr
xy
(
t
)
,
(4.14
b
)
where
D
xx
(
t
)=
1
2
d
d
t
X
h
X
h
,
(4.15
a
)
D
xy
(
t
)=
1
2
d
d
t
X
h
y
,
(4.15
b
)
and
D
corr
xx
(
t
)=
〈
d
X
h
d
t
∫
t
0
̇
γy
(
τ
)d
τ
〉
,
(4.16
a
)
D
corr
xy
(
t
)=
1
2
〈
d
y
d
t
∫
t
0
̇
γy
(
τ
)d
τ
〉
.
(4.16
b
)
We refer to the term omitted by the previous authors as a ‘correction’ simply for
clarity – both terms can be equally important for the shear-induced diffusivity. It can
also be seen as a correction compared to the Brownian case in the absence of shear
where such a term is absent.
Calculating the self-diffusivities from (4.13) is now straightforward for any
simulation. We chose to express these diffusivities as functions of the non-affine
displacement because we wanted to stress the importance of the missing term in
the calculations of Morris & Brady (1996) and Foss & Brady (1999). It is also
straightforward to calculate the diffusivities directly from (4.6) and (4.7):
D
xx
(
t
)=
1
2
d
d
t
X
h
X
h
+
〈
d
X
h
d
t
∫
t
0
̇
γy
(
τ
)d
τ
〉
=
1
2
d
d
t
x
(
t
)
x
(
t
)
−
̇
γ
x
(
t
)
y
(
t
)
,
(4.17)
and
D
xy
(
t
)=
1
2
d
d
t
X
h
y
+
1
2
〈
d
y
d
t
∫
t
0
̇
γy
(
τ
)d
τ
〉
=
1
2
d
d
t
x
(
t
)
y
(
t
)
−
1
2
̇
γ
y
(
t
)
y
(
t
)
.
(4.18)
Shear-induced self-diffusion in non-colloidal suspensions
303
In other words, all the information needed for the calculation of the diffusivities is
simply the overall mean-square displacements
xx
,
xy
and
yy
. In an experimental
set-up, equations (4.13) or (4.17) and (4.18) can be used without further modification.
The only requirement is that the particle trajectories must be followed closely enough
so that the derivatives and integrals of the mean-square displacements are calculated
with sufficient accuracy. Note also, that in general this would be superior to trying to
extract from an experiment (or simulation) the linear-in-time terms of the advectively
dominated mean-square displacements in (4.1).
Finally, it is important to realize that the time-dependent diffusivities
D
xx
(
t
)and
D
xy
(
t
) defined in (4.13) or (4.17) and (4.18) are not necessarily identical at long
times to what we would obtain by simply looking at the terms in the mean-square
displacements given by (4.1) that grow linearly in time. This can be appreciated by
supposing that
D
yy
(
t
) in (4.13
b
) is independent of time and integrating to find
yy
=2
D
yy
t
+
a
yy
,
where
a
yy
is a constant of integration. The constant
a
yy
would correspond to the
intercept as
t
→
0 of the straight line in figure 6 and is a manifestation that the
motion is not diffusive at short times. Now, using this in (4.18), assuming
D
xy
is a
constant, and integrating gives
xy
=(2
D
xy
+
̇
γa
yy
)
t
+
D
yy
̇
γt
2
+
b
xy
,
where
b
xy
is another constant and so forth for
xx
. Thus, a measurement of the
linear in time dependence of
xy
(assuming this was possible given the domination
of the
D
yy
̇
γt
2
term) would yield an apparent
xy
diffusivity of
D
xy
=
D
xy
+
1
2
̇
γa
yy
.
So, is
D
xy
or
D
xy
the correct diffusivity?
From the perspective of the long-time behaviour of the mean-square displacement,
we would naturally choose
D
xy
. From the perspective of a (possibly time-dependent)
diffusivity in the Fokker–Planck equation (4.12),
D
xy
is correct. In other words, since
the Fokker–Planck equation involves the time-rate of change of the moments (e.g.
(d
yy
/
d
t
)
/
2, which are in general functions of time, only approaching a constant in
the long-time limit), the mean-square displacements at long time given by the solution
of the Fokker–Planck equation are not equal to
D
xy
etc., but to (e.g.)
D
xy
,andso
the two perspectives are actually in complete agreement. This apparent discrepancy
illustrates the importance of distinguishing between a diffusivity,
D
xy
, and a mean-
square-displacement linear in time,
D
xy
+
̇
γa
yy
/
2, when advective motion couples to
diffusion. We take the definition of the diffusivity (time-dependent or not) as that
which appears naturally in the Fokker–Planck equation, with the understanding that
we must solve this equation to give the resulting mean-square displacements and
they may be different from the underlying diffusivities when advective effects are
operative.
Figure 10 shows the
D
xx
self-diffusivity (as calculated by either (4.13) or (4.17) –
in a simulation where the exact positions are known at each instant in time, the two
approaches are identical) as a function of time for a system of
N
= 1000 particles
and a volume fraction of 0
.
35. The curve has been averaged over all particles and
over 30 independent configurations. The two separate contributions to the diffusivity,
(4.14
a
), are also shown. Both contributions become constant after the correlation time,
304
A. Sierou and J. F. Brady
0.2
0.1
0
–0.1
12
10
8
6
4
2
0
γ
t
D
xx
D
*
xx
D
xx
corr
.
γ
a
2
.
D
xx
Figure 10.
The self-diffusivity in the longitudinal direction,
D
xx
, is plotted as a function of
strain for a system of
N
= 1000,
φ
=0
.
35. The values of the two terms contributing to
D
xx
according to (4.14
a
) are also given.
–0.10
–0.08
–0.06
–0.04
–0.02
0
0.02
12
10
8
6
4
2
0
γ
t
D
D
D
xy
xy
xy
corr
*
.
γ
.
a
2
D
xy
Figure 11.
The cross-diffusivity,
D
xy
, is plotted as a function of strain for a system of
N
= 1000,
φ
=0
.
35. The values of the two terms contributing to
D
xy
according to (4.14
b
)are
also given.
verifying that they both should be included in the calculation of the self-diffusivity.
The same behaviour is seen in figure 11, this time for
D
xy
.Bothcorrectionsto
D
xx
and
D
xy
are negative, resulting in a smaller absolute value of the
D
xx
diffusivity and
a larger absolute value of the
D
xy
diffusivity, since the non-affine term,
D
xy
, is also
negative (the physical origin of the negative
D
xy
is discussed in Foss & Brady 1999).
For this relatively high-volume fraction, the effect of the correction to the value of
D
xx
is small. The effect on the value of
D
xy
is rather significant, however, since the
non-affine contribution,
D
xy
, is relatively small. The correction to the diffusivities for
a much lower-volume fraction of
φ
=0
.
15 is shown in figure 12. The absolute values
of both corrections are now significantly larger, in agreement with the fact that it now
Shear-induced self-diffusion in non-colloidal suspensions
305
–0.20
–0.16
–0.12
–0.08
–0.04
0
0.04
D
corr
20
15
10
5
0
γ
t
φ
= 0.15,
N
= 1000
correction for
D
xy
correction for
D
xx
.
γ
a
2
.
Figure 12.
The corrective terms,
D
corr
xx
and
D
corr
xy
,forasystemof
N
= 1000 and a relatively
small-volume fraction
φ
=0
.
15 are plotted as a function of strain.
–0.12
–0.10
–0.08
–0.06
–0.04
–0.02
0
D
corr
0.6
0.5
0.4
0.3
0.2
0.1
0
φ
D
xx
D
xy
.
γ
a
2
corr
corr
Figure 13.
The corrective terms,
D
corr
xx
and
D
corr
xy
,forasystemof
N
= 1000 are plotted as a
function of the volume fraction.
takes a much longer time for the particle velocities to become uncorrelated, increasing
the time limit for the integrals in
D
corr
xx
and
D
corr
xy
; it is also apparent that it takes a
longer time for the corrective term to reach a constant value.
For the two-particle problem the non-affine term of the mean-square displacement
is well behaved and gives a finite contribution to the diffusivity as
φ
→
0. The
contribution
D
corr
, which leads to the singular behaviour observed by Acrivos
et al.
(1992), can be very important, especially for small-volume fractions, and can
dominate the overall behaviour of the mean-square displacement. Figure 13 shows the
corrections to the self-diffusivities as a function of volume fraction. Both corrections
are increasing (in magnitude) functions of the volume fraction for small volume
306
A. Sierou and J. F. Brady
fractions, and decreasing (in magnitude) functions of
φ
as the volume fraction is
further increased. This behaviour can be explained from the functional form (4.16).
The correlation time is a monotonically decreasing function of the volume fraction
since it depends directly on the number of collisions between particles. Therefore, the
time over which the integrals of (4.16) are calculated decreases as
φ
increases. On the
other hand, for low-volume fractions, the particle’s non-affine
y
displacement, upon
which the corrections to the diffusivity directly depend, are smaller –
D
yy
is small for
small-volume fractions. The combination of those two factors results in the observed
non-monotonic behaviour for the corrective contribution to the self-diffusivities.
4.2.
The volume fraction dependence
Before presenting the results for
D
xx
, we shall briefly discuss some limitations
concerning the calculation of the diffusivity in the flow-direction for small volume
fractions. It has been shown by Acrivos
et al.
(1992) that in the two-particle limit,
the behaviour of the mean-square displacement in the
x
-direction is singular and
more particles are needed to remove this singularity. Even when the singularity is
not present, however, as is the case in our simulations (the singularity is a result of
particles coming from infinity interacting with a test particle – in a periodic system, the
finite size of the simulation cell introduces a cutoff), the resulting
x
-displacement for
a given two-body interaction can be large. For the diffusive regime to be established,
each particle must sample a randomly distributed number of different displacements
with a zero average. In a dilute periodic system, this implies that each particle, after it
interacts with a given second particle and before it again interacts with this particle’s
image, must sample a sufficient number of interactions with other, third particles.
This can give a strict condition on the minimum size of the unit cell, and thus on
the size of the system simulated. Although the same condition must be satisfied for
concentrated systems, the magnitude of each individual displacement is significantly
smaller and their frequency significantly larger at high
φ
.
To demonstrate the severity of this limitation, consider two particles suspended in
a periodic unit cell of finite length
L
, initially at positions with relative displacements
X
0
,
Y
0
, and for simplicity
Z
0
= 0 (as shown in figure 14). If initially in positions
with different
y
-coordinates, the two particles will interact and at the end of the
interaction will return at the same
y
-positions, but there will be a net displacement
in the
x
-direction (Acrivos
et al.
1992). After the end of the interaction, the particles
will continue moving undisturbed until the ‘faster’ particle reaches the end of the
periodic domain, at which time it will re-enter the cell on the opposite face. For a
periodic domain, the maximum distance between two particles has to be less than
L
, so eventually (the time will depend both on
L
and the difference in the affine
undisturbed velocities between the two particles, which is proportional to
Y
0
)the
two particles will interact again. This second interaction will be identical to the first,
since the particles’ vertical positions did not change as a result of the first interaction.
This behaviour will continue indefinitely and at each ‘collision’ the displacement of
each particle is always exactly the same. It is clear that such a scenario would never
evolve into diffusive-like motion, even at infinite time, since only one displacement is
sampled for each particle, and of course
X
=0.
In a dilute system the same problem is present: a given particle, after interacting
with an initial second particle, has a very high probability that it will not collide
with any third particle until it exits the periodic cell and interacts with the same
particle for the second time. Such a behaviour leads to particle velocities that are
strongly correlated for aphysically long times and mean-square displacements that