7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Motion of Heavy Ellipsoids Under Linear Shear in Stokes Flow
F. Lundell* and A. Carlssont
Linn6 Flow Centre, KTH Mechanics, Royal Institute of Technology, S100 44 Stockholm, Sweden
t Department of Mechanical Engineering, University of British Columbia, Vancouver, British Columbia V6T 1Z4, Canada
fredrik@mech.kth.se and allanfc@ interchange.ubc.ca
Keywords: Ellipsoid, Shear Flow, Stokes flow, Particle inertia
Abstract
Motions of ellipsoids with nonnegligible particle inertia in linear shear flow is studied under Stokes flow conditions.
This constitutes a fundamental system that is used as a basis for simulations and analysis of flows with heavy
nonspherical particles. The torque on the ellipsoid is given analytically Jeffery [Proc. Roy. Society 102 (715),
1922] and this torque is coupled with the equations of motion for the ellipsoid. The resulting equations are studied
theoretically in terms of stability of the periodic solutions representing rotation around the principal axes of the
ellipsoid. Furthermore, particle motions are obtained numerically. It is known that rotation around the shortest axis
is inherently unstable for inertiafree prolate triaxial ellipsoids and that this instability leads to chaotic motions.
We show that this instability is stabilized at a critical level of particle inertia (the exact level depends on particle
geometry). The simulations of particle motion reveal that particle inertia induces a drift towards rotation around
the shortest axis for all ellipsoids studied. As this motion is approached, the behaviour depends on (i) whether the
ellipsoid is prolate and inherently unstable or not and (ii) on the level of particle inertia. If the particle geometry is
inherently stable, it eventually rotates around the shortest axis, which is aligned with the vorticity axis. If its geometry
is inherently unstable, the particle enters a chaotic motion for particle inertia below the critical one. For supercritical
particle inertia, the motion is similar to that of inherently stable particles. Quantitative results for particles with an
aspect ratio of 10 between the longest and shortest axis are presented.
Introduction
Ellipsoids are often used as a model for nonspherical
particles such as flakes or fibers, which appear in sev
eral biological, environmental and processing situations.
Here, motions of ellipsoids in linear shear flow is stud
ied. The basis for modelling the particle motion is the
analytical work by Jeffery (1922) who determined the
torque in linearly varying Stokes flow for all ellipsoid
axes lengths, orientations and rotation velocities. He
also determined the motion of inertiafree spheroids (el
lipsoids with two axis of equal length) by setting the
torque to zero and study the resulting equations. This
results in the socalled Jeffery orbits, which is a family
of kayaking motions.
For the present analysis, we introduce two coordinate
systems, an inertial system x, y, z and one fixed in the el
lipsoid x', y', z'. The ellipsoid is assumed to be defined
by
'/2 '2
'2 +Y k2
x
' '
()2
2 } '
where 1 is the length of the longest axis of the ellipsoid.
The external flow driving the particle motion is
u = y.
where u is the velocity in the xdirection and j is the
shear rate.
The motion of the particle motion is then governed by
(i) the ellipsoid geometry defined by kb and ke, (ii) the
Reynolds number Re# = 12/v, measuring the effect of
fluid inertia, and (iii) the Stokes number St = KRe ,
where K = pe/pf where pe and pf are the densities of
the particle and fluid, respectively. The Stokes number
measures the effect of particle inertia.
The effects of particle and fluid inertia in this system,
mainly on spheroids, has been studied extensively with
different methods by e.g. Ding & Aidun (2000); Qi &
Luo (2003); Subramanian & Koch (2006); Altenbach
et al. (2007); Lundell & Carlsson (2010). When it
comes to triaxial ellipsoids (kb / ke), the literature
is not as rich but important contributions are given by
Hinch & Leal (1979) and Yarin, Gottlieb & Roisman
(1997), who both studied the stability of the motion of
triaxial ellipsoids. Considering the equations of Jeffery
(1922), Hinch & Leal (1979) showed that the rotation of
prolate (fiber like) triaxial particles around the smallest
axis and oblate (flake like) triaxial particles around the
longest axis can be unstable if the two smallest prolatee)
or largest (oblate) axes are close enough in length. In the
present parametrization, this means that rotation with the
z and z' axes aligned is unstable for kc < kb < k"1,c
c,cit 
and k1Jr < kc < 1 in the limit of zero Re and St.
The critical values kIC o and k~ i are functions of kc
and kb, respectively. Yarin, Gottlieb & Roisman (1997)
showed that this instability leads to chaotic motions of
prolate triaxial particles if its initial orientation is close
to the unstable one.
Effects of nonzero particle inertia on prolate ellip
soids was studied by Lundell & Carlsson (2010) and
drastic effects, both on the rotation rate and the motion in
space was found. Here, that work is continued by study
ing how particle inertia affects the instability of prolate
triaxial ellipsoids. First, particle inertia is found to be
stabilizing. Consequences for particle motion are then
determined by numerical integration of the equations of
motion.
Models
The numerical models of the present work are the same
as the ones used by Lundell & Carlsson (2010) and are
obtained by coupling the law of angular momentum for
the particle with the torques given by Jeffery (1922).
The models are summarized below. Two models are
used: one for rotation with the z and z' axes aligned and
one for free threedimensional rotation of the particle.
Rotation With the z and z' Axes Aligned. For this
case, the angle p is introduced as the angle between the
x and x' axis. The law of angular momentum for rotation
around a fixed axis together with the torque of Jeffery
(1922) then gives an equation of motion for p as:
S k(kbKo + K)..
20
1i+ 1
1+b (2
sin2 p
(1 ) (1)
2>+ 1
where the dot is differentiation with respect to time,
Ko K(kk, f), K K(kb2, kb2) and
K(Cl, C2)
dA
= __(2)
S(1 + A){(1 + A)(C1 + A)(C2 + )}/2 (2)
In the present work, this model will be used to obtain
periodic solutions for which the stability is determined
by Floquet analysis (to be described later).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Rotation in Three Dimensions.
Shivarama & Fahrentgold (2'k14) derived a
constraintfree quaternionbased formulation of the
equations of motion for a rotating (and translating) body
and the rotation part of this formulation will be used
below. The orientation of a body is now described by a
quaternion, i.e. a vector of four numbers with norm 1,
e = [eo, e, e, es]T.
The physical interpretation of the quaternion is that if
the body has been rotated by an angle y around an axis
defined by a unit vector in the inertial frame of reference,
b' [b[, b', b'], the corresponding quaternion is
e [cos(q/2), bW sin(q/2), b' sin(q/2), Wb sin(q/2)]T.
The rotation matrix R relating a vector a' in the inertial
system to the components in the body fixed system, a
(i.e. a' Ra), is then
R EGT
where
eI e[ e3 e2
E = e2 e3 e0 e1
e3 e2 eI eo
eC eC e3 e2
G e2 e3 6o 61 (6)
e3 e2 eC e0
Now introduce the angular momentum vector h = Jc,
where J and cw are the momentofinertia matrix and ro
tational velocity vector, respectively, in the bodyfixed
coordinates. The nondimensional equations of motion
for an ellipsoid in creeping shear can then be written as:
1 16T
e GTw, h + 16T
2 3St
where Q 2GGT and 16iT/3St is the non
dimensional torque, obtained from Jeffery (1922):
(Kakb + Kc) 1x
[(k k)D2,3 + (k + k)(V2,3
(KcKc + Ka)1x
[(2 1)D3,1+ (2 + 1)(V3,1 
(K, + Kkb)lX
[(1 k)D1,2 + (1+ )(Vi,2 
 C1)]
C2)]
c3)]
(8)
where K K(k 2, k 2 ) and D, V are the rates of
deformation and vorticity in the body fixed system.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
This model for threedimensional rotation will be
used to (i) determine the stability of periodic solutions
obtained with the fixed axis model and (ii) to study the
motion of selected triaxial ellipsoids.
Floquet Analysis. The Floquet analysis (e.g. Joseph
(1976)) is performed for periodic solutions of p with
period T, p(t) = (t + T) and t > to. The value of
to, at which the solutions have reached periodicity after
initial transients due to inertia, is determined by the cri
terion that the time for the last five periods must differ
less than 0.05%. The stability of these solutions is then
studied by writing equation 7 as
x f(x)
where x = [C1, C2,3,, C, C1, e 2, Ce]T and derive the
following matrix equation:
X Of((t), (t)) (p(t))X. (9)
Ox
Note that x(t) = [0, 0,, cos(p/2), 0, 0, sin(p/2)]T
in this case. Equation () is integrated for one period of
x, which means that the integration is performed from
to to to + 2T with the initial condition X(to) = I. Flo
quet theory then states that the periodic motion given by
p(t) p(t + T) is then unstable if
A = eig[X(to + 2T)] ,,, > 1.
Results
Stability of Rotation Around Shortest Axis. In figure
1, contours of the amplitude of the largest eigenvalue
of X(to + 2T) is shown for kc 0.1 and different kb
and St. The contours start at 2 and increases towards
the peak. The contour at A = 2 is shown rather then
the nominal stability limit of A = 1, since the stable
solutions gives A = 1. Due to numerical inaccuracies,
this value is often slightly higher than 1 why the contour
A = 1 is meaningless. However, the gradient of A from
the unstable region towards the stable region is large,
why the somewhat arbitrarily chosen contour A = 'E I
very close to the stability boundary.
For light particles (low St), the range 0 < kb < 0.24
is unstable, in good agreement with the results of Yarin,
Gottlieb & Roisman (1997). As St is increased, A de
creases and eventually the periodic solutions obtained
from equation 1 is stabilized. Stabilization sets in at high
kb, but also the major portion of the unstable range of kb
is stabilized for St in the range 100 to 250 for the par
ticular choice of kc = 0.1.
Ellipsoid Motions. The effect of this instability and sta
bilization will now be demonstrated by looking at nu
merical solutions of equation 7 for the parameter cases
Figure 1: Magnitude of the largest Floquet eigenvalue A
as a function of kb and St for kc = 0.1. The
contours are, from blue to red: 2, 6, 15, 35,
55, 75, 95 and 105. Particle motions for the
parameter combinations marked with 0 are
shown in figures 3 6 while the ones shown
with x are shown in the appendix (figures 8
13). The markers have the same colours the
respective curves in the figures to come.
marked with in figure 1. For completeness, the cases
marked with x are provided in the appendix.
The cases studied in detail are kc = 0.1 and kb 0.14
and 0.22 and St 1, 10, 100 and 1000. The initial con
dition is that the particle is at rest with the x' and z' axes
in the xz plane and an angle of 15 between the x and x'
axes. The initial configurations are illustrated in figure 2.
The trace of the z'axis is shown in figures 3 (kb 0.14)
and 5 (kb = 0.22). The coloured 3D curves shows the
trace for increasing time and the projections on the co
ordinate planes are shown with grey curves. The colours
are the same as the respective marker in figure 1. In fig
ures 4 (kb 0.14) and 6 (kb 0.22) the zvalue of the
x'axis (denoted xz) is shown as a function of time for
the Stvalues studied. If xz 0, the longest axis of the
ellipsoid (the z'axis) is in the .r/plane
It is seen that for low St and the lower value of kb,
figure 3 (a,b), there is a kayaking motion with a trend
of xz towards zero (the blue and red curves in figure 4).
For higher St, 3 (c) and the green curve in figure 4 the
extremes of xz decreases initially but eventually starts to
vary intermittently from positive to negative values. This
motion is similar to the chaotic motion of light particles
found by Yarin, Gottlieb & Roisman (1997). The trace
of the end point (figure 3 c) is kayaking with varying
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
 
 4
0
1
0
1 1
x
(a) 1
y 0
(a) 1
0.5
0
0.5
1
1
(b) 1
Y 0
1
1
(c) 1
S0
Figure 2: Ellipsoid, shear flow and initial configuration
for the results shown in figures 3 6. kc = 0.1
and (a): kb = 0.14, (b): kb = 0.22. The angle
between the x' and xaxes is 15. Projections
of the axes of the ellipsoids are shown with
red lines.
(d) 1 
amplitude on both sides of z = 0 in a chaotic manner.
At St = 1000, in the stable region in figure 1, figure
3 (d) shows that there is no chaotic motion. Instead, the
ellipsoid rotates around a fixed axis that approaches the
zaxis, similar to the motions found for prolate spheroids
at high St byLundell & Carlsson (2010).
For the particle with kb 0.22 and kc 0.1 studied
in figures 5 and 6, the motions are similar. However, the
motion is stable already at St 100 as seen in figure 5
(c) and the green curve with a rapid approach to xz = 0
in figure 6. This stabilization was predicted in figure 1.
y 0
1
1
0
1 1
x
Figure 3: Orbit of the x'axis of the ellipsoid with kb
0.14, kc 0.10 in figure 2 (a). The Stokes
number is 1, 10, 100, 1000 from (a) to (d).
1
0
(b) 1
0.5
0
0.5
1
1
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a) 1
y 0
1 1
1 n
5 10
t/T
Jeff
Figure 4: Time history of the zcoordinate of the x'axis
as a function of time for the orbits in figure 3
(kb = 0.14, kc = 0.10). St is 1 (blue), 10
(red), 100 (green) and 1000 (magenta).
Conclusions
The effect of particle inertia on triaxial ellipsoids in
Stokes flow has been studied by Floquet stability analy
sis and simulation of particle motion. Particle inertia is
seen to introduce a drift towards rotation with the short
est axis aligned with the vorticity axis. This state is un
stable for inertiafree prolate triaxial ellipsoids (as shown
previously by Hinch & Leal and Yarin et al.). We show
that this rotation is stabilized for sufficiently strong parti
cle inertia. The resulting particle motion can be summa
rized as follows for particles that are within the unstable
region:
At low St, the drift towards an unstable state results
in a chaotic motion after an initial transient.
The length of the transient depends on particle ge
ometry, initial condition and Stokes number (it en
ters the chaotic motion faster for increasing Stokes
numbers).
At St above the critical St the rotation around the
shortest axis is stable and the particle approaches
and asymptotically reaches this rotation with the
shortest axis aligned with the vorticity axis, again
faster for high St.
For particles whose geometry are not unstable, parti
cle inertia gives a drift towards rotation with the shortest
axis aligned with the vorticity axis. This state is reached
asymptotically and faster for high St than for low.
This means that the ultimate motion of particles de
pend on St and geometry. For prolate ellipsoids that are
within the range of unstable kb for the kc in question,
a chaotic motion will occur eventually if St is below
(b) l
Y
1
1
(c) 1
Y 0
(d) 1 
Y 0
S 0
1 1
x
Figure 5: Orbit of the x'axis of the ellipsoid with kb
0.22, kc 0.10 in figure 2 (b). The Stokes
number is 1, 10, 100, 1000 from (a) to (d).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a) 1
0.5 j'
1
0 5 10
t/T
Jeff
15
0.5
0
0.5
1.
1
Figure 6: Time history of the zcoordinate of the x'axis
as a function of time for the orbits in figure 5
(kb = 0.22, k, = 0.10). St is 1 (blue), 10
(red), 100 (green) and 1000 (magenta).
the critical value. For other geometries, and higher St,
the particle will end up rotating around its shortest axis,
which is aligned with the vorticity axis.
(b) 1
0.5
Acknowledgements
This work is supported by the Swedish Research Coun
cil. Discussions with Ame Nordmark and Harry
Dankowicz regarding quaternions and Floquet analysis
have been of critical importance for this work.
0
1
z
1
x
1
References
Altenbach, H., Naumenko, K., Pylypenko, S. and Ren
ner, B. Influence of rotary inertia on the fiber dynam
ics in homogeneous creeping flows, Z. Angew. Math.
Mech., Vol. 87, pp. 8193, 2007
Ding, E.J. and Aidun, C.K. The dynamics and scaling
law for particles suspended in shear flow with inertia. J.
Fluid Mech. Vol. 423, pp. 317344, 2000
Hinch, E.J. and Leal, L.G. Rotation of smal non
axisymmetric particles in a simple shear flow. J. Fluid
Mech. Vol. 92, pp. 591608, 1979
Jeffery, G.B. The motion of ellipsoidal particles im
mersed in a viscous fluid. Proc. Roy. Soc. A, Vol. 102,
pp. 161179, 1922
Joseph, D.D. Stability of Fluid Motions, Springer, 1976
Lundell, F. and Carlsson, A. Heavy ellipsoids in creep
ing shear flow: Transitions of the particle rotation rate
and orbit shape. Phys. Rev. E. Vol. 81, 016323, 2010
(c)1
0.5
0
0.5
1 L 0 1
1 1
z X
Figure 7: Ellipsoid, shear flow and initial configuration
for the results shown in figures 8 13. kc
0.1; (a): kb 0.10, (b): kb 0.71 and (c):
kb 1. The angle between the x' and xaxes
is 15 in (a) and 75 in (b,c).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a) 1
y 0
1
1
(b)
Y o
1
1
(c) 1
Y o
1
1
(d)
Y 0
1
1
0
1 1
x
Figure 8: Orbit of the x'axis of the spheroid with kb
kc 0.1 in figure 7 (a). The Stokes number
is 1, 10, 100, 1000 from (a) to (d).
Figure 9: Time history of the zcoordinate of the x'axis
as a function of time for the orbits in figure 8
(kb =kc= 0.1). St is 1 (blue), 10 (red), 100
(green) and 1000 (magenta).
Qi, D. and Luo, LS. Rotational and orientational be
haviour of threedimensional spheroidal particles in
Couette flows. J. Fluid Mech. Vol. 477, pp. 201213,
2003
Shivarama, R. and Fahrenthold, E.P. Hamilton's equa
tions with Euler parameters for rigid body dynamics
modeling. J. Dyn. Sys. Meas. Control, Vol. 126, pp.
124130,2004
Subramanian, G. and Koch, D.L. Inertial effects on the
orientation of nearly spherical particles in simple shear
flow. J. Fluid Mech. Vol. 557, pp. 257296, 2006
Yarin, A.L., Gottlieb, O. and Roisman, I.V Chaotic ro
tation of triaxial ellipsoids in simple shear flow. J. Fluid
Mech. Vol. 340, pp. 83100, 1997
Appendix
In this appendix, motion of ellipsoids with kc = 0.1 and
kb 0.1, 0.71 and 1 at St 1, 10, 100 and 1000 are
shown for completeness. These parameter combinations
are marked with x in figure 1. The initial configura
tions are shown in figure 7. The axisymmetric prolate
spheroid (kb kc 0.1) is starting with the same ori
entation of the axes as the prolate triaxial ellipsoids re
ported in figures 26. For the oblate spheroid (kb 1,
kc = 0.1) and the oblate triaxial ellipsoid (kb 0.71,
kc 0.1), the initial orientation has the x' and z' axes
in the ./plane and an angle of 15 between the x' and
z axes. This initial configuration was chosen since this
is close to the configuration that has been shown to be
unstable for oblate triaxial ellipsoids by Hinch & Leal
(1979) and Yarin, Gottlieb & Roisman (1997), i.e. rota
tion around the longest axis.
0 5 10 15
t/T
Jeff
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The motion of the prolate spheroid is similar to that
reported by Lundell & Carlsson (2010). At low St there
is a kayaking with a drift towards rotation around the (a) 1
vorticity axis. As St is increased, the drift is faster (fig 
ure 8 ac). For very high St (figure 8 d), the approach to
rotation around the zaxis changes character as in figure
3 (d) and 5 (d). Since the oblate spheroid is axisymmet *
ric, the vairations of the kayaking motion seen for the
prolate triaxial ellipsoids is not present in this case and 1 1
the time history of z, is smoother in figure 9. 1 0
Turning to the oblate spheroid in figures 10 and 11 1 1
and the oblate triaxial ellipsoid in figures 12 and 13, the
effect of particle inertia and the instability, are less dra
matic albeit substantial. For these cases, the z'axis (the
shortest axis) is used to track the motion rather than the (b) 1
x'axis used before. Furthermore, the time histories are
given for the y coordinate of the z' axis, denoted zy. It
is seen that particle inertia induces a drift of the z'axis Y 0
towards the zaxis. The approach is kayaking for low St
in figures 10 (a,b) and 12 (a,b). At intermediate St, fig
ures 10 (c) and 12 (c) show the approach is rather direct.
Finally, there is a spiralling approach at the highest St 1
studied in figure 10 (d) and 12 (d). The transition from 0 1 1
the kayaking to the spiralling approach is similar albeit z x
not identical to the transition from kayaking to inclined
rotation for prolate spheroids and prolate ellipsoids re
ported above and by Lundell & Carlsson (2010).
(c)1
Y O
1 1
1 X
1 1
z x
(d) 1
Y o
1 1
1 0 0
1 1
z x
Figure 10: Orbit of the i'axis of the oblate spheroid
with kb = 1, c = 0.10 in figure 7 (b). The
Stokes number is 1, 10, 100, 1000 from (a)
to (d).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a) 1
yo
1
1
(b)
Y
1
1
(c) 1
Y o
5 10
tT
Jeff
Figure 11: Time history of the zcoordinate of the x'
axis as a function of time for the orbits in
figure 10 (kb = 1, kc = 0.10). St is 1 (blue),
10 (red), 100 (green) and 1000 (magenta).
(d) 1
Y 0
1
1
0
1 1
z x
Figure 12: Orbit of the x'axis of the oblate ellipsoid
with kb = 0.71, kc = 0.10 in figure 2 (c).
The Stokes number is 1, 10, 100, 1000 from
(a) to (d).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
5 10
t/T
Jeff
Figure 13: Time history of the zcoordinate of the x'
axis as a function of time for the orbits in
figure 12 ( kb = 0.71, kc = 0.10). St is 1
(blue), 10 (red), 100 (green) and 1000 (ma
genta).
0.5
1
0
