A KINETIC MODEL FOR THE SEDIMENTATION OF ROD–LIKE PARTICLES
CHRISTIANE HELZEL
!
,FELIX OTTO
†
,
AND
ATHANASIOS E.TZAVARAS
‡
Abstract.
We consider a kinetic multiscale model,describing suspensions of rodlik
e particles,which couples a
microscopic Smoluchowski equation  for the distribution of rod or
ientations  to a macroscopic Stokes or NavierStokes
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 conﬁrmthe wavelength selection.
By looking at special ﬂow conﬁgurations 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 twodimensional
KellerSegel model.
Key words.
sedimentation,rodlike particles,linear stability,asymptotic exp
ansions
AMS subject classiﬁcations.
76T20,74A25,35B40
1.Introduction.
Complex ﬂuids are ﬂuids with a suspended microstructure suc
h as polymers
with a long chain structure.A simpliﬁed model is that of susp
ensions of rodlike particles.In such
systems multiscale interactions can cause interesting phe
nomena.The macroscopic ﬂow leads to a
change of the orientation and shape (in the case of ﬂexible pa
rticles) of the suspended microstructure.
In turn,this produces a ﬂuid stress which interacts with and
modiﬁes the macroscopic ﬂow.
Here we consider the sedimentation of dilute suspensions of
rodlike 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 selflocomotion [20].In
contrast to spherical particles,rodlike 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 nonBrownian 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 rodlike molecules neglecting the e!ects of
gravity,see [9],[16].
The sedimentation of ﬁbers 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 wellstirred
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 ﬂip.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 rodlike particles (an
d other orientable particles) has also been
studied via numerical simulations of multiscale models.G
ustavsson and Tornberg [8,21] used a very
detailed description of rodlike 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,RuhrUniversity Bochum,Germany
christiane.helzel@rub.de
,
†
Institute of Applied Mathematics,University of Bonn,Germany
otto@iam.unibonn.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 twodimensional numerical studies.All numerical re
sults conﬁrm the basic experimental ﬁndings:
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 nonBrownian sedimentation of spherical part
icles,there is a mathematically related
divergence of ﬂuctuations at long wavelength [3].For that p
hysical system,a small stable density
stratiﬁcation has been proposed as a mechanism which suppre
sses this long wavelength divergence [13].
In the case of sedimenting rodlike particles it also leads t
o a wavelength selection at the level of linear
stability analysis [18].In experiments and numerical simu
lations,this stratiﬁcation 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 rodl
ike particles.We ﬁrst present an analysis of
the thermodynamic consistency of the kinetic model and purs
ue a nondimensional analysis.There are
three nondimensional parameters in the problem.The Reyno
lds number which quantiﬁes the relative
importance of inertial forces versus viscous forces and two
additional nondimensional 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 inﬂuence 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 ﬂow is modeled by the Stokes equation,lin
ear stability theory does not predict a
wavelength selection mechanism.However,if the macroscop
ic ﬂow is modeled by the NavierStokes
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 nonBrownian particles.
While we do not have a theoretical justiﬁcation 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 ﬂow conﬁguration and
a di!usive scaling of the equations the
density follows a nonlinear concentration mechanismthat i
s captured by the wellknown twodimensional
KellerSegel 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 multiscale 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 rodlike 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 rodlike pa
rticles was considered.Here we extend this
model to account for the e!ects of gravity.
We consider inﬂexible rodlike 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 rodlike 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 rodlike
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 ﬂow
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 ﬂow
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 ﬂow do
main are denoted by
&
x
and
&
x
∙
.The total number of rodlike 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 speciﬁc weight of the
suspension that generally can not be compensated by a hydros
tatic pressure and thus trigger a ﬂuid
motion (buoyancy).The macroscopic ﬂow is described by the N
avierStokes equation.Let
&
f
be the
density of the ﬂuid 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 redeﬁning 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 rodlike particles.
We summarize the ﬁnal 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
]
deﬁned
in (2.11) satisﬁes 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 NavierStokes 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.Nondimensionalization.
We ﬁrst 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
nondimensionalization 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 ﬂow.
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 ﬂuid
$
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 nondimensionalization of the elastic stress tensor le
ads to the expression
ˆ
%
=
&
S
d
"
1
(
d
n
%
n
$
I
)
ˆ
fd
n
.
(2.24)
Finally,we nondimensionalize the NavierStokes equatio
n.We obtain
X
2
Tµ
&
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 nondimensional numbers emerging:First,a
Reynolds number based on the microscopic
velocity
X/T
and the microscopic length scale
X
,
Re
:=
&
f
X
2
Tµ
=
&
f
µ
!
m
0
g
#
"
"
2
1
D
r
.
Second the nondimensional 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 nondimensional 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 nondimensional form of the NavierStokes 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 nondimensional 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 NavierStokes equation in the equivalent f
orm (2.6),then the nondimensionalized
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.Multiscale mechanism for instability and cluster form
ation.
The multiscale mecha
nismthat leads to the instability and the formation of clust
ers was ﬁrst 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 ﬂow with ﬂow 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 ﬂow 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 ﬁgure 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 nonvan
ishing angle between ﬂow 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 ﬁeld,(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
rodlike 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 inﬁnitesimal growth rate is independent of
the wave number

!

in the limit

!
 *
0.Saintillan et al.[18] showed that a small density stratiﬁ
cation
can lead to a ﬁnite 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 ﬂow.
Here we show that the use of the NavierStokes equation leads
to a wavelength selection mechanism
which depends on the dimensionless parameters
Re
,
(
and
'
.
In this section we restrict considerations to the spatially
twodimensional case (one horizontal
direction + direction of gravity).In this conﬁguration 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 ﬁrst
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 deﬁned 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 ﬁne 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 micromacro 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 NavierStokes equation.
We now study the linear stability for the system (3.3)(3.5)
together with the NavierStokes 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 nonmovin
g frame (instead of transforming the
NavierStokes 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
)
$
4ˆ
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
)
$
4ˆ
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 micromacro model with the NavierStokes 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 conﬁguration
.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 speciﬁed 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 ﬁts 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 KellerSegel model.
4.Sedimentation of nonBrownian 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.
nonBrownian 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 nondimensionalize 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 nondimensional 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 nondimensional 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 justiﬁcation
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 twodimensional 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 deﬁned 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 ﬁrst 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 nonBrownian particles in a Stokes ﬂ
ow (with
$
= 0
.
01
,
0
.
1
,
1
).
no wavelength selection mechanism for the micromacro 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
vierStokes 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 nonBrownian 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 nonBrownian particles in a ﬂow 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 ﬂuid with ﬂuid 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 speciﬁed 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 onedimensional nonlinear model for
$
= 1
and
Re
= 1
.
5.The hyperbolic and di!usive scalings.Relation to the Kel
lerSegel 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 ﬂows the
hyperbolic scaling produces a trivial behavior,and it is th
en natural to consider the di!usive scaling.
Such a situation occurs for twodimensional rectilinear ﬂo
ws of suspensions.We will show that the
collective behavior in the di!usive limit is described by th
e KellerSegel 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 rodlike 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 KellerSegel 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 ﬁrst 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 ﬂows with a vertical velocity
ﬁeld 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 twodimensional case  is
motivated by experimental observations of long clusters wi
th higher particle density.
Indeed,for ﬂows 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 KellerSegel 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
satisﬁes
$
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
satisﬁes
$
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
satisﬁes 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 KellerSegel
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 rodlike
Brownian particles in low Reynolds number ﬂow.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 nonBrownian particles.
By restricting considerations to second moments,we derive
a linear model and study the linear
stability.We show that a nonzero 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
inﬂuence 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 KellerSegel 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
satisﬁes certain elementary properties that are extensive
ly used in this article:
Let
F
be a vectorvalued function and
f
,
g
be scalarvalued 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 satisﬁes 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 scalarvalued function
f
,the surface gradient and the LaplaceBeltrami 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 vectorvalued 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 ﬁnal 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
FP7REGPOt20091 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 ﬁbers.
J.Fluid
Mech.
,468:205–237,2002.
[3] R.E.Caﬂish 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ﬁber mixtures.
Physical Separation in Science and Engineering
,2007.
34
[8] K.Gustavsson and A.K.Tornberg.Gravity induced sedimentation of
slender ﬁbers.
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 ﬁber 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
ﬁbers.
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 ﬂuctuations in a stably stratiﬁed 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 ﬁbers.
J.Fluid Mech.
,575:307–332,2007.
[15] B.Metzger,J.E.Butlerz,and E.Guazzelli.On wavelength selection by stra
tiﬁcation in the instability of settling
ﬁbers.
Phys.Fluids
,19:098105,2007.
[16] F.Otto and A.Tzavaras.Continuity of velocity gradients in su
spensions of rodlike 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 ﬁbers.
Phys.Fluids
,17:033301,2005.
[18] D.Saintillan,S.G.Shaqfeh,and E.Darve.The e!ects of strati
ﬁcation 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 ﬂuctuations 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 selflocomoting rods.
Phys.Rev.
Lett.
,99:058102,2007.
[21] A.K.Tornberg and K.Gustavsson.A numerical method for simulati
ons of rigid ﬁber suspensions.
J.Comput.Phys.
,
215:172–196,2006.
[22] J.Wang and A.Layton.Numerical simulations of ﬁber sedimentat
ion in NavierStokes ﬂows.
J.Comput.Phys.
,
5:61–83,2009.
Σχόλια 0
Συνδεθείτε για να κοινοποιήσετε σχόλιο