Complex fluids are fluids with a suspended microstructure sucha sp o l y m e r s with a long chain structure. A simplified model is that of suspensions of rod-like particles. In such systems multiscale interactions can cause interesting phenomena. The macroscopic flow leads to a change of the orientation and shape (in the case of flexible particles) of the suspended microstructure. In turn, this produces a fluid stress which interacts with and modifies the macroscopic flow. Here we consider the sedimentation of dilute suspensions of rod-like particles under gravity. This problem exhibits interesting pattern formations and has been studied by several authors in theoretical, numerical and experimental work. Understanding sedimentation is of importance in recycling processes

gayoldMechanics

Feb 21, 2014 (3 years and 6 months ago)

121 views

A KINETIC MODEL FOR THE SEDIMENTATION OF ROD–LIKE PARTICLES
CHRISTIANE HELZEL
!
,FELIX OTTO

,
AND
ATHANASIOS E.TZAVARAS

Abstract.
We consider a kinetic multi-scale model,describing suspensions of rod-lik
e particles,which couples a
microscopic Smoluchowski equation - for the distribution of rod or
ientations - to a macroscopic Stokes or Navier-Stokes
equation via an elastic stress tensor and a gravitational forcing term.
A reciprocal coupling of phenomena on these
di!erent scales leads to the formation of clusters.The buoyancy force creat
es a macroscopic velocity gradient that causes
the microscopic particles to align in such a way that their sedimentatio
n reinforces the formation of clusters of higher
particle density.
We perform a linear stability analysis and show that linear theory can
predict a wavelength selection mechanism for
the cluster width,provided that the Reynolds number is larger than zer
o.An asymptotic analysis of the largest eigenvalue
around
Re
= 0 explains this response,while numerical simulations of the nonlin
ear model confirmthe wavelength selection.
By looking at special flow configurations and at parameter scalings,we der
ive simpler models which describe the aggregate
behavior of the suspension in various regimes.In particular,in a di!u
sive type of scaling,we obtain the two-dimensional
Keller-Segel model.
Key words.
sedimentation,rod-like particles,linear stability,asymptotic exp
ansions
AMS subject classifications.
76T20,74A25,35B40
1.Introduction.
Complex fluids are fluids with a suspended microstructure suc
h as polymers
with a long chain structure.A simplified model is that of susp
ensions of rod-like particles.In such
systems multiscale interactions can cause interesting phe
nomena.The macroscopic flow leads to a
change of the orientation and shape (in the case of flexible pa
rticles) of the suspended microstructure.
In turn,this produces a fluid stress which interacts with and
modifies the macroscopic flow.
Here we consider the sedimentation of dilute suspensions of
rod-like particles under gravity.This
problem exhibits interesting pattern formations and has be
en studied by several authors in theoretical,
numerical and experimental work.Understanding sedimenta
tion is of importance in recycling processes
of the paper industry [7],and related models arise also in th
e description of self-locomotion [20].In
contrast to spherical particles,rod-like particles have a
n orientation,what adds complexity to the
sedimentation process.For instance,the sedimentation ve
locity of a single particle depends on its
orientation.Furthermore,a particle falls not only in the g
ravity direction but also sidewards.By
restricting considerations to rigid rods we eliminate defo
rmation e!ects such as bending of polymers.
While previous experimental and numerical work studies the
sedimentation of non-Brownian particles,
we mainly concentrate on the case of Brownian particles.
The kinetic model considered here is an extension of the Doi m
odel for dilute suspensions of rod-
like particles described in Doi and Edwards [5].In previous
work we considered the Doi model for
suspensions of rod-like molecules neglecting the e!ects of
gravity,see [9],[16].
The sedimentation of fibers has been studied in experimental
work for instance by Guazzelli and
coworkers [10,11,14].These experiments reveal the follow
ing scenarios:Starting from a well-stirred
suspension packets of particles form after some time.These
packets seem to have a mesoscopic equilib-
rium width,suggesting that the density of particles acquir
es variations of a characteristic length scale.
Within such a cluster,individual particles are aligned wit
h the direction of gravity during most of the
time;occasionally they flip.The average settling speed in a
suspension is larger than the sedimentation
speed of a single particle oriented in the direction of gravi
ty.
In recent years the sedimentation of rod-like particles (an
d other orientable particles) has also been
studied via numerical simulations of multi-scale models.G
ustavsson and Tornberg [8,21] used a very
detailed description of rod-like particles in a dilute susp
ension based on a slender body approximation.
They were able to simulate suspensions with up to a few hundre
d particles and a domain size of the
order of a few particle length.Butler and Shaqfeh [2] used a l
ower order slender body description.
Saintillan et al.[17] accelerated this algorithm using fas
t summation techniques.This allowed them to
!
Department of Mathematics,Ruhr-University Bochum,Germany
christiane.helzel@rub.de
,

Institute of Applied Mathematics,University of Bonn,Germany
otto@iam.uni-bonn.de

Department of Applied Mathematics,University of Crete,Heraklion,Greece
and Institute for Applied and Computational Mathematics,FORTH
,Heraklion,Greece,
tzavaras@tem.uoc.gr
1
2
simulate several thousand particles.Wang and Layton [22] u
sed the immersed boundary method for
their two-dimensional numerical studies.All numerical re
sults confirm the basic experimental findings:
Packet formation and alignment in gravity direction (at lea
st for higher concentration).Note that
the models used in these simulations are of more microscopic
nature than the model considered here.
Instead of a number density function for the rod orientation
in every point of the domain,those authors
consider a large number of individual particles.Moreover,
Brownian e!ects are not considered,which
amounts to setting the absolute temperature
!
= 0 in our model.
Koch and Shaqfeh [12] performed a linear stability analysis
for a related kinetic model by neglecting
Brownian e!ects.They found a purely imaginary continuum sp
ectrum,but also a point spectrum with
positive real part.This predicts instability of density mo
dulations for horizontal waves with a su"ciently
small wavenumber but not a wavelength selection mechanism.
The introduction of Brownian e!ects in
terms of translational di!usion does also not provide a wave
length selection mechanism at the level of
linear stability analysis,see Saintillan et al.[19].
In the case of non-Brownian sedimentation of spherical part
icles,there is a mathematically related
divergence of fluctuations at long wavelength [3].For that p
hysical system,a small stable density
stratification has been proposed as a mechanism which suppre
sses this long wavelength divergence [13].
In the case of sedimenting rod-like particles it also leads t
o a wavelength selection at the level of linear
stability analysis [18].In experiments and numerical simu
lations,this stratification may be induced by
the boundary conditions at the horizontal container walls [
18].However,the predicted scaling of the
unstable wavelength does not agree with experiments [15].
Our objective in this work is on one hand to address the issue o
f initial wavelength selection,and
on the other to discuss the subsequent mechanism of aggregat
ion.We employ a kinetic model in the
spirit of Doi theory focusing on dilute suspensions of rod-l
ike particles.We first present an analysis of
the thermodynamic consistency of the kinetic model and purs
ue a non-dimensional analysis.There are
three non-dimensional parameters in the problem.The Reyno
lds number which quantifies the relative
importance of inertial forces versus viscous forces and two
additional non-dimensional parameters which
measure the importance of the rate of work due to elastic forc
e vs.buoyancy as well as buoyancy vs.
viscous force.The influence of these parameters on linear st
ability theory is studied.In Section 3
we present a linear stability analysis which includes also t
he e!ects of rotational di!usion.We show:
If the macroscopic flow is modeled by the Stokes equation,lin
ear stability theory does not predict a
wavelength selection mechanism.However,if the macroscop
ic flow is modeled by the Navier-Stokes
equation (i.e.
Re >
0),we obtain a wavelength selection mechanism.This behavi
or is explained by
doing a (singular) perturbation analysis of the largest eig
envalue in Reynolds number.
In Section 4,we perform the stability analysis for the sedim
entation of non-Brownian particles.
While we do not have a theoretical justification for the use of
the linear model in the case of non-
Brownian particles,we show that the predictions of the line
ar model are in good agreement with
numerical simulations of the nonlinear model.
Finally,in Section 5 we derive the aggregate behavior of the
suspension in various parameter
scalings.We show that for a particular flow configuration and
a di!usive scaling of the equations the
density follows a nonlinear concentration mechanismthat i
s captured by the well-known two-dimensional
Keller-Segel model.Related derivations of macroscopic mo
dels as limits of kinetic models can be found
in [4,6]
We also present numerical simulations for the coupled nonli
near multi-scale model which indicate
the initial wavelength selection as well as capture the subs
equent concentration.
2.The mathematical model.
We describe a kinetic model for sedimentation in dilute susp
en-
sions of rod-like Brownian particles.Models of this type we
re introduced by Doi and Edwards,see [5,
Ch.8].In [16] a related model for suspensions of rod-like pa
rticles was considered.Here we extend this
model to account for the e!ects of gravity.
We consider inflexible rod-like particles of thickness
b
which is much smaller than the particles
length
l
.Our considerations are restricted to the dilute regime whi
ch is characterized by the relation
"
!
l
!
3
where
"
is the number density.Note that in contrast to the model cons
idered in [16],the
number density is not constant here.The orientation of a rod
-like particle is characterized by
n
"
S
d
!
1
where
d
is the dimension of the macroscopic physical space#
#
R
d
.
3
g
n
u
x
1
x
2
x
3
x
Fig.2.1
.
Basic notation for rod-like molecule which is falling sidewar
ds.
We denote by
e
2
the unit vector in the upward direction and by
g
=
$
g
e
2
the acceleration of
gravity with gravitational constant
g
.Furthermore,let
m
0
denote the mass of an individual rod-like
particle and
G
=
$
m
0
g
e
2
be the force of gravity on a single particle.Some of our basic
notation is
depicted in Figure 2.1.
In a quiescent suspension (i.e.macroscopic velocity
u
= 0) each particle sediments at a speed
depending on its orientation
n
according to
d
x
dt
=
1
#
||
(
G

n
)
n
+
1
#
"
(
G
$
(
G

n
)
n
)
=
!
1
#
||
n
%
n
+
1
#
"
(
I
$
n
%
n
)
"
G
=
1
#
"
(
n
%
n
+
I
)
G
where
#
||
and
#
"
are the frictional coe"cients in the tangential and the norm
al direction which satisfy
#
"
= 2
#
||
.In particular,a particle with a horizontal orientation se
diments slower than a particle with
a vertical orientation and a particle of oblique orientatio
n moves also sideways.
In addition a macroscopic velocity gradient causes a rotati
on according to the equation
d
n
dt
=
P
n
!
&
x
un
where
P
n
!
&
x
un
:=
&
x
un
$
(
&
x
un

n
)
n
is the projection of the vector
&
x
un
on the tangent space
at
n
.
We next include the e!ects of rotational and translational B
rownian motion and account for the
macroscopic mean flow
u
(
x
,t
).The model then becomes the system of stochastic di!erenti
al equations
d
x
=
u
dt
+
!
1
#
||
n
%
n
+
1
#
"
(
I
$
n
%
n
)
"
G
dt
+
#
2
k
B
!
#
||
n
%
n
+
2
k
B
!
#
"
(
I
$
n
%
n
)
dW
d
n
=
P
n
!
&
x
un
dt
+
#
2
k
B
!
#
r
dB
(2.1)
where
W
is the translational Brownian motion and
B
is the rotational Brownian motion,
#
r
is the
rotational friction coe"cient,
k
B
is the Boltzmann constant and
!
the absolute temperature.
The above stochastic equations may be equivalently express
ed via the Smoluchowski equation for
4
the evolution of the
local orientational distribution function
:
$
t
f
+
&
x

$!
u
+
1
#
"
(
n
%
n
+
I
) (
$
m
0
g
e
2
)
"
f
%
+
&
n

(
P
n
!
&
x
un
f
)
=
k
B
!
#
r
$
n
f
+
k
B
!
#
"
&
x

(
n
%
n
+
I
)
&
x
f.
(2.2)
Here
f
(
x
,t,
n
)
d
n
describes the number of particles per unit volume at macrosc
opic position
x
and
time
t
with orientations in the element centered at
n
and of volume
d
n
.The second term on the left
hand side of (2.2) models transport of the center of mass of th
e particles due to the macroscopic flow
velocity and due to gravity.The last term on the left hand sid
e models the rotation of the axis due
to a macroscopic velocity gradient
&
x
u
.The terms on the right hand side describe rotational as well
as translational di!usion.The gradient,divergence and La
placian on the sphere are denoted by
&
n
,
&
n

and $
n
,while the gradient and divergence in the macroscopic flow do
main are denoted by
&
x
and
&
x

.The total number of rod-like particles is
&
!
&
S
d
"
1
f
(
x
,t,
n
)
d
n
d
x
=
&
!
&
S
d
"
1
f
(
x
,
0
,
n
)
d
n
d
x
=
N,
i.e.
f
has dimensions of number density.It is convenient to rewrit
e the Smoluchowski equation in the
form
$
t
f
+
&
x

(
u
f
) +
&
n

(
P
n
!
&
x
un
f
)
=
D
r
$
n
f
+
D
"
&
x

(
I
+
n
%
n
)
'
&
x
f
+
1
k
B
!
f
&
x
U
(
,
(2.3)
where
D
r
:=
k
B
!
"
r
and
D
"
:=
k
B
!
"
!
stand for the rotational and translational di!usion coe"ci
ents and
U
(
x
) =
m
0
g
(
x

e
2
) is the potential of the gravity force
G
=
$&
U
=
$
m
0
g
e
2
.
As can be seen from (2.3),a velocity gradient
&
x
u
distorts an isotropic distribution
f
which leads
to an increase in entropy.Thermodynamic consistency requi
res (compare with Section 2.1),that this
is balanced by a stress tensor
%
(
x
,t
) given by
%
(
x
,t
):=
k
B
!
&
S
d
"
1
(
d
n
%
n
$
I
)
f
(
x
,t,
n
)
d
n
.
(2.4)
Local variations in the density
m
0
)
S
d
"
1
fd
n
lead to spatial variations of the specific weight of the
suspension that generally can not be compensated by a hydros
tatic pressure and thus trigger a fluid
motion (buoyancy).The macroscopic flow is described by the N
avier-Stokes equation.Let
&
f
be the
density of the fluid which is assumed to be constant.The balan
ce laws of mass and momentum have
the form
&
f
(
$
t
u
+(
u
∙ &
x
)
u
) =
µ
'
x
u
$&
x
p
+
&
x

%
$
&
f
g
e
2
$
'
&
S
d
"
1
fd
n
(
m
0
g
e
2
&
x

u
= 0
(2.5)
The term
&
f
g
e
2
can be incorporated to the pressure.For the linear stabilit
y analysis it will be
convenient to express the momentum equation in the equivale
nt form
&
f
(
$
t
u
+(
u
∙ &
x
)
u
) =
µ
$
x
u
$&
x
p
#
+
&
x

%
+
!
N
V
$
&
S
d
"
1
f d
n
"
m
0
g
e
2
&
x

u
= 0
(2.6)
by redefining the pressure,
p
#
=
p
+
&
f
g
e
2

x
+
m
0
Ng
V
e
2

x
,
5
to account for the hydrostatic pressures,where
V
is the volume occupied by the suspension and
N
the
total number of rod-like particles.
We summarize the final model:
$
t
f
=
$&
x

(
u
f
)
$&
n

(
P
n
!
&
x
un
f
) +
D
r
$
n
f
+
D
"
&
x

(
I
+
n
%
n
)
!
&
x
f
+
1
k
B
!
m
0
g
e
2
f
"
(2.7)
%
(
x
,t
) =
k
B
!
&
S
d
"
1
(
d
n
%
n
$
I
)
f
(
x
,t,
n
)
d
n
(2.8)
&
x

u
= 0
(2.9)
&
f
(
$
t
u
+(
u
∙ &
x
)
u
) =
µ
$
x
u
$&
x
p
+
&
x

%
$
'
&
S
d
"
1
f d
n
(
m
0
g
e
2
(2.10)
2.1.Thermodynamic consistency of the model.
To show thermodynamic consistency of the
model we use the free energy functional
A
[
f
]:=
&
!
&
S
2
(
k
B
!f
ln
f
+
fU
(
x
))
d
n
d
x
,
(2.11)
where
U
(
x
) =
m
0
g
xe
2
is the gravitational potential.
Proposition 2.1.
For
f
satisfying the Smoluchowski equation
(2.7)
,the free energy
A
[
f
]
defined
in (2.11) satisfies the identity
$
t
A
[
f
] +
D
r
k
B
!
&
!
&
S
2
f
|&
n
ln
f
|
2
d
n
d
x
+
D
"
k
B
!
&
!
&
S
2
&
x
*
ln
f
+
1
k
B
!
U
+

*
I
+
n
%
n
+
&
x
*
ln
f
+
1
k
B
!
U
+
d
n
d
x
=
&
!
&
x
u
:
%d
x
+
&
!
m
0
g
e
2
!
&
S
2
fd
n
"

u
d
x
(2.12)
Moreover,the total energy
E
[
u
,f
] =
&
!
!
1
2
&
f
|
u
|
2
+
&
S
2
'
(
k
B
!
)
f
ln
f
+
fU
(
x
)
(
d
n
"
d
x
(2.13)
of the system
(2.7)
-
(2.10)
dissipates.
Proof:
We di!erentiate (2.11) with respect to
t
,
$
t
A
[
f
] =
&
!
&
S
2
'
k
B
!
(1 +ln
f
) +
U
(
x
)
(
f
t
d
n
d
x
(2.14)
and use (2.7) to express the various contributions.
The contribution of the transport term
$&
x

(
u
f
) gives
I
tr
=
$
&
!
&
S
2
'
k
B
!
(1 +ln
f
) +
U
(
x
)
(
&
x

(
u
f
)
d
n
d
x
=
&
!
&
S
2
'
k
B
!
&
x
f
+
f
&
x
U
(

u
d
n
d
x
(2
.
9)
=
&
!
m
0
g
e
2
!
&
S
2
fd
n
"

u
d
x
6
Next,consider the contribution of the drift term
$&
n

(
P
n
!
&
x
un
f
).We obtain:
I
dr
=
$
&
!
&
S
2
'
k
B
!
(1 +ln
f
) +
U
(
x
)
(
&
n

(
P
n
!
&
x
un
f
)
d
n
d
x
(
B.
1)
=
&
!
&
S
2
k
B
!
&
n
ln
f

P
n
!
(
&
x
un
f
)
d
n
d
x
=
&
!
&
x
u
:
k
B
!
&
S
2
n
%&
n
f d
n
d
x
(
B.
3)
,
(2.8)
=
&
!
&
x
u
:
% d
x
.
The contribution of rotational di!usion is:
I
rd
=
&
!
&
S
2
'
k
B
!
(1 +ln
f
) +
U
(
x
)
(
D
r
$
n
f d
n
d
x
=
&
!
&
S
2
'
k
B
!
(1 +ln
f
) +
U
(
x
)
(
D
r
&
n

(
f
&
n
ln
f
)
d
n
d
x
=
$
(
k
B
!
)
2
#
r
&
!
&
S
2
|&
n
ln
f
|
2
f d
n
d
x
Finally,the contribution of the last term in (2.7),modelin
g translational friction and translational
di!usion,reads
I
tdf
=
&
!
&
S
2
'
k
B
!
(1 +ln
f
) +
U
(
x
)
(
D
"
&
x

(
I
+
n
%
n
)
'
&
x
f
+
f
k
B
!
&
x
U
(
=
$
(
k
B
!
)
2
#
"
&
S
2
&
!
&
x
!
ln
f
+
1
k
B
!
U
"

(
I
+
n
%
n
)
f
&
x
!
ln
f
+
1
k
B
!
U
"
Combining all these contributions together yields (2.12).
Next,we multiply the Navier-Stokes equation (2.10) by
u
and integrate over#to obtain the balance
of the kinetic energy
d
dt
&
!
1
2
&
f
|
u
|
2
d
x
+
µ
&
!
&
x
u
:
&
x
u
d
x
=
$
&
!
&
x
u
:
% d
x
$
&
!
u

m
0
g
e
2
!
&
S
2
fd
n
"
d
x
.
(2.15)
Combining (2.12) and (2.15) leads to the balance of total ene
rgy.In particular,it follows that the total
energy dissipates.
2.2.Non-dimensionalization.
We first list the dimensions of the terms that appear in the equ
a-
tions.The units of mass,length and time are denoted by
M
,
L
and
T
.We also monitor the dependence
on the number of particles
N
.

v
:velocity
,
L
T
-

m
0
g
:mass
(
acceleration
,
ML
T
2
-
;

k
B
!
:energy = force
(
length
.
ML
2
T
2
/

#
"
:translational friction orthogonal to rod = force/velocit
y
,
M
T
-
D
"
=
k
B
!
"
!
.
L
2
T
/

#
r
:rotational friction = torque/rotatational velocity
.
ML
2
T
/
D
r
=
k
B
!
"
r
,
1
T
-

µ
=
force/area
velocity gradient
,
M
LT
-

p
:pressure = force/area
,
M
LT
2
-
7

f
:number of particles/volume
,
N
L
3
-

%
)
(
k
B
!
)
f
,
MN
LT
2
-
Now we consider a change of scale of the form
t
=
T
ˆ
t,
x
=
X
ˆ
x
,
u
=
X
T
ˆ
u
,f
=
N
V
ˆ
f,p
=
µ
T
ˆ
p,%
= (
k
B
!
)
N
V
ˆ
%.
(2.16)
We have used two length scales in (2.16):a length scale
X
that will be selected in the course of the
non-dimensionalization process,and a length scale
L
standing for the size of the macroscopic domain
and entering only through the volume occupied by the suspens
ion
V
=
O
(
L
3
).
In these new units the Smoluchowski equation takes the form
$
ˆ
t
ˆ
f
+
&
ˆ
x

'
ˆ
u
ˆ
f
(
+
&
n

'
P
n
!
&
ˆ
x
ˆ
un
ˆ
f
(
=
TD
r
$
n
ˆ
f
+
T
X
2
D
"
&
ˆ
x

(
I
+
n
%
n
)
&
ˆ
x
ˆ
f
+
D
"
X
T
m
0
g
k
B
!
&
ˆ
x

'
(
I
+
n
%
n
)
e
2
ˆ
f
(
(2.17)
We select the time scale
T
to be the time scale of the rotational Brownian motion,i.e.w
e set
T
=
1
D
r
.
(2.18)
Furthermore,we select the length scale
X
via
X
=
m
0
g
#
"
D
r
.
(2.19)
This selection corresponds to a velocity scale
X
T
=
m
0
g
#
"
,
(2.20)
i.e.the velocity scale is proportional to the motion of a sin
gle rod falling due to gravity in a friction
dominated flow.
Using (2.18)-(2.20) in Equation (2.17) we obtain the kineti
c equation in the dimensionless form
$
ˆ
t
f
+
&
ˆ
x

(
ˆ
u
f
) +
&
ˆ
n

(
P
n
!
&
x
ˆ
un
f
)
e
2
f
$&
ˆ
x

((
I
+
n
%
n
)
e
2
f
)
= $
n
f
+
'
&
ˆ
x

(
I
+
n
%
n
)
&
ˆ
x
f,
(2.21)
with the dimensionless parameter
'
:=
k
B
!#
"
D
r
(
m
0
g
)
2
.
(2.22)
Remark 2.1.
The parameter
'
measures the relation of the work of elastic forces vs.the wo
rk of
buoyancy.
Recall the equation (2.15) for the rate of change of the kinet
ic energy of the fluid
$
t
&
&
f
1
2
|
u
|
2
d
x
=
$
&
µ
&
x
u
:
&
x
u
d
x
0
12
3
rate of energy dissipation
$
&
!
&
u
:
%
+
!
&
fd
n
"
m
0
g
u

e
2
"
d
x
0
12
3
rate of work of microscopic forces
.
(2.23)
The rate of work of microscopic forces is the sum of the rate of
work of elastic forces and the rate of
work of buoyancy.Thus we obtain:
rate of work of elastic forces
rate of work of buoyancy
=
4
4
4
4
4
&
x
u
:
%
*
)
fd
n
+
m
0
g
u

e
2
4
4
4
4
4
)
1
T
(
k
B
!
)
|
f
|
X
T
m
0
g
|
f
|
=
k
B
!#
"
D
r
(
m
0
g
)
2
=
'
8
The non-dimensionalization of the elastic stress tensor le
ads to the expression
ˆ
%
=
&
S
d
"
1
(
d
n
%
n
$
I
)
ˆ
fd
n
.
(2.24)
Finally,we non-dimensionalize the Navier-Stokes equatio
n.We obtain
X
2

&
f
(
$
ˆ
t
ˆ
u
+(
ˆ
u
∙ &
ˆ
x
)
ˆ
u
)
= $
ˆ
x
ˆ
u
$&
ˆ
x
ˆ
p
+
XT
µ
k
B
!
X
N
V
&
ˆ
x

ˆ
%
$
N
V
m
0
g
XT
µ
!
&
S
d
"
1
ˆ
f d
n
"
e
2
.
There are three non-dimensional numbers emerging:First,a
Reynolds number based on the microscopic
velocity
X/T
and the microscopic length scale
X
,
Re
:=
&
f
X
2

=
&
f
µ
!
m
0
g
#
"
"
2
1
D
r
.
Second the non-dimensional parameter
(
:=
N
V
m
0
g
XT
µ
which is proportional to the number density
"
=
N
V
and turns out to describe the ratio between buoyancy
and viscous forces.The third non-dimensional quantity is a
composite,expressed in the form
XT
µ
k
B
!
X
N
V
=
!
N
V
m
0
g
XT
µ
"
k
B
!
m
0
gX
=
(
k
B
!#
"
D
r
(
m
0
g
)
2
=
('
The non-dimensional form of the Navier-Stokes equation is
Re
(
$
ˆ
t
ˆ
u
+(
ˆ
u
∙ &
ˆ
x
)
ˆ
u
) = $
ˆ
x
ˆ
u
$&
ˆ
x
ˆ
p
+
('
&
ˆ
x

%
$
(
!
&
S
d
"
1
ˆ
f d
n
"
e
2
.
(2.25)
Remark 2.2.
The parameter
(
measures the relation between the rate of work of buoyancy vs
.the
rate of work of viscous force.
From (2.23) we express
rate of work of buoyancy
rate of work of viscous force
=
4
4
4
4
4
*
)
fd
n
+
m
0
g
ue
2
µ
|&
u
|
2
4
4
4
4
4
)
N
V
m
0
g
X
T
µ
1
T
2
=
N
V
m
0
g
XT
µ
=
(.
(2.26)
We summarize the non-dimensional form of the equations (dro
pping the hats)
$
t
f
+
&
x

(
u
f
) +
&
n

(
P
n
!
&
x
un
f
)
$&
x

((
I
+
n
%
n
)
e
2
f
)
= $
n
f
+
'
&
x

(
I
+
n
%
n
)
&
x
f
%
=
&
S
d
"
1
(
d
n
%
n
$
I
)
f d
n
Re
(
$
t
u
+(
u
∙ &
x
)
u
) = $
x
u
$&
x
p
+
('
&
x

%
$
(
!
&
S
d
"
1
f d
n
"
e
2
&
x

u
= 0
(2.27)
If we express the Navier-Stokes equation in the equivalent f
orm (2.6),then the non-dimensionalized
form is given by
Re
(
$
t
u
+(
u
∙ &
x
)
u
) = $
x
u
$&
x
p
+
('
&
x

%
+
(
!
1
$
&
S
d
"
1
f d
n
"
e
2
&
x

u
= 0
(2.28)
9
2.3.Multi-scale mechanism for instability and cluster form
ation.
The multi-scale mecha-
nismthat leads to the instability and the formation of clust
ers was first explained by Koch and Shaqfeh,
see [12].
In our kinetic model,the function
&
(
x
,t
) =
&
S
d
!
1
f
(
x
,t,
n
)
d
n
measures the density of rod like particles.By virtue of the b
uoyancy term in the Stokes equation,a
density modulation (as indicated in Figure 2.2 (a)) trigger
s a modulated shear flow with flow direction
e
2
(see Fig.2.2 (b)).
By virtue of the microscopic drift term on the sphere
&
n

(
P
n
!
&
x
un
f
) this shear destroys the
uniform distribution
f
in
n
.For moderate local Deborah numbers the distribution
f
slightly concen-
trates in a direction at 45 degrees between the flow direction
and the shear direction as shown in Fig.
2.2 (c).For larger shear rates the distribution
f
concentrates more pronounced in direction of gravity.
In this figure we plotted the average microscopic orientatio
n
)
r
,where
)
is the largest eigenvalue of
%
and
r
is the corresponding eigenvector.
By virtue of the term
$&
x

((
I
+
n
%
n
)
e
2
f
) this nonuniform distribution in
n
implies that
particles on average fall in a direction which is at a non-van
ishing angle between flow direction and
shear direction.Hence this term acts as a horizontal drift t
erm for the modulated density
&
.In fact,
it reinforces the original horizontal modulation of the den
sity
&
since particles with an orientation as
shown in Fig.2.2 (c) move towards the center.
(a)
(b)
(c)
(d)
Fig.2.2
.
Illustration of the concentration mechanism(a) initial density
modulation,(b) velocity field,(c) microscopic
orientation:plot of
!
r
,where
!
is the largest eigenvalue of
"
and
r
is the corresponding eigenvector,(d) increased density
modulation at later time.
3.Linear stability analysis.
Experimental studies for the sedimentation of suspensions
with
rod-like particles (see [10,11,14]) reveal the formations
of packets of particles which seem to have a
mesoscopic equilibrium width.In our formulation of the pro
blem,this means that the density
&
should
acquire variations of a characteristic length scale.Our go
al is to give a theoretical explanation of this
phenomenon.
Several authors have performed a linear stability analysis
([12,19]),which however did not predict a
wavelength selection mechanism.Instead they found that th
e infinitesimal growth rate is independent of
the wave number
|
!
|
in the limit
|
!
| *
0.Saintillan et al.[18] showed that a small density stratifi
cation
can lead to a finite wavelength selection.However,the predi
cted wavelength does not agree with
experiments,see [15].In those papers the Stokes equation w
as used to model the macroscopic flow.
Here we show that the use of the Navier-Stokes equation leads
to a wavelength selection mechanism
which depends on the dimensionless parameters
Re
,
(
and
'
.
In this section we restrict considerations to the spatially
two-dimensional case (one horizontal
direction + direction of gravity).In this configuration we l
ook at the Smoluchowski equation on
S
1
and
use the notation
x
= (
x,y
)
t
and
u
= (
u,v
)
t
.We linearize the Smoluchowski equation around the state
u
0
=
0
and
f
0
= 1
/
2
*
and consider small perturbations of
f
0
which we denote by
f
and perturbations
of
u
0
which we denoted by
u
= (
u,v
)
t
.The evolution of
f
is described by the equation
$
t
f
+
&
n

(
P
n
!
&
x
un
)
1
2
*
$
(
D
(
n
)
e
2
)
∙ &
x
f
= $
n
f
+
'
&
x

D
(
n
)
&
x
f
(3.1)
10
with
n
=
!
cos
!
sin
!
"
,D
(
n
) = (
I
+
n
%
n
) =
!
cos
2
!
+1 cos
!
sin
!
cos
!
sin
!
sin
2
!
+1
"
.
Now we consider a closure at the level of second moments
&
(
x
,t
):=
&
2
#
0
f
(
x
,!,t
)
d!,c
(
x
,t
):=
1
2
&
2
#
0
cos(2
!
)
f d!,s
(
x
,t
):=
1
2
&
2
#
0
sin(2
!
)
f d!
(3.2)
which allows us to derive evolution equations for
&
,
c
and
s
,see Appendix A for details.Note that
&
now describes the perturbation of the density.This closure
makes sense in the regime where
(
is small.
In this situation,the velocity gradient caused by modulati
ons in density is small and the rotational
di!usion term in the Smoluchowski equation will lead to a fas
t damping of higher order moments.
The linear model has the form
$
t
&
=
3
2
&
y
$
c
y
+
s
x
+
'
!
3
2
(
&
xx
+
&
yy
) +2
s
xy
+
c
xx
$
c
yy
"
(3.3)
$
t
c
=
$
1
8
&
y
+
3
2
c
y
+
1
4
(
u
x
$
v
y
)
$
4
c
+
'
'
3
2
(
c
xx
+
c
yy
) +
1
8
(
&
xx
$
&
yy
)
(
(3.4)
$
t
s
=
1
8
&
x
+
3
2
s
y
+
1
4
(
u
y
+
v
x
)
$
4
s
+
'
!
1
4
&
xy
+
3
2
(
s
xx
+
s
yy
)
"
.
(3.5)
In a moving coordinate frame of the form (
x
!
,t
#
) = (
x
+
3
2
e
2
t,t
) we obtain the system (dropping
the primes)
$
t
&
=
s
x
$
c
y
+
'
!
3
2
(
&
xx
+
&
yy
) +
c
xx
$
c
yy
+2
s
xy
"
(3.6)
$
t
c
=
$
1
8
&
y
+
1
4
(
u
x
$
v
y
)
$
4
c
+
'
'
1
8
(
&
xx
$
&
yy
) +
3
2
(
c
xx
+
c
yy
)
(
(3.7)
$
t
s
=
1
8
&
x
+
1
4
(
u
y
+
v
x
)
$
4
s
+
'
!
1
4
&
xy
+
3
2
(
s
xx
+
s
yy
)
"
.
(3.8)
3.1.The Stokes equation.
We consider (3.6)-(3.8) together with the Stokes equation (
Re
= 0)
$
u
$
p
x
+
('
(2
c
x
+2
s
y
) = 0
$
v
$
p
y
+
('
(2
s
x
$
2
c
y
) =
(&
u
x
+
v
y
= 0
.
(3.9)
Since here
&
is the perturbation of the density,we have on the right hand s
ide the term
(&
instead of
(
(
&
$
1) (compare with (2.28)).
We rewrite (3.6)-(3.8) as a system for
&
,
+
=
s
x
$
c
y
and
,
=
s
y
+
c
x
:
$
t
&
=
+
+
'
!
3
2
$
&
+
+
y
+
,
x
"
$
t
+
=
1
8
$
&
+
1
4
$
v
$
4
+
+
'
!
1
8
$(
&
y
) +
3
2
$
+
"
$
t
,
=
1
4
$
u
$
4
,
+
'
!
1
8
$(
&
x
) +
3
2
$
,
"
.
The terms $
u
and $
v
are replaced using the Stokes equation,i.e.
$
u
=
p
x
$
('
(2
c
x
+2
s
y
) =
p
x
$
('
2
,
$
v
=
p
y
$
('
(2
s
x
$
2
c
y
) +
(&
=
p
y
$
('
2
+
+
(&.
11
We obtain the system
$
t
&
=
+
+
'
!
3
2
$
&
+
+
y
+
,
x
"
$
t
+
=
1
8
$
&
+
1
4
(
p
y
$
('
2
+
+
(&
)
$
4
+
+
'
!
1
8
$(
&
y
) +
3
2
$
+
"
$
t
,
=
1
4
(
p
x
$
('
2
,
)
$
4
,
+
'
!
1
8
$(
&
x
) +
3
2
$
,
"
.
After Fourier transformation we obtain with
!
= (
-
1
,-
2
)
t
$
t
ˆ
&
(
!
,t
) =
ˆ
+
(
!
,t
) +
'
!
$
3
2
|
!
|
2
ˆ
&
(
!
,t
) +
i-
2
ˆ
+
(
!
,t
) +
i-
1
ˆ
,
(
!
,t
)
"
$
t
ˆ
+
(
!
,t
) =
$
1
8
|
!
|
2
ˆ
&
(
!
,t
) +
1
4
i-
2
ˆ
p
(
!
,t
)
$
1
2
('
ˆ
+
(
!
,t
) +
1
4
(
ˆ
&
(
!
,t
)
$
4
ˆ
+
(
!
,t
) +
'
!
$
1
8
|
!
|
2
i-
2
ˆ
&
(
!
,t
)
$
3
2
|
!
|
2
ˆ
+
(
!
,t
)
"
$
t
ˆ
,
(
!
,t
) =
1
4
i-
1
ˆ
p
$
1
2
('
ˆ
,
(
!
,t
)
$
4
ˆ
,
(
!
,t
)
+
'
!
$
1
8
|
!
|
2
i-
1
ˆ
&
(
!
,t
)
$
3
2
|
!
|
2
ˆ
,
(
!
,t
)
"
.
(3.10)
We want to calculate ˆ
p
from the Stokes equations.For this we di!erentiate the first
equation of (3.9)
with respect to
x
,the second equation with respect to
y
and add the two equations.We obtain
$(
u
x
+
v
y
)
0
12
3
=0
$
$
p
+2
('
(
,
x
+
+
y
) =
(&
y
.
Fourier transformation gives us
ˆ
p
=
$
('
2
|
!
|
2
'
i-
1
ˆ
,
+
i-
2
ˆ
+
(
+
(
i-
2
|
!
|
2
ˆ
&.
(3.11)
In (3.10) we need expressions for
i-
1
ˆ
p/
4 and
i-
2
ˆ
p/
4,which have the form
1
4
i-
1
ˆ
p
=
('
2
|
!
|
2
'
-
2
1
ˆ
,
+
-
1
-
2
ˆ
+
(
$
(
4
-
1
-
2
|
!
|
2
ˆ
&
1
4
i-
2
ˆ
p
=
('
2
|
!
|
2
'
-
1
-
2
ˆ
,
+
-
2
2
ˆ
+
(
$
(
4
-
2
2
|
!
|
2
ˆ
&
(3.12)
Using the expressions (3.12) in (3.10) we obtain the system
$
t
ˆ
+
=
ˆ
+
!
$
4
$
3
2
'
|
!
|
2
$
1
2
('
-
2
1
|
!
|
2
"
+
ˆ
,
!
1
2
('
-
1
-
2
|
!
|
2
"
+ ˆ
&
!
$
1
8
|
!
|
2
+
1
4
(
-
2
1
|
!
|
2
$
1
8
'i-
2
|
!
|
2
"
$
t
ˆ
,
=
ˆ
+
!
1
2
('
-
1
-
2
|
!
|
2
"
+
ˆ
,
!
$
4
$
3
2
'
|
!
|
2
$
1
2
('
-
2
2
|
!
|
2
"
+ ˆ
&
!
$
1
4
(
-
1
-
2
|
!
|
2
$
1
8
'i-
1
|
!
|
2
"
$
t
ˆ
&
=
ˆ
+
(1 +
'i-
2
) +
ˆ
,
(
'i-
1
) + ˆ
&
!
$
3
2
'
|
!
|
2
"
(3.13)
We have rewritten the problem into a linear ODE system of the f
orm
$
t
U
(
!
,t
) =
A
(
!
,(,'
)
U
(
!
,t
)
,
12
with
U
(
!
,t
) =
'
ˆ
+,
ˆ
,,
ˆ
&
(
T
and
A
(
!
,(,'
) equal to
5
6
7
$
4
$
3
2
'
|
!
|
2
$
1
2
('
$
2
1
|
!
|
2
1
2
('
$
1
$
2
|
!
|
2
$
1
8
|
!
|
2
+
1
4
(
$
2
1
|
!
|
2
$
1
8
'i-
2
|
!
|
2
1
2
('
$
1
$
2
|
!
|
2
$
4
$
3
2
'
|
!
|
2
$
1
2
('
$
2
2
|
!
|
2
$
1
4
(
$
1
$
2
|
!
|
2
$
1
8
'i-
1
|
!
|
2
1 +
'i-
2
'i-
1
$
3
2
'
|
!
|
2
8
9
:
With
.
(
A
) we denote the spectral abscissa,which is defined by
.
(
A
):= max
1
$
j
$
3
Re(
)
j
)
,
where
)
j
,
j
= 1
,
2
,
3 are the eigenvalues of
A
.Now
U
(
!
,t
) remains bounded provided
.
(
A
)
+
0 and
U
(
!
,t
)
*
0 as
t
*,
if
.
(
A
)
<
0.
In Figure 3.1 we show plots of the positive part of the spectra
l abscissa for di!erent values of
(
and
'
.Since we can not give analytical expressions of the eigenva
lues of the matrix
A
,we calculated (using
postive part of spectral abscissa for
δ
=0.01,
γ
=0
ξ
1
ξ
2


0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
1
2
3
4
5
6
x 10
−4
postive part of spectral abscissa for
δ
=0.01,
γ
=0.01
ξ
1
ξ
2


0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
1
2
3
4
5
6
x 10
−4
postive part of spectral abscissa for
δ
=0.01,
γ
=0.1
ξ
1
ξ
2


0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
1
2
3
4
5
6
x 10
−4
postive part of spectral abscissa for
δ
=0.1,
γ
=0
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
1
2
3
4
5
6
x 10
−3
postive part of spectral abscissa for
δ
=0.1,
γ
=0.01
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
1
2
3
4
5
6
x 10
−3
postive part of spectral abscissa for
δ
=0.1,
γ
=0.1
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
1
2
3
4
5
6
x 10
−3
postive part of spectral abscissa for
δ
=1,
γ
=0
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.01
0.02
0.03
0.04
0.05
0.06
postive part of spectral abscissa for
δ
=1,
γ
=0.01
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.01
0.02
0.03
0.04
0.05
0.06
postive part of spectral abscissa for
δ
=1,
γ
=0.1
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.01
0.02
0.03
0.04
0.05
0.06
Fig.3.1
.
Positive part of spectral abscissa for the matrix A vs.
#
1
and
#
2
for di!erent values of
$
and
%
.
MATLAB) the positive part of the spectral abscissa on a fine me
sh in the
-
1
-
-
2
-plane.We summarize
our observations in Remark 3.1.
Remark 3.1.
On the level of linear instability,there is no wavelength se
lection mechanism for the
coupled micro-macro model with the macroscopic Stokes equa
tion.
An increase of the parameter
(
has a destabilizing e!ect,while an increase of the paramete
r
'
has
a stabilizing e!ect.
13
3.2.The Navier-Stokes equation.
We now study the linear stability for the system (3.3)-(3.5)
together with the Navier-Stokes equation
Re$
t
u
=
u
xx
+
u
yy
$
p
x
+
('
2(
c
x
+
s
y
) (3.14)
Re$
t
v
=
v
xx
+
v
yy
$
p
y
+
('
2(
s
x
$
c
y
)
$
(&.
(3.15)
Since we linearize around the velocity
u
0
= 0,the e!ect of the convective term is neglected.Note that
we here use the linear moment closure system in the non-movin
g frame (instead of transforming the
Navier-Stokes equation to the moving coordinate system).
A Fourier transformation of (3.3)-(3.5),(3.14),(3.15) le
ads to the system
$
t
ˆ
&
=
3
2
i-
2
ˆ
&
+
i-
1
ˆ
s
$
i-
2
ˆ
c
+
'
!
$
3
2
|
!
|
2
ˆ
&
$
-
2
1
ˆ
c
+
-
2
2
ˆ
c
$
2
-
1
-
2
ˆ
s
"
(3.16)
$
t
ˆ
c
=
$
1
8
i-
2
ˆ
&
+
3
2
i-
2
ˆ
c
+
1
4
(
i-
1
ˆ
u
$
i-
2
ˆ
v
)
$

c
+
'
!
1
8
*
$
-
2
1
ˆ
&
+
-
2
2
ˆ
&
+
$
3
2
|
!
|
2
ˆ
c
"
(3.17)
$
t
ˆ
s
=
1
8
i-
1
ˆ
&
+
3
2
i-
2
ˆ
s
+
1
4
(
i-
2
ˆ
u
+
i-
1
ˆ
v
)
$

s
+
'
*
$
1
4
-
1
-
2
ˆ
&
$
3
2
|
!
|
2
ˆ
s
+
(3.18)
$
t
ˆ
u
=
1
Re
*
$|
!
|
2
ˆ
u
$
i-
1
ˆ
p
+2
('
(
i-
1
ˆ
c
+
i-
2
ˆ
s
)
+
(3.19)
$
t
ˆ
v
=
1
Re
*
$|
!
|
2
ˆ
v
$
i-
2
ˆ
p
+2
('
(
i-
1
ˆ
s
$
i-
2
ˆ
c
)
$
(
ˆ
&
+
(3.20)
We eliminate the pressure term ˆ
p
using the expression
ˆ
p
=
2
('
|
!
|
2
*
2
-
1
-
2
ˆ
s
+
-
2
1
ˆ
c
$
-
2
2
ˆ
c
+
+
(
|
!
|
2
i-
2
ˆ
&.
(3.21)
Equations (3.19) and (3.20) can thus be written in the form
$
t
ˆ
u
=
1
Re
!
(
-
1
-
2
|
!
|
2
ˆ
&
$|
!
|
2
ˆ
u
+2
('-
2
i
!
1
$
2
-
2
1
|
!
|
2
"
ˆ
s
+
4
('i-
1
-
2
2
|
!
|
2
ˆ
c
"
(3.22)
$
t
ˆ
v
=
1
Re
!
$
(
-
2
1
|
!
|
2
ˆ
&
$|
!
|
2
ˆ
v
+2
('-
1
i
!
1
$
2
-
2
2
|
!
|
2
"
ˆ
s
$
4
('-
2
i-
2
1
|
!
|
2
ˆ
c
"
(3.23)
We have rewritten the problem into a linear ODE system of the f
orm
$
t
U
(
!
,t
) =
A
(
!
,(,',Re
)
U
(
!
,t
)
,
with
U
(
!
,t
) = (ˆ
&,
ˆ
c,
ˆ
s,
ˆ
u,
ˆ
v
)
t
and
A
equal to
5
6
6
6
6
6
7
!
3
2
!
|
!
|
2
+
3
2
i"
2
!
i"
2
+
!
(
"
2
2
!
"
2
1
)
i"
1
!
2
!"
1
"
2
0 0
1
8
!
!
i"
2
+
!
(
"
2
2
!
"
2
1
)
"
!
4
!
3
2
!
|
!
|
2
+
3
2
i"
2
0
1
4
i"
1
!
1
4
i"
2
1
8
i"
1
!
1
4
!"
1
"
2
0
!
4
!
3
2
!
|
!
|
2
+
3
2
i"
2
1
4
i"
2
1
4
i"
1
!"
1
"
2
Re
|
!
|
2
4
!#i"
1
"
2
2
Re
|
!
|
2
2
!#i"
2
Re
#
1
!
2
"
2
1
|
!
|
2
$
!
|
!
|
2
Re
0
!
!"
2
1
Re
|
!
|
2
!
4
!#i"
2
"
2
1
Re
|
!
|
2
2
!#i"
1
Re
#
1
!
2
"
2
2
|
!
|
2
$
0
!
|
!
|
2
Re
8
9
9
9
9
9
:
In Figures 3.2,3.3 we plot the positive part of the spectral a
bscissa for di!erent choices of the
parameter values
Re
,
(
and
'
.
Remark 3.2.
For the coupled micro-macro model with the Navier-Stokes eq
uation,the linear
stability analysis predicts a wavelength selection which d
epends on
Re
,
(
and
'
.An increase of
Re
decreases the wavelength of the most unstable configuration
.An increase of
(
has a destabilizing e!ect
and decreases the wavelength of the most unstable wave.An in
crease of
'
has a stabilizing e!ect and
increases the wavelength of the most unstable wave.
14
(a)
positive part of the spectral abscissa (Re=0.001,
δ
=0.1,
γ
=0)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
1
2
3
4
5
6
x 10
−3
(b)
positive part of the spectral abscissa (Re=0.001,
δ
=0.1,
γ
=0.1)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
1
2
3
4
5
6
x 10
−3
(c)
positive part of the spectral abscissa (Re=0.1,
δ
=0.1,
γ
=0)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
5.5
x 10
−3
(d)
positive part of the spectral abscissa (Re=0.1,
δ
=0.1,
γ
=0.1)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
−3
(e)
positive part of the spectral abscissa (Re=1,
δ
=0.1,
γ
=0)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
−3
(f)
positive part of the spectral abscissa (Re=1,
δ
=0.1,
γ
=0.1)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
1
1.5
2
2.5
3
x 10
−3
Fig.3.2
.
Positive part of the spectral abscissa for
$
= 0
.
1
and di!erent values of the parameter values
Re
and
%
.
3.3.Asymptotic expansion of the relevant eigenvalue.
It is expected that the most unstable
wavelength lies in the nearly horizontal direction.This mo
tivates to consider the special case
-
2
= 0
and (for simplicity)
'
= 0.We give an asymptotic expansion of the unstable eigenval
ue of the linearized
systemin
Re
.This asymptotic expansion captures the singular e!ect of t
he dependence on the Reynolds
number and explains why the inclusion of inertial e!ects lea
ds to a wavelength selection.
Proposition 3.1.
Consider horizontal waves in the regime
'
= 0
.For
Re
= 0
,the coupled system
is linearly stable provided that the eigenvalue
)
0
=
1
4
!
$
8 +
;
64 +4
(
$
2
-
2
1
"
15
(a)
positive part of the spectral abscissa (Re=0.001,
δ
=1,
γ
=0)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.01
0.02
0.03
0.04
0.05
0.06
(b)
positive part of the spectral abscissa (Re=0.001,
δ
=1,
γ
=0.1)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
0.055
(c)
positive part of the spectral abscissa (Re=0.1,
δ
=1,
γ
=0)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
0.055
(d)
positive part of the spectral abscissa (Re=0.1,
δ
=1,
γ
=0.1)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
(e)
positive part of the spectral abscissa (Re=1,
δ
=1,
γ
=0)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
(f)
positive part of the spectral abscissa (Re=1,
δ
=1,
γ
=0.1)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.005
0.01
0.015
0.02
0.025
0.03
Fig.3.3
.
Positive part of the spectral abscissa for
$
= 1
and di!erent values of the parameter values
Re
and
%
.
is smaller than zero,otherwise it is unstable.
For
Re >
0
(and
Re
!
1
) the change of this eigenvalue can be described by an asympto
tic expansion
of the form
)
=
)
0
+(
Re
)
)
1
+
...
with
)
1
=
$
(
2
$
2
1
(2
%
!
$
2
1
)
'
8 +
<
64 +4
(
$
2
-
2
1
(
2
+4
-
2
1
.
16
Proof:
Our starting point is the linear system which is obtained via
Fourier transformation of the
linear moment closure.For
-
2
= 0,
'
= 0,this system has the form (dropping all the hats)
$
t
&
=
i-
1
s
$
t
c
=
1
4
i-
1
u
$
4
c
$
t
s
=
1
8
i-
1
&
+
1
4
i-
1
v
$
4
s
(3.24)
Using the notation
x
= (
&,c,s
)
T
,
y
= (
u,v
)
T
,we rewrite (3.24) into the form
˙
x
=
A
x
+
B
y
,
(3.25)
with
A
=
5
7
0 0
i-
1
0
$
4 0
1
8
i-
1
0
$
4
8
:
,B
=
5
7
0 0
1
4
i-
1
0
0
1
4
i-
1
8
:
.
From the Stokes equation (
Re
= 0) we obtain
$
-
2
1
u
= 0
$
-
2
1
v
$
(&
= 0
,
which we rewrite into the form
0 =
C
y
+
D
x
,
(3.26)
with
C
=
!
$
-
2
1
0
0
$
-
2
1
"
,D
=
!
0 0 0
$
(
0 0
"
.
Using (3.26),we obtain
y
=
$
C
!
1
D
x
.Introducing this in (3.25) gives
˙
x
=
*
A
$
BC
!
1
D
+
x
.
(3.27)
The eigenvalues of (
A
$
BC
!
1
D
) are
=
$
4
,
1
4
!
$
8
$
;
64 +4
(
$
2
-
2
1
"
,
1
4
!
$
8 +
;
64 +4
(
$
2
-
2
1
">
.
(3.28)
The last eigenvalue (which we denote by
)
0
) is larger than zero provided that
( >
0 and
-
2
1
is small
enough.Thus,for horizontal waves and
'
= 0,the moment closure system coupled with the Stokes
equation is most unstable for waves with wave number
-
1
*
0,compare with Figure 3.1 (left column).
The eigenvector corresponding to the eigenvalue
)
0
has the form
x
0
=
!
2
i-
1
2
(
$
-
2
1
!
8 +
;
64 +4
(
$
2
-
2
1
"
,
0
,
1
"
T
.
(3.29)
Now we consider the case
Re >
0,for which the linearized coupled system can be expressed i
n the
form
!
˙
x
˙
y
"
=
!
A B
1
Re
D
1
Re
C
"!
x
y
"
.
(3.30)
17
Our goal is to give an asymptotic expansion for eigenvalues o
f the matrix arising on the right hand side
of (3.30),which is valid for small values of
Re
.In particular we wish to understand how the eigenvalue
)
0
changes with
Re
.We make the ansatz
A
x
&
+
B
y
&
=
)
&
x
&
C
y
&
+
D
x
&
=
/)
&
y
&
,
(3.31)
with
)
&
=
)
0
+
/)
1
+
...
x
&
=
x
0
+
/
x
1
+
...
y
&
=
y
0
+
/
y
1
+
....
We obtain
O
(1):(as considered above)
A
x
0
+
B
y
0
=
)
0
x
0
C
y
0
+
D
x
0
= 0
(3.32)
O
(
/
):
A
x
1
+
B
y
1
=
)
0
x
1
+
)
1
x
0
C
y
1
+
D
x
1
=
)
0
y
0
,
(3.33)
which we rewrite into the form
!
A
$
)
0
I B
D C
"!
x
1
y
1
"
=
!
)
1
x
0
)
0
y
0
"
.
(3.34)
The left null space of the matrix on the left hand side of (3.34
) is given by
u
T
=
!
$
8
$
;
64 +4
(
$
2
-
2
1
,
0
,
$
4
i-
1
,
0
,
1
"
.
(3.35)
We multiply both sides of Equation (3.34) from the left with
u
T
and obtain
0 =
u
T

!
)
1
x
0
)
0
y
0
"
.
(3.36)
Here
y
0
=
$
C
!
1
D
x
0
=
$
?
0
2
%i
$
1
(8+
-
64+4
%
!
2
$
2
1
)
2
%
!
$
2
1
@
.
Using the expressions for
)
0
,
x
0
,
y
0
and
u
in Equation (3.36),we can now calculate
)
1
,which has the
form
)
1
=
$
(
2
$
2
1
(2
%
!
$
2
1
)
'
8 +
<
64 +4
(
$
2
-
2
1
(
2
+4
-
2
1
.
(3.37)
In Figure 3.4 we show a plot of
)
calculated by the asymptotic expansion (solid line) and com
pare
it with the positive part of the spectral abscissa (calculat
ed at discrete points) for purely horizontal
waves.
18
0
0.5
1
1.5
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
positive part of the spectral abscissa (Re=10
−3
,
δ
= 1,
γ
=0)
Fig.3.4
.
Solid line:eigenvalue calculated from asymptotic expansion,’+
’-symbol:positive part of spectral abscissa
for
Re
= 0
.
001
and
$
= 1
.
3.4.Numerical simulations of the nonlinear model.
Here we show numerical simulations for
the system
$
t
f
+
$
!
*
v
x
cos
2
!f
+
$
$
x
(sin
!
cos
!f
) =
$
!!
f
$
t
v
=
v
xx
+
(
!
1
$
&
S
1
fd!
"
,
(3.38)
i.e.we solve the coupled model in a simple setting with
u
=
5
7
0
v
(
x,t
)
0
8
:
,f
=
f
(
x,!,t
)
,
We use
(
= 1 (and
'
= 0,
Re
= 1).
Initial values are specified on an interval of length
L
= 500,periodicity is assumed.The initial
values for
f
have the form
f
(
x,
0
,!
) =
1
2
*
+
p
(
x
)
x
"
[0
,L
)
,
.
!
(3.39)
where
p
(
x
) is a small perturbation with size of the order 5
(
10
!
7
,see Figure 3.5.The velocity
v
is
initially set to zero.
We discretize the macroscopic domain by 8000 grid cells.At e
ach grid point of the macroscopic
domain we discretize the Smoluchowski equation on
S
1
using 100 grid cells.
In Figure 3.6 we show numerical simulations for the nonlinea
r initial value problem.In 3.6 (a) we
show a plot of
f
as a function of
x
(horizontal axis) and
!
(vertical axis),(b) shows a plot of
&
.The
results of a discrete Fourier transformation of
&
$
1 as a function of
|
!
|
=
|
-
1
|
are shown in Figure 3.6
(c),(d).The numerical result obtained for the nonlinear mo
del fits well to the prediction of the linear
theory,compare with Figure 3.3 (e).Instabilities are only
observed within the domain predicted by the
linear theory.Furthermore,the most unstable waves are obs
erved for wave number
|
-
1
|/
0
.
5.
In Figure 3.7 we show plots of the density
&
at later times.We observe regions with high concen-
tration which are separated by regions with very small conce
ntration of particles.This is the nonlinear
regime,the solution structure can no longer be predicted by
a linearized model.In Section 5 we show
that for a particular scaling of our model the evolution of
&
can be described by a Keller-Segel model.
4.Sedimentation of non-Brownian particles.
In experimental and numerical studies that
can be found in the literature (e.g.[14,17]),the authors di
scuss the sedimentation for suspensions of
19
0
50
100
150
200
250
300
350
400
450
500
1
1
1
1
1
1
1
1
1
density at initial time
0
10
20
30
40
50
60
0
1
2
3
4
5
6
7
8
x 10
−8
|
ξ
|
fft(
ρ
−1) vs. wave number at initial time
Fig.3.5
.
Initial values of density
&
and
fft
(
&
!
1)
.
(a)
(b)
0
50
100
150
200
250
300
350
400
450
500
0.998
0.9985
0.999
0.9995
1
1.0005
1.001
1.0015
1.002
1.0025
ρ
at time t=250
(c)
0
10
20
30
40
50
60
0
0.5
1
1.5
x 10
−4
|
ξ
|
fft(
ρ
−1) vs. wave number at time t=250
(d)
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0
1
2
x 10
−4
|
ξ
|
zoom of fft(
ρ
−1) vs. wave number at time t=250
Fig.3.6
.
Numerical results for the initial value problem (3.38),(3.3
9) at time
t
= 250
.(a)
f
(
',x
)
,(b)
&
(
x
)
,(c)
fft
(
&
!
1)
vs.
|
#
1
|
,(d) closeup view of results from Fourier transformation.
non-Brownian particles.For our kinetic model,this case ca
n be considered by setting the absolute
20
(a)
0
50
100
150
200
250
300
350
400
450
500
0.85
0.9
0.95
1
1.05
1.1
1.15
1.2
1.25
ρ
at time t=400
(b)
0
50
100
150
200
250
300
350
400
450
500
0.5
1
1.5
2
2.5
3
3.5
ρ
at time t=450
(c)
0
50
100
150
200
250
300
350
400
450
500
0
5
10
15
20
25
30
35
40
ρ
at time t=500
(d)
0
5
10
15
20
25
30
35
40
45
50
0
5
10
15
20
25
30
35
40
zoom of
ρ
at time t=500
Fig.3.7
.
Numerical results for the initial value problem (3.38) at lat
er times.Figure (d) showes a closeup view of
the solution from (c)
temperature equal to zero.We then obtain the system
$
t
f
=
$&
x

(
u
f
)
$&
n

(
P
n
!
&
x
un
f
) +
m
0
g
#
"
&
x

(
D
(
n
)
e
2
f
)
&
f
(
$
t
u
+(
u
∙ &
x
)
u
) =
µ
$
u
$&
x
p
+
!
N
V
$
&
S
2
fd
n
"
m
0
g
e
2
0 =
&
x

u
.
(4.1)
In order to non-dimensionalize the system we consider a scal
ing of the form
t
=
T
ˆ
t,
x
=
X
ˆ
x
,
u
=
X
T
ˆ
u
,f
=
N
V
ˆ
f,p
=
µ
T
ˆ
p.
where
V
=
O
(
L
3
) is the volume occupied by the suspension,and
L
is the size of the macroscopic
domain.Time and space scales are selected according to
T
=
m
0
#
"
and
X
=
m
2
0
g
#
2
"
,
(4.2)
21
which corresponds to the velocity scale (2.20).We obtain th
e non-dimensional equations (dropping all
the hats)
$
t
f
=
$&
x

(
u
f
)
$&
n

(
P
n
!
&
x
un
f
) +
&
x

(
D
(
n
)
e
2
f
)
Re
(
$
t
u
+(
u
∙ &
x
)
u
) = $
x
u
$&
x
p
+
(
!
1
$
&
S
2
fd
n
"
e
2
0 =
&
x

u
,
(4.3)
with non-dimensional parameters
Re
=
&
f
X
2
µT
and
(
=
N
V
m
0
gTX
µ
.
In analogy to the previous section,we now present a linear st
ability analysis for this slightly sim-
pler model.By neglecting di!usion there is no justification
to restrict considerations to second order
moments,since there is no fast damping of higher moments.Ho
wever,our numerical simulations of the
nonlinear model will show good agreement with the wavelengt
h selection mechanism predicted by the
linear model.
We consider the spatially two-dimensional case and lineari
ze around the state
u
0
= 0 and
f
0
=
1
2
#
,
f
and
u
now denote small perturbations of
f
0
and
u
0
.The linearized equations have the form
$
t
f
=
$&
n

(
P
n
!
&
x
un
)
1
2
*
+(
D
(
n
)
e
2
)
∙ &
x
f
Re$
t
u
= $
x
u
$&
x
p
$
(
!
&
S
2
fd
n
"
e
2
.
(4.4)
We consider a moment closure on the level of
&
,
c
,
s
as defined in (3.2) and obtain the evolution
equations
$
t
&
=
s
x
+
3
2
&
y
$
c
y
$
t
c
=
3
2
c
y
$
1
8
&
y
+
1
4
(
u
x
$
v
y
)
$
t
s
=
1
8
&
x
+
3
2
s
y
+
1
4
(
u
y
+
v
x
)
.
(4.5)
By changing to a moving coordinate system,we get
$
t
&
=
s
x
$
c
y
$
t
c
=
$
1
8
&
y
+
1
4
(
u
x
$
v
y
)
$
t
s
=
1
8
&
x
+
1
4
(
u
y
+
v
x
)
.
(4.6)
We first consider the system (4.6) together with the Stokes eq
uation
$
u
$
p
x
= 0
$
v
$
p
y
=
(&
u
x
+
v
y
= 0
.
(4.7)
Using
+
=
s
x
$
c
y
,we can rewrite (4.6) into the form
$
t
&
=
+
$
t
+
=
1
8
$
&
+
1
4
$
v,
and by using (4.7) we obtain
$
t
&
=
+
$
t
+
=
1
8
$
&
+
1
4
(
(&
+
p
y
)
.
22
Using Fourier transformation and elimination of the pressu
re term,we derive the linear ode system
$
t
!
ˆ
&
ˆ
+
"
=
?
0 1
$
1
8
|
!
|
2
+
1
4
(
$
2
1
|
!
|
2
0
@
!
ˆ
&
ˆ
+
"
.
(4.8)
In Figure 4.1,we plot contour lines in the
-
1
$
-
2
-plane for the positive part of the spectral abscissa
for the matrix on the right hand side of equation (4.8).At the
level of linear stability analysis there is
positive part of the spectral abscissa (
δ
=0.01)
ξ
1
ξ
2


0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
positive part of the spectral abscissa (
δ
=0.1)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.02
0.04
0.06
0.08
0.1
0.12
0.14
positive part of the spectral abscissa (
δ
=1)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Fig.4.1
.
Positive part of spectral abscissa for the matrix arising in t
he linear stability analysis of the system
modeling the sedimentation of non-Brownian particles in a Stokes fl
ow (with
$
= 0
.
01
,
0
.
1
,
1
).
no wavelength selection mechanism for the micro-macro mode
l (4.3) with
Re
= 0.An increase of the
parameter
(
has a destabilizing e!ect.
We now consider system (4.5) together with the linearized Na
vier-Stokes equation
Re$
t
u
= $
u
$
p
x
Re$
t
v
= $
v
$
p
y
$
(&
0 =
u
x
+
v
y
.
(4.9)
Fourier transformation and elimination of the pressure lea
ds to the linear ode system
$
t
5
6
6
6
6
7
ˆ
&
ˆ
c
ˆ
s
ˆ
u
ˆ
v
8
9
9
9
9
:
=
5
6
6
6
6
6
7
3
2
i-
2
$
i-
2
i-
1
0 0
$
1
8
i-
2
3
2
i-
2
0
1
4
i-
1
$
1
4
i-
2
1
8
i-
1
0
3
2
i-
2
1
4
i-
2
1
4
i-
1
%$
1
$
2
Re
|
!
|
2
0 0
$
1
Re
|
!
|
2
0
$
%$
2
1
Re
|
!
|
2
0 0 0
$
1
Re
|
!
|
2
8
9
9
9
9
9
:
5
6
6
6
6
7
ˆ
&
ˆ
c
ˆ
s
ˆ
u
ˆ
v
8
9
9
9
9
:
(4.10)
In Figure 4.2,we plot the positive part of the spectral absci
ssa for the matrix on the right hand side of
equation (4.10) for di!erent values of
Re
and
(
.We observe a wavelength selection mechanism which
depends on
Re
and
(
.
Proposition 4.1.
Consider horizontal waves,i.e.
-
2
= 0
.For
Re
= 0
,the coupled system
describing the sedimentation of non-Brownian particles is
linearly stable provided that the eigenvalue
)
0
=
A
(
4
$
-
2
1
8
is smaller than zero,otherwise it is unstable.
For
Re >
0
,(
Re
!
1
) the change of this eigenvalue can be expressed by an asympto
tic expansion
of the form
)
=
A
(
4
$
-
2
1
8
$
(
Re
)
(
8
-
2
1
+
....
(4.11)
23
positive part of the spectral abscissa (Re = 0.001,
δ
=0.01)
ξ
1
ξ
2


0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
positive part of the spectral abscissa (Re = 0.001,
δ
=0.1)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.02
0.04
0.06
0.08
0.1
0.12
0.14
positive part of the spectral abscissa (Re = 0.001,
δ
=1)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
positive part of the spectral abscissa (Re = 0.1,
δ
=0.01)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0
0.05
0.1
0.15
0.2
0.25
0.3
0.005
0.01
0.015
0.02
0.025
0.03
positive part of the spectral abscissa (Re = 0.1,
δ
=0.1)
ξ
1
ξ
2


0
0.1
0.2
0.3
0.4
0.5
0
0.1
0.2
0.3
0.4
0.5
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0.11
positive part of the spectral abscissa (Re = 0.1,
δ
=1)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
positive part of the spectral abscissa (Re = 1,
δ
=0.01)
ξ
1
ξ
2


0
0.05
0.1
0.15
0.2
0.25
0.3
0
0.05
0.1
0.15
0.2
0.25
0.3
1
2
3
4
5
6
7
8
9
10
11
x 10
−3
positive part of the spectral abscissa (Re = 1,
δ
=0.1)
ξ
1
ξ
2


0
0.1
0.2
0.3
0.4
0.5
0.6
0
0.1
0.2
0.3
0.4
0.5
0.6
0.01
0.02
0.03
0.04
0.05
0.06
positive part of the spectral abscissa (Re = 1,
δ
=1)
ξ
1
ξ
2


0
0.2
0.4
0.6
0.8
1
1.2
1.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.05
0.1
0.15
0.2
0.25
0.3
Fig.4.2
.
Positive part of spectral abscissa for the matrix arising in t
he linear stability analysis of the systemmodeling
the sedimentation of non-Brownian particles in a flow with
Re >
0
.Di!erent values of
Re
and
$
are considered.
The proof can be carried out in analogy to Proposition 3.1.Us
ing the same notation,the eigenvalues
of the matrix
A
$
BC
!
1
D
are now
{
0
,
$
A
(
4
$
-
2
1
8
,
A
(
4
$
-
2
1
8
}
.
The eigenvector corresponding to the last eigenvalue (whic
h we denote with
)
0
) is
x
=
?
2
-
2
-
1
i
<
2
(
$
-
2
1
,
0
,
1
@
T
.
The left null space of the matrix
!
A
$
)
0
I B
D C
"
is the vector
!
$
;
4
(
$
2
-
2
1
,
0
,
$
4
-
1
i,
0
,
1
"
.
Using all these terms,one can carry out the calculation lead
ing to (4.11).
24
For horizontal waves we now derive an approximative express
ion for the most unstable wavenumber
as a function of the parameter values
Re
and
(
.We rewrite the leading terms of (4.11) in the form
)
=
-
(
2
A
1
$
-
2
1
2
(
$
Re
(
8
-
2
1
=
-
(
2
-
1
$
x
$
Re
16
1
x
,
with
x
=
$
2
1
2
%
.For
x
!
1,we can use the approximation
)
/
-
(
2
'
1
$
x
2
(
$
Re
16
1
x
=
˜
).
˜
)
attains a maximum for
x
=
1
2
Re
1
/
2
(
!
1
/
4
,which corresponds to the most unstable wavenumber
-
1
=
Re
1
/
4
(
3
/
8
.
(4.12)
In Figure 4.3 we show plots of the spectral abscissa for horiz
ontal waves (solid line) and the approxi-
mation of the largest eigenvalue as calculated from (4.11) (
plus symbols).For small values of
Re
(here
we used
Re
= 10
!
4
and
Re
= 10
!
3
) the expansion (4.11) is a good approximation of the spectra
l
abscissa.The dashed line indicates the approximation of th
e most unstable wavenumber as calculated
from (4.12).
(a)
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
Spectral abscissa (Re=10
−4
,
δ
= 0.1)
(b)
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
Spectral abscissa (Re=10
−3
,
δ
= 0.1)
Fig.4.3
.
Positive part of spectral abscissa for horizontal wave over w
avenumber (solid line),
!
(
#
1
)
as calculated
from (4.11) (plus symbol),most unstable wavenumber as approx
imated by (4.12) (dashed line);(a)
Re
= 10
"
4
,
$
= 0
.
1
,
(b)
Re
"
3
,
$
= 0
.
1
.
Now we compare the predicted most unstable wavelength with e
xperimental studies from [10].
The experiments where performed using glas rods with partic
le density
&
p
= 2
.
25
g
cm
3
,particle length
l
p
= 0
.
108
cm
and diameter
d
p
= 0
.
0102
cm
.The particles are suspended in a fluid with fluid viscosity
µ
= 5
.
5
g
cms
and density
&
f
= 1
.
07
g
cm
3
.The particle volume fraction in the suspension is 0
.
48%.From
those data we calculate the volume of a single particle
V
p
=
*
'
d
p
2
(
2
l
p
= 8
.
825

10
!
6
cm
3
,
the mass of
one particle
m
0
=
&
p
V
p
= 1
.
98562

10
!
5
g,
the number density
N
V
= 0
.
0048
/V
p
= 543
.
9103
cm
!
3
,
and the
translational friction coe"cient
#
"
=
4
#µl
p
ln(
l
p
/d
p
)
= 3
.
1632
g
s
(using the formula from [5,Section 8.3]).Our
characteristic length and time scales are
X
= 3
.
8655

10
!
8
cm
= 3
.
5792

10
!
7
l
p
and
T
= 6
.
2773

10
!
6
s
.
The nondimensional parameters can now be calculated and we o
btain the values
Re
= 4
.
6308
(
10
!
11
and
(
= 4
.
6729

10
!
13
.Via (4.12) the predicted most unstable wavenumber is
-
1
= 6
.
2017

10
!
8
which
corresponds to a predicted selected wavelength of the size
X
$
1
= 5
.
7713
l
p
.This is in good qualitative
agreement with experimental studies which show the formati
on of clusters with a typical width of the
size of a few particle lengths,see [10,11,14].An increase o
f the number of particles (within the
25
dilute regime) leads to an increase of the parameter
(
and via (4.12) to a decrease of the predicted
most unstable wavelength.In experimental studies such a de
pendence of the cluster width from the
concentration of particles was also observed.Furthermore
,the length of the particles as well as the
aspect ration a!ect the most unstable wavelength.
Finally,we show numerical simulations for the nonlinear mo
del
$
t
f
+
$
!
*
v
x
cos
2
!f
+
$
$
x
(sin
!
cos
!f
) = 0
$
t
v
$
v
xx
=
!
1
$
&
S
1
fd!
"
,
with initial values as specified in (3.39).We choose
L
= 500 and discretize the macroscopic domain
[0
,L
) with 8000 grid cells.The Smoluchowski equation is discret
ized using 100 grid cells.In Figure 4.4
we plot the density
&
and the result of the Fourier transformation of the data vs.w
ave number at time
t
= 35.The observed wavelength selection (with a wave number o
f about 0.7) is in good agreement
with the prediction of the linear stability,compare with Fi
gure 4.2 (
(
= 1,
Re
= 1).
(a)
0
50
100
150
200
250
300
350
400
450
500
0.9975
0.998
0.9985
0.999
0.9995
1
1.0005
1.001
1.0015
1.002
1.0025
ρ
at time t=35
(b)
0
0.5
1
1.5
0
1
2
x 10
−4
|
ξ
|
zoom of fft(
ρ
−1) vs. wave number at time t=35
Fig.4.4
.
Simulations of the one-dimensional nonlinear model for
$
= 1
and
Re
= 1
.
5.The hyperbolic and di!usive scalings.Relation to the Kel
ler-Segel model.
In this
section we derive the hyperbolic and di!usive limits for the
system (2.27) with
Re
= 0.The limiting
behavior in the hyperbolic scaling will be described by a Bou
ssinesq type system.For certain flows the
hyperbolic scaling produces a trivial behavior,and it is th
en natural to consider the di!usive scaling.
Such a situation occurs for two-dimensional rectilinear flo
ws of suspensions.We will show that the
collective behavior in the di!usive limit is described by th
e Keller-Segel model
It is expedient to view the scaling limits from the perspecti
ve of the question of describing the
aggregate behavior of a suspension starting from rest.The f
unction
&
(
x
,t
) =
&
S
d
"
1
f
(
x
,t,
n
)
d
n
measures the density of rod-like particles.As discussed in
the previous section,linear stability theory
predicts an instability and a wavelength selection mechani
sm for modulations of
&
.
For such a problem it is natural to ask how the aggregate respo
nse of the system is described in
long times.Our analysis of the di!usive scaling will show th
at it is described by the Keller-Segel model.
5.1.The hyperbolic limit.
With the goal to derive equations describing the aggregate b
ehavior
of the suspension,we turn to calculate the hyperbolic and di
!usive limits.
We first rescale the model (2.27) (with
Re
= 0) in the
hyperbolic scaling
,
x
=
1
(
ˆ
x
,t
=
1
(
ˆ
t,
u
=
ˆ
u
,p
= ˆ
p.
26
The scaled equations (after dropping the hats) are
($
t
f
+
(
u
∙ &
x
f
$
(D
(
n
)
e
2
∙ &
x
f
+
(
&
n

(
P
n
!
&
x
un
f
)
= $
n
f
+
(
2
'
&
x

D
(
n
)
&
x
f
$
(
2
$
x
u
+
(
&
x
p
$
(
2
'
&
x

%
=
$
(
!
&
S
2
fd
n
"
e
2
(
&
x

u
= 0
(5.1)
where
D
(
n
) =
I
+
n
%
n
.
We introduce the ansatz
f
=
(f
0
+
(
2
f
1
+
...
u
=
u
0
+
(
u
1
+
...
p
=
(p
0
+
(
2
p
1
+
...
to the system (5.1) and obtain equations for the various orde
rs of the expansion:
O
(
(
) $
n
f
0
= 0
(5.2)
O
(
(
2
)
$
t
f
0
+
u
0
∙ &
x
f
0
$
D
(
n
)
e
2
∙ &
x
f
0
+
&
n

(
P
n
!
&
x
u
0
n
f
0
) = $
n
f
1
(5.3)
O
(
(
2
)
$
$
x
u
0
+
&
x
p
0
=
$
!
&
S
2
f
0
d
n
"
e
2
O
(
(
)
&
x

u
0
= 0
It follows from (5.2) that
f
0
is independent of
n
and thus
f
0
(
t,
x
,
n
) =
1
4
*
&
S
2
f
0
d
n
=
1
4
*
&
0
(
t,
x
)
Then integrating (5.3) over the sphere,we deduce that
&
0
=
)
S
2
f
0
d
n
and
u
0
satisfy the Boussinesq
system
$
t
&
0
+
&
x

!
u
0
&
0
$
*
1
4
*
&
S
2
D
(
n
)
d
n
+
e
2
&
0
"
= 0
$
$
x
u
0
+
&
x
p
0
=
$
&
0
e
2
&
x

u
0
= 0
(5.4)
5.2.The di!usive limit.
Next,we consider rectilinear flows with a vertical velocity
field obeying
the ansatz
u
(
t,x,z
) = (0
,v
(
t,x,z
)
,
0)
T
,f
=
f
(
t,x,z
)
and depending only on the horizontal variables.This restri
ction - to the two-dimensional case - is
motivated by experimental observations of long clusters wi
th higher particle density.
Indeed,for flows as above the hydrodynamic scale produces a t
rivial equation and one may consider
the behavior in the di!usive scaling.We will show that in the
di!usive scale the aggregate behavior of
the system is described by the Keller-Segel system,which is
known to provide a nonlinear concentration
mechanism.
In the sequel we will monitor the quantities
f
(
t,
x
,
n
) =
ˆ
f
(
t,
x
+
t
e
2
1
4
*
&
S
2
D
(
n
)
d
n
,
n
)
u
(
t,
x
) =
ˆ
u
(
t,
x
+
t
e
2
1
4
*
&
S
2
D
(
n
)
d
n
)
(5.5)
27
Then
ˆ
f
and
ˆ
u
satisfy
$
t
ˆ
f
+
ˆ
u
∙ &
x
ˆ
f
$
!
D
(
n
)
$
1
4
*
&
S
2
D
(
n
)
d
n
"
e
2
∙ &
x
ˆ
f
+
&
n

'
P
n
!
&
x
ˆ
un
ˆ
f
(
= $
n
ˆ
f
+
'
&
x

D
(
n
)
&
x
ˆ
f.
$
$
x
ˆ
u
+
&
x
ˆ
p
$
'(
&
x

ˆ
%
=
$
(
!
&
S
2
ˆ
fd
n
"
e
2
(5.6)
We rescale the model using the
di!usive scaling
,i.e.
x
=
1
(
ˆ
x
,t
=
1
(
2
ˆ
t,
u
=
ˆ
u
,p
= ˆ
p.
(5.7)
The scaled equations (dropping the hats) have the form
(
2
$
t
f
+
(
u
∙ &
x
f
$
(
!
D
(
n
)
$
1
4
*
&
S
2
D
(
n
)
d
n
"
e
2
∙ &
x
f
+
(
&
n

(
P
n
!
&
x
un
f
) = $
n
f
+
(
2
'
&
x

D
(
n
)
&
x
f
$
(
2
$
x
u
+
(
&
x
p
$
(
2
'
&
x

%
=
$
(
!
&
S
2
fd
n
"
e
2
&
x

u
= 0
We introduce the ansatz
f
(
t,x,z,
n
) =
(f
0
+
(
2
f
1
+
...
u
(
t,x,z
) =
u
0
+
(
u
1
+
...
=
5
7
0
v
0
0
8
:
+
(
5
7
0
v
1
0
8
:
+
...
p
=
(p
0
+
(
2
p
1
+
...
to the above system and collect the terms of the same order,th
us arriving at
O
(
(
) $
n
f
0
= 0
(5.8)
O
(
(
2
)
$
!
D
(
n
)
$
1
4
*
&
S
2
D
(
n
)
d
n
"
e
2
∙ &
x
f
0
+
&
n

(
P
n
!
&
x
u
0
n
f
0
) = $
n
f
1
(5.9)
O
(
(
3
)
$
t
f
0
$
!
D
(
n
)
$
1
4
*
&
S
2
D
(
n
)
d
n
"
e
2
∙ &
x
f
1
*
P
n
!
*
&
x
u
1
n
f
0
+
&
x
u
0
n
f
1
++
= $
n
f
2
+
'
&
x

D
(
n
)
&
x
f
0
(5.10)
Here,we used the fact that due to the form of the ansatz all ter
ms of the type
u
i

(
&
x
f
)
j
= 0 vanish.
The same procedure applied to the Stokes system yields
O
(
(
2
)
$
$
x
u
0
+
&
x
p
0
=
$
!
&
S
2
f
0
d
n
"
e
2
O
(
(
)
&
x

u
0
= 0
Equation (5.8) implies that
f
0
is independent on
n
,that is
f
0
=
1
4
*
&
S
2
f
0
d
n
=
1
4
*
&
0
(
t,x,z
)
Next,integration of (5.10) over the sphere gives that
&
0
satisfies
$
t
&
0
=
&
x

&
S
2
!
D
(
n
)
$
1
4
*
&
S
2
D
(
n
)
d
n
"
e
2
f
1
d
n
0
12
3
=:
I
1
+
'
&
x

1
4
*
&
S
2
D
(
n
)
d
n
0
12
3
=:
I
2
&
x
&
0
(5.11)
28
where the terms
I
1
and
I
2
are computed in terms of
f
1
solving (5.9) for
f
0
=
1
4
#
&
0
.
It remains to compute the terms
I
1
and
I
2
.Observe now that we have the identities:
&
S
2
(
n
%
n
$
1
3
I
)
d
n
(B.4)
=
$
1
6
&
S
2
'
n
(
n
%
n
$
1
3
I
)
d
n
= 0
I
2
:=
1
4
*
&
S
2
D
(
n
)
d
n
=
1
4
*
&
S
2
(
n
%
n
+
I
)
d
n
=
4
3
I
(5.12)
D
(
n
)
$
1
4
*
&
S
2
D
(
n
)
d
n
=
n
%
n
$
1
3
I
These,in conjunction with (5.8),(5.9) and (B.7),imply tha
t
f
1
satisfies
$
n
f
1
=
$
*
D
(
n
)
$
1
4
*
&
S
2
D
(
n
)
d
n
+
e
2

1
4
*
&
x
&
0
+
1
4
*
&
n

(
P
n
!
&
x
u
0
n
)
&
0
=
$
*
n
%
n
$
1
3
I
+
e
2

1
4
*
&
x
&
0
$
3
4
*
&
0
(
n
∙ &
x
u
0
n
) (5.13)
Next,we compute
I
1
I
1
:=
&
S
2
!
D
(
n
)
$
1
4
*
&
S
2
D
(
n
)
d
n
"
e
2
f
1
d
n
=
&
S
2
*
n
%
n
$
1
3
I
+
e
2
f
1
d
n
(B.4)
=
$
1
6
&
S
2
'
n
*
n
%
n
$
1
3
I
+
e
2
f
1
d
n
(5.13)
=
1
24
*
&
S
2
*
n
%
n
$
1
3
I
+
e
2
.
*
n
%
n
$
1
3
I
+
e
2
∙ &
x
&
0
+3
&
0
(
n
∙ &
x
u
0
n
)
/
d
n
=
1
24
*
&
S
2
5
7
n
1
n
2
n
2
2
$
1
3
n
2
n
3
8
:
.
(3
&
0
v
0
x
+
&
0
x
)
n
1
n
2
+(3
&
0
v
0
z
+
&
0
z
)
n
2
n
3
/
d
n
Observe that,due to symmetry considerations,the integral
s
&
S
2
n
1
n
3
2
$
1
3
n
1
n
2
d
n
= 0
&
S
2
n
3
2
n
3
$
1
3
n
2
n
3
d
n
= 0
&
S
2
n
1
n
2
2
n
3
d
n
= = 0
,
while the remaining integrals are computed via spherical co
ordinates
&
S
2
n
2
1
n
2
2
d
n
=
&
#
0
sin
5
!d!
&
2
#
0
sin
2
0
cos
2
0d0
=
4
*
15
&
S
2
n
2
2
n
2
3
d
n
=
&
#
0
sin
3
!
cos
2
!d!
&
2
#
0
sin
2
0d0
=
4
*
15
.
We conclude that
I
1
=
1
90
(3
&
0
v
0
x
+
&
0
x
,
0
,
3
&
0
v
0
z
+
&
0
z
)
T
and that
&
0
satisfies the equation
$
t
&
0
=
1
30
!
$
x
!
1
3
&
0
x
+
&
0
v
0
x
"
+
$
z
!
1
3
&
0
z
+
&
0
v
0
z
""
+
'
4
3
$
(
x,z
)
&
0
=
1
30
&
(
x,z
)

*
&
0
&
(
x,z
)
v
0
+
+
1
3
!
4
'
+
1
30
"
'
(
x,z
)
&
0
29
We summarize the result:
Proposition 5.1.
Let
f
(
t,x,z,
n
)
,
u
(
t,x,z
) = (0
,v
(
t,x,z
)
,
0)
be a solution of
(5.6)
.In the
di!usive scaling
(5.7)
the leading order terms
f
0
=
1
4
#
&
0
(
t,x,z
)
and
v
0
(
t,x,z
)
satisfy the Keller-Segel
model
$
t
&
0
=
1
30
&
(
x,z
)

*
&
0
&
(
x,z
)
v
0
+
+
1
3
*
4
'
+
1
30
+
'
(
x,z
)
&
0
$
(
x,z
)
v
0
=
&
0
.
(5.14)
Conclusions.
We present a kinetic model for the sedimentation of dilute su
spensions of rod-like
Brownian particles in low Reynolds number flow.This model de
scribes concentration and the formation
of clusters with higher particle density.Such phenomena ha
ve previously been observed in experimental
and numerical studies of non-Brownian particles.
By restricting considerations to second moments,we derive
a linear model and study the linear
stability.We show that a non-zero Reynolds number leads to a
wavelength selection mechanism.This
indicates that inertial e!ects (which have often been exclu
ded in related mathematical models) can
influence the pattern formation on the macroscopic scale.
Furthermore,we show that at large times the macroscopic beh
avior of the system can be described
by a Keller-Segel equation.
Appendix A.Linear moment closure.
Here we give the details on the derivation of the two-
dimensional linear moment closure system used in Section 3.
Using the relations
sin
2
!
=
1
2
(1
$
cos(2
!
))
,
cos
2
!
=
1
2
(1 +cos(2
!
))
D
(
n
)
e
2
&
x
f
=
1
2
(
$
x
f
sin(2
!
) +
$
y
f
(3
$
cos(2
!
)))
&
n

(
P
n
!
&
x
un
) =
$
!
*
n
"
∙ &
x
un
+
=
$
!
!!
$
sin
!
cos
!
"
∙ &
x
u
!
cos
!
sin
!
""
= (
v
y
$
u
x
) cos(2
!
)
$
(
u
y
+
v
x
) sin(2
!
)
&
x

D
(
n
)
&
x
f
=
1
2
(3 +cos(2
!
))
$
xx
f
+sin(2
!
)
$
xy
f
+
1
2
(3
$
cos(2
!
))
$
yy
f
we rewrite the linearized Smoluchowski equation (3.1) into
the form
$
t
f
=
1
2
(sin(2
!
)
$
x
f
+(3
$
cos(2
!
))
$
y
f
)
0
12
3
(
I
)
+
1
2
*
(
u
x
$
v
y
) cos(2
!
) +
1
2
*
(
u
y
+
v
x
) sin(2
!
)
0
12
3
(
II
)
+
$
2
!
f
0123
(
III
)
+
'
!
1
2
(3 +cos(2
!
))
$
xx
f
+sin(2
!
)
$
xy
f
+
1
2
(3
$
cos(2
!
))
$
yy
f
"
0
12
3
(
IV
)
.
(A.1)
We consider a closure at the level of second moments
&
(
x
,t
):=
&
2
#
0
f
(
x
,!,t
)
d!,c
(
x
,t
):=
1
2
&
2
#
0
cos(2
!
)
f d!,s
(
x
,t
):=
1
2
&
2
#
0
sin(2
!
)
f d!
(A.2)
and derive evolution equations for
&
,
c
and
s
.
30
Evolution equation for
&
:.
The evolution of
&
can be expressed by
$
t
&
=
&
2
#
0
$
t
f d!.
We replace
$
t
f
by the terms on the right hand side of (A.1).Note that there is
no contribution from
(II) and (III).
Contribution from (I):
&
2
#
0
1
2
(sin(2
!
)
$
x
f
+(3
$
cos(2
!
))
$
y
f
)
d!
=
s
x
+
3
2
&
y
$
c
y
Contribution from (IV):
&
2
#
0
'
!
3 +cos(2
!
)
2
$
xx
f
+sin(2
!
)
$
xy
f
+
3
$
cos(2
!
)
2
$
yy
f
"
d!
=
'
!
3
2
(
&
xx
+
&
yy
) +2
s
xy
+
c
xx
$
c
yy
"
Thus we obtain
$
t
&
=
3
2
&
y
$
c
y
+
s
x
+
'
!
3
2
(
&
xx
+
&
yy
) +2
s
xy
+
c
xx
$
c
yy
"
(A.3)
Evolution equation for
c
:.
$
t
c
=
1
2
&
2
#
0
cos(2
!
)
$
t
f d!.
Contribution from (I):
1
4
&
2
#
0
cos(2
!
) (sin(2
!
)
$
x
f
+(3
$
cos(2
!
))
$
y
f
)
d!
=
1
2
&
2
#
0
sin(4
!
)
$
x
f d!
+
3
2
c
y
$
1
4
&
2
#
0
cos
2
(2
!
)
0
12
3
=
1
2
(1+cos(4
!
))
$
y
f d!
˙=
3
2
c
y
$
1
8
&
y
In this approximation (notation ˙=) we neglect terms that in
volve integrals of the form
)
2
#
0
sin(4
!
)
q d!
and
)
2
#
0
cos(4
!
)
q d!
,respectively.
Contribution from (II):
1
2
1
2
*
&
2
#
0
cos(2
!
) [(
u
x
$
v
y
) cos(2
!
) +(
u
y
+
v
x
) sin(2
!
)]
d!
=
1
2
1
2
*
5
6
7
&
2
#
0
(
u
x
$
v
y
) cos
2
(2
!
)
0
12
3
=
1
2
(1+cos(4
!
))
d!
+
&
2
#
0
(
u
y
+
v
x
)
1
2
sin(4
!
)
d!
8
9
:
˙=
1
4
(
u
x
$
v
y
)
Contribution from (III):
1
2
&
2
#
0
cos(2
!
)
*
$
2
!
f
+
d!
=
1
2
&
2
#
0
2sin(2
!
)
$
!
f d!
=
$
&
2
#
0
2cos(2
!
)
f d!
=
$
4
c
31
Contribution from (IV):
1
2
'
&
2
#
0
cos(2
!
)
$
3 +cos(2
!
)
2
$
xx
f
+sin(2
!
)
$
xy
f
+
3
$
cos(2
!
)
2
$
yy
f
%
d!
=
'
.
3
2
c
xx
+
1
4
&
2
#
0
cos
2
(2
!
)
0
12
3
=
1
2
(1+cos(4
!
))
$
xx
f d!
+
1
4
&
2
#
0
sin(4
!
)
$
xy
f d!
+
3
2
c
yy
$
1
4
&
2
#
0
cos
2
(2
!
)
$
yy
f d!
/
˙=
'
!
3
2
(
c
xx
+
c
yy
) +
1
8
(
&
xx
$
&
yy
)
"
We obtain:
$
t
c
=
$
1
8
&
y
+
3
2
c
y
+
1
4
(
u
x
$
v
y
)
$
4
c
+
'
'
3
2
(
c
xx
+
c
yy
) +
1
8
(
&
xx
$
&
yy
)
(
(A.4)
Evolution equation for
s
:.
$
t
s
=
1
2
&
2
#
0
sin(2
!
)
$
t
f d!.
Contribution from (I):
1
4
&
2
#
0
sin(2
!
) [sin(2
!
)
$
x
f
+(3
$
cos(2
!
))
$
y
f
]
d!
=
1
4
&
2
#
0
sin
2
(2
!
)
0
12
3
=
1
2
(1
!
cos(4
!
))
$
x
f d!
+
3
2
s
y
$
1
8
&
2
#
0
sin(4
!
)
$
y
f d!
˙=
1
8
&
x
+
3
2
s
y
Contribution from (II):
1
2
1
2
*
&
2
#
0
sin(2
!
) [(
u
x
$
v
y
) cos(2
!
) +(
u
y
+
v
x
) sin(2
!
)]
d!
=
1
4
1
2
*
(
u
x
$
v
y
)
&
2
#
0
sin(4
!
)
d!
+
1
2
1
2
*
(
u
y
+
v
x
)
&
2
#
0
sin
2
(2
!
)
d!
˙=
1
4
(
u
y
+
v
x
)
Contribution from (III):
1
2
&
2
#
0
sin(2
!
)(
$
2
!
h
)
d!
=
$
&
2
#
0
cos(2
!
)
$
!
f d!
=
$
2
&
2
#
0
sin(2
!
)
f d!
=
$
4
s
Contribution from (IV):
1
2
'
&
2
#
0
sin(2
!
)
$
1
2
(3 +cos(2
!
))
$
xx
f
+sin(2
!
)
$
xy
f
+
1
2
(3
$
cos(2
!
))
$
yy
f
%
d!
˙=
'
$
3
2
s
xx
+
1
2
&
2
#
0
sin
2
(2
!
)
$
xy
f d!
+
3
2
s
yy
%
˙=
'
!
3
2
(
s
xx
+
s
yy
) +
1
4
&
xy
"
32
We obtain
$
t
s
=
1
8
&
x
+
3
2
s
y
+
1
4
(
u
y
+
v
x
)
$
4
s
+
'
!
1
4
&
xy
+
3
2
(
s
xx
+
s
yy
)
"
(A.5)
Appendix B.Properties of the di!erential operators in spher
ical coordinates.
The operator
&
n
satisfies certain elementary properties that are extensive
ly used in this article:
Let
F
be a vector-valued function and
f
,
g
be scalar-valued functions,then
&
S
2
(
&
n

F
)
fd
n
=
$
&
S
2
F

(
&
n
f
$
2
n
f
)
d
n
(B.1)
&
S
2
(
&
n
∙ &
n
f
)
gd
n
=
&
S
2
(
&
n
∙ &
n
g
)
fd
n
(B.2)
&
S
2
n
%&
n
fd
n
=
&
S
2
&
n
f
%
n
d
n
=
&
S
2
(3
n
%
n
$
id)
fd
n
(B.3)
The components of the tensor 3
n
%
n
$
id are the surface spherical harmonics of order 2.That is
they are harmonic polynomials on
R
3
of order 2,restricted to
S
2
.The surface spherical harmonics are
eigenfunctions of the Laplacian on
S
2
with corresponding eigenvalue
$
1
(
1
+1),where
1
is the order [1,
App.E].Hence
'
n
(3
n
i
n
j
$
(
ij
) =
$
6(3
n
i
n
j
$
(
ij
)
.
(B.4)
It is convenient to use spherical coordinates in proving suc
h formulas,see [1,App.A.6 and E.6].
A point
P
with Cartesian coordinates (
n
x
,n
y
,n
z
) is expressed in spherical coordinates via
n
x
=
r
sin
!
cos
0,n
y
=
r
sin
!
sin
0,n
z
=
r
cos
!
where 0
<!< *
,0
+
0 <
2
*
.Let
e
r
,
e
!
,
e
'
be the orthonormal coordinate system associated to
spherical coordinates and attached at
P
.It satisfies the derivative formulas
(
e
r
(r
= 0
,
(
e
r
(!
=
e
!
,
(
e
r
('
=
e
'
sin
!,
(
e
!
(r
= 0
,
(
e
!
(!
=
$
e
r
,
(
e
!
('
=
e
'
cos
!,
(
e
"
(r
= 0
,
(
e
"
(!
= 0
,
(
e
"
('
=
$
e
r
sin
!
$
e
!
cos
!.
(B.5)
We visualize the sphere
S
2
as embedded in the Euclidean space.The surface gradient
&
n
is related
to the (ambient) gradient operator
&
through
&
n
=
r
(id
$
n
%
n
)
∙ &
=
e
!
$
$!
+
e
'
1
sin
!
$
$0
and,for a scalar-valued function
f
,the surface gradient and the Laplace-Beltrami operator ar
e given
respectively by
&
n
f
=
e
!
$f
$!
+
e
'
1
sin
!
$f
$0
'
n
f
=
&
n
∙ &
n
f
=
1
sin
!
$
$!
!
sin
!
$f
$!
"
+
1
sin
2
!
$
2
f
$0
2
The divergence of a vector-valued function
F
=
F
r
e
r
+
F
!
e
!
+
F
'
e
'
has the form
&
n

F
=
!
e
!
$
$!
+
e
'
1
sin
!
$
$0
"

*
F
r
e
r
+
F
!
e
!
+
F
'
e
'
+
(B.5)
=
1
sin
!
$
$!
(sin
!F
!
) +
1
sin
!
$F
'
$0
+2
F
r
33
It is now easy to compute (B.1)-(B.3).Observe that
&
S
2
(
&
n

F
)
fd
n
=
&&
!
1
sin
!
$
$!
(sin
!F
!
) +
1
sin
!
$F
'
$0
+2
F
r
"
f
sin
!d!d0
=
$
&&
!
$
2
F
r
f
+
F
!
$f
$!
+
1
sin
!
F
'
$f
$0
"
sin
!d!d0
=
$
&
S
2
F

(
&
n
f
$
2
n
f
)
d
n
which gives (B.1),and (B.2) follows by applying (B.1) twice
:
&
S
2
(
&
n
∙ &
n
f
)
gd
n
=
$
&
S
2
&
n
f

(
&
n
g
$
2
n
g
)
d
n
=
&
S
2
f
(
&
n
∙ &
n
g
)
d
n
.
Using integration by parts,we obtain the chain of identitie
s
&
S
2
n
%&
n
fd
n
=
&&
e
r
%
'
e
!
$f
$!
+
e
'
1
sin
!
$f
$0
(
sin
!d!d0
=
$
&&
.
$
$!
(
e
r
%
e
!
) +
cos
!
sin
!
e
r
%
e
!
+
1
sin
!
$
$0
(
e
r
%
e
'
)
/
f
sin
!d!d0
(
B.
5)
=
$
&&
.
e
!
%
e
!
+
e
'
%
e
'
$
2
e
r
%
e
r
/
f
sin
!d!d0
=
&
S
2
(3
n
%
n
$
id)
fd
n
(B.6)
Since the final equation in (B.6) is a symmetric tensor,(B.3)
follows.
Note that,for any 3
(
3 matrix
2
,the equation holds
&
n

(
P
n
!
2
n
) = tr
2
$
3
n

2
n
,
(B.7)
where tr stands for the trace operator.Indeed,using (B.5) w
e obtain
&
n

(
P
n
!
2
n
) =
!
e
!
$
$!
+
e
'
1
sin
!
$
$0
"

*
2
e
r
$
(
e
r

2
e
r
)
e
r
+
(
B.
5)
=
e
!

2
e
!
+
e
'

2
e
'
$
2
e
r

2
e
r
= tr
2
$
3
n

2
n
(B.8)
Acknowledgements.
This joint research was partially supported by the ACMAC pro
ject of the
FP7-REGPOt-2009-1 program of the European Commission.FO p
artially supported by the German
Science Foundation through SFB 611.AET partially supporte
d by the National Science Foundation.
REFERENCES
[1] R.B.Bird,Ch.F.Curtiss,R.C.Armstrong,and O.Hassager.
Dynamics of Polymeric Liquids,Vol.2,Kinetic
Theory
.Wiley Interscience,1987.
[2] J.E.Butler and E.S.G.Shaqfeh.Dynamic simulation of the inh
omogeneous sedimentation of rigid fibers.
J.Fluid
Mech.
,468:205–237,2002.
[3] R.E.Caflish and J.H.C.Luke.Variance in the sedimentation speed of
a suspension.
Phys.Fluids
,28:759–760,1985.
[4] F.A.C.C.Chalub,P.A.Markowich,B.Perthame,and C.Schmeiser
.Kinetic models for chemotaxis and their drift-
di!usion limits.
Monatsh.Math.
,142:123–141,2004.
[5] M.Doi and S.F.Edwards.
The Theory of Polymer Dynamics
.Oxford University Press,1986.
[6] J.Dolbeault,P.Markowich,D.Oelz,and C.Schmeiser.Non linear
di!usions as limit of kinetic equations with
relaxation collision kernels.
Arch.Rational Mech.Anal.
,186:133–158,2007.
[7] M.Feist,H.Nirschl,J.Wagner,G.Hirsch,and S.Schabel.Experi
mental results for the settling behaviour of
particle-fiber mixtures.
Physical Separation in Science and Engineering
,2007.
34
[8] K.Gustavsson and A.K.Tornberg.Gravity induced sedimentation of
slender fibers.
Phys.Fluids
,21:123301,2009.
[9] C.Helzel and F.Otto.Multiscale simulations for suspensions of rod
-like molecules.
J.Comput.Phys.
,216(1):52–75,
2006.
[10] B.Herzhaft and
´
E.Guazzelli.Experimental study of the sedimentation of a dilute fiber su
spension.
Physical Review
Letters
,77:290–293,1996.
[11] B.Herzhaft and
´
E.Guazzelli.Experimental study of the sedimentation of dilute and semi
-dilute suspensions of
fibers.
J.Fluid Mech.
,384:133–158,1999.
[12] D.L.Koch and E.S.G.Shaqfeh.The instability of a dispersion
of sedimenting spheroids.
J.Fluid Mech.
,209:521–542,
1989.
[13] J.H.L.Luke.Decay of velocity fluctuations in a stably stratified su
spension.
Phys.Fluids
,12:1619–1621,2000.
[14] B.Metzger,J.E.Butlerz,and E.Guazzelli.Experimental investigatio
n of the instability of a sedimenting suspension
of fibers.
J.Fluid Mech.
,575:307–332,2007.
[15] B.Metzger,J.E.Butlerz,and E.Guazzelli.On wavelength selection by stra
tification in the instability of settling
fibers.
Phys.Fluids
,19:098105,2007.
[16] F.Otto and A.Tzavaras.Continuity of velocity gradients in su
spensions of rod-like molecules.
Comm.Math.Phys.
,
277(3):729–758,2008.
[17] D.Saintillan,S.Darve,and S.G.Shaqfeh.Asmooth particle-
mesh ewald algorithmfor stokes suspension simulations:
The sedimentation of fibers.
Phys.Fluids
,17:033301,2005.
[18] D.Saintillan,S.G.Shaqfeh,and E.Darve.The e!ects of strati
fication on the wave number selection in the instability
of sedimenting spheroids.
Phys.Fluids
,18:121503,2006.
[19] D.Saintillan,S.G.Shaqfeh,and E.Darve.The grows of concen
tration fluctuations in dilute dispersions of orientable
and deformable particles under sedimentation.
J.Fluid Mech.
,553:347–388,2006.
[20] D.Saintillan and M.J.Shelley.Orientational order and inst
ability in suspensions of self-locomoting rods.
Phys.Rev.
Lett.
,99:058102,2007.
[21] A.K.Tornberg and K.Gustavsson.A numerical method for simulati
ons of rigid fiber suspensions.
J.Comput.Phys.
,
215:172–196,2006.
[22] J.Wang and A.Layton.Numerical simulations of fiber sedimentat
ion in Navier-Stokes flows.
J.Comput.Phys.
,
5:61–83,2009.