2
SINGLE PARTICLE MOTION
2.1 INTRODUCTION
This chapter will briefly review the issues and problems involved in con-
structing the equations of motion for individual particles, drops or bubbles
moving through a fluid. For convenience we shall use the generic name
par-
ticle
to refer to the finite pieces of the disperse phase or component. The
analyses are implicitly confined to those circumstances in which the interac-
tions between neighboring particles are negligible. In very dilute multiphase
flows in which the particles are very small compared with the global dimen-
sions of the flow and are very far apart compared with the particle size, it
is often sufficient to solve for the velocity and pressure,
u
i
(
x
i
,t
)and
p
(
x
i
,t
),
of the continuous suspending fluid while ignoring the particles or disperse
phase. Given this solution one could then solve an equation of motion for
the particle to determine its trajectory. This chapter will focus on the con-
struction of such a particle or bubble equation of motion.
The body of fluid mechanical literature on the subject of flows around
particles or bodies is very large indeed. Here we present a summary that
focuses on a spherical particle of radius,
R
, and employs the following com-
mon notation. The components of the translational velocity of the center
of the particle will be denoted by
V
i
(
t
). The velocity that the fluid would
have had at the location of the particle center in the absence of the particle
will be denoted by
U
i
(
t
). Note that such a concept is difficult to extend to
the case of interactive multiphase flows. Finally, the velocity of the particle
relative to the fluid is denoted by
W
i
(
t
)=
V
i
−
U
i
.
Frequently the approach used to construct equations for
V
i
(
t
)(or
W
i
(
t
))
given
U
i
(
x
i
,t
) is to individually estimate all the fluid forces acting on the
particle and to equate the total fluid force,
F
i
,to
m
p
dV
i
/dt
(where
m
p
is
the particle mass, assumed constant). These fluid forces may include forces
52
due to buoyancy, added mass, drag, etc. In the absence of fluid acceleration
(
dU
i
/dt
= 0) such an approach can be made unambiguously; however, in
the presence of fluid acceleration, this kind of heuristic approach can be
misleading. Hence we concentrate in the next few sections on a fundamental
fluid mechanical approach, that minimizes possible ambiguities. The classical
results for a spherical particle or bubble are reviewed first. The analysis is
confined to a suspending fluid that is incompressible and Newtonian so that
the basic equations to be solved are the continuity equation
∂u
j
∂x
j
=0
(2.1)
and the Navier-Stokes equations
ρ
C
∂u
i
∂t
+
u
j
∂u
i
∂x
j
=
−
∂p
∂x
i
+
ρ
C
ν
C
∂
2
u
i
∂x
j
∂x
j
(2.2)
where
ρ
C
and
ν
C
are the density and kinematic viscosity of the suspending
fluid. It is assumed that the only external force is that due to gravity,
g
.
Then the actual pressure is
p
=
p
−
ρ
C
gz
where
z
is a coordinate measured
vertically upward.
Furthermore, in order to maintain clarity we confine our attention to
rectilinear relative motion in a direction conveniently chosen to be the
x
1
direction.
2.2 FLOWS AROUND A SPHERE
2.2.1 At high Reynolds number
For steady flows about a sphere in which
dU
i
/dt
=
dV
i
/dt
=
dW
i
/dt
=0,it
is convenient to use a coordinate system,
x
i
, fixed in the particle as well as
polar coordinates (
r, θ
) and velocities
u
r
,u
θ
as defined in figure 2.1.
Then equations 2.1 and 2.2 become
1
r
2
∂
∂r
(
r
2
u
r
)+
1
r
sin
θ
∂
∂θ
(
u
θ
sin
θ
) = 0
(2.3)
and
ρ
C
∂u
r
∂t
+
u
r
∂u
r
∂r
+
u
θ
r
∂u
r
∂θ
−
u
2
θ
r
=
−
∂p
∂r
(2.4)
+
ρ
C
ν
C
1
r
2
∂
∂r
r
2
∂u
r
∂r
+
1
r
2
sin
θ
∂
∂θ
sin
θ
∂u
r
∂θ
−
2
u
r
r
2
−
2
r
2
∂u
θ
∂θ
53
Figure 2.1.
Notation for a spherical particle.
ρ
C
∂u
θ
∂t
+
u
r
∂u
θ
∂r
+
u
θ
r
∂u
θ
∂θ
+
u
r
u
θ
r
=
−
1
r
∂p
∂θ
(2.5)
+
ρ
C
ν
C
1
r
2
∂
∂r
r
2
∂u
θ
∂r
+
1
r
2
sin
θ
∂
∂θ
sin
θ
∂u
θ
∂θ
+
2
r
2
∂u
r
∂θ
−
u
θ
r
2
sin
2
θ
The Stokes streamfunction,
ψ
, is defined to satisfy continuity automatically:
u
r
=
1
r
2
sin
θ
∂ψ
∂θ
;
u
θ
=
−
1
r
sin
θ
∂ψ
∂r
(2.6)
and the inviscid potential flow solution is
ψ
=
−
Wr
2
2
sin
2
θ
−
D
r
sin
2
θ
(2.7)
u
r
=
−
W
cos
θ
−
2
D
r
3
cos
θ
(2.8)
u
θ
=+
W
sin
θ
−
D
r
3
sin
θ
(2.9)
φ
=
−
Wr
cos
θ
+
D
r
2
cos
θ
(2.10)
where, because of the boundary condition (
u
r
)
r
=
R
= 0, it follows that
D
=
−
WR
3
/
2. In potential flow one may also define a velocity potential,
φ
,such
that
u
i
=
∂φ/∂x
i
. The classic problem with such solutions is the fact that
the drag is zero, a circumstance termed D’Alembert’s paradox. The flow is
symmetric about the
x
2
x
3
plane through the origin and there is no wake.
The real viscous flows around a sphere at large Reynolds numbers,
54
Figure 2.2.
Smoke visualization of the nominally steady flows (from
left to right) past a sphere showing, at the top, laminar separation at
Re
=2
.
8
×
10
5
and, on the bottom, turbulent separation at
Re
=3
.
9
×
10
5
.
Photographs by F.N.M.Brown, reproduced with the permission of the Uni-
versity of Notre Dame.
Re
=2
WR/ν
C
>
1, are well documented. In the range from about 10
3
to
3
×
10
5
, laminar boundary layer separation occurs at
θ
∼
=
84
◦
and a large
wake is formed behind the sphere (see figure 2.2). Close to the sphere the
near-wake
is laminar; further downstream transition and turbulence occur-
ring in the shear layers spreads to generate a turbulent
far-wake.
As the
Reynolds number increases the shear layer transition moves forward until,
quite abruptly, the turbulent shear layer reattaches to the body, resulting
in a major change in the final position of separation (
θ
∼
=
120
◦
)andinthe
form of the turbulent wake (figure 2.2). Associated with this change in flow
55
Figure 2.3.
Drag coefficient on a sphere as a function of Reynolds number.
Dashed curves indicate the drag crisis regime in which the drag is very
sensitive to other factors such as the free stream turbulence.
pattern is a dramatic decrease in the drag coefficient,
C
D
(defined as the
drag force on the body in the negative
x
1
direction divided by
1
2
ρ
C
W
2
πR
2
),
from a value of about 0
.
5 in the laminar separation regime to a value of
about 0
.
2 in the turbulent separation regime (figure 2.3). At values of
Re
less than about 10
3
the flow becomes quite unsteady with periodic shedding
of vortices from the sphere.
2.2.2 At low Reynolds number
At the other end of the Reynolds number spectrum is the classic Stokes
solution for flow around a sphere. In this limit the terms on the left-hand side
of equation 2.2 are neglected and the viscous term retained. This solution
has the form
ψ
=sin
2
θ
−
Wr
2
2
+
A
r
+
Br
(2.11)
u
r
=cos
θ
−
W
+
2
A
r
3
+
2
B
r
(2.12)
u
θ
=
−
sin
θ
−
W
−
A
r
3
+
B
r
(2.13)
where
A
and
B
are constants to be determined from the boundary conditions
on the surface of the sphere. The force,
F
,onthe
particle
in the
x
1
direction
56
is
F
1
=
4
3
πR
2
ρ
C
ν
C
−
4
W
R
+
8
A
R
4
+
2
B
R
2
(2.14)
Several subcases of this solution are of interest in the present context. The
first is the classic Stokes (1851) solution for a solid sphere in which the no-slip
boundary condition, (
u
θ
)
r
=
R
= 0, is applied (in addition to the kinematic
condition (
u
r
)
r
=
R
= 0). This set of boundary conditions, referred to as the
Stokes boundary conditions, leads to
A
=
−
WR
3
4
,B
=+
3
WR
4
and
F
1
=
−
6
πρ
C
ν
C
WR
(2.15)
The second case originates with Hadamard (1911) and Rybczynski (1911)
who suggested that, in the case of a bubble, a condition of zero shear stress
on the sphere surface would be more appropriate than a condition of zero
tangential velocity,
u
θ
. Then it transpires that
A
=0
,B
=+
WR
2
and
F
1
=
−
4
πρ
C
ν
C
WR
(2.16)
Real bubbles may conform to either the Stokes or Hadamard-Rybczynski
solutions depending on the degree of contamination of the bubble surface,
as we shall discuss in more detail in section 3.3. Finally, it is of interest to
observe that the potential flow solution given in equations 2.7 to 2.10 is also
asubcasewith
A
=+
WR
3
2
,B
=0
and
F
1
= 0
(2.17)
However, another paradox, known as the Whitehead paradox, arises when
the validity of these Stokes flow solutions at small (rather than zero)
Reynolds numbers is considered. The nature of this paradox can be demon-
strated by examining the magnitude of the neglected term,
u
j
∂u
i
/∂x
j
,in
the Navier-Stokes equations relative to the magnitude of the retained term
ν
C
∂
2
u
i
/∂x
j
∂x
j
. As is evident from equation 2.11, far from the sphere the
former is proportional to
W
2
R/r
2
whereas the latter behaves like
ν
C
WR/r
3
.
It follows that although the retained term will dominate close to the body
(provided the Reynolds number
Re
=2
WR/ν
C
1), there will always be a
radial position,
r
c
,givenby
R/r
c
=
Re
beyond which the neglected term will
exceed the retained viscous term. Hence, even if
Re
1, the Stokes solution
is not uniformly valid. Recognizing this limitation, Oseen (1910) attempted
to correct the Stokes solution by retaining in the basic equation an approxi-
mation to
u
j
∂u
i
/∂x
j
that would be valid in the far field,
−
W∂u
i
/∂x
1
.Thus
57
the Navier-Stokes equations are approximated by
−
W
∂u
i
∂x
1
=
−
1
ρ
C
∂p
∂x
i
+
ν
C
∂
2
u
i
∂x
j
∂x
j
(2.18)
Oseen was able to find a closed form solution to this equation that satisfies
the Stokes boundary conditions approximately:
ψ
=
−
WR
2
r
2
sin
2
θ
2
R
2
+
R
sin
2
θ
4
r
+
3
ν
C
(1 + cos
θ
)
2
WR
1
−
e
Wr
2
ν
C
(1
−
cos
θ
)
(2.19)
which yields a drag force
F
1
=
−
6
πρ
C
ν
C
WR
1+
3
16
Re
(2.20)
It is readily shown that equation 2.19 reduces to equation 2.11 as
Re
→
0.
The corresponding solution for the Hadamard-Rybczynski boundary con-
ditions is not known to the author; its validity would be more question-
able since, unlike the case of Stokes boundary conditions, the inertial terms
u
j
∂u
i
/∂x
j
are not identically zero on the surface of the bubble.
Proudman and Pearson (1957) and Kaplun and Lagerstrom (1957) showed
that Oseen’s solution is, in fact, the first term obtained when the method
of matched asymptotic expansions is used in an attempt to patch together
consistent asymptotic solutions of the full Navier-Stokes equations for both
the near field close to the sphere and the far field. They also obtained the
next term in the expression for the drag force.
F
1
=
−
6
πρ
C
ν
C
WR
1+
3
16
Re
+
9
160
Re
2
ln
Re
2
+0(
Re
2
)
(2.21)
The additional term leads to an error of 1% at
Re
=0
.
3 and does not,
therefore, have much practical consequence.
The most notable feature of the Oseen solution is that the geometry of
the streamlines depends on the Reynolds number. The downstream flow is
not
a mirror image of the upstream flow as in the Stokes or potential flow
solutions. Indeed, closer examination of the Oseen solution reveals that,
downstream of the sphere, the streamlines are further apart and the flow
is slower than in the equivalent upstream location. Furthermore, this effect
increases with Reynolds number. These features of the Oseen solution are
entirely consistent with experimental observations and represent the initial
development of a wake behind the body.
The flow past a sphere at Reynolds numbers between about 0
.
5andseveral
thousand has proven intractable to analytical methods though numerical so-
58
lutions are numerous. Experimentally, it is found that a recirculating zone
(or vortex ring) develops close to the rear stagnation point at about
Re
=30
(see Taneda 1956 and figure 2.4). With further increase in the Reynolds
number this recirculating zone or wake expands. Defining locations on the
surface by the angle from the front stagnation point, the separation point
moves forward from about 130
◦
at
Re
= 100 to about 115
◦
at
Re
= 300. In
the process the wake reaches a diameter comparable to that of the sphere
when
Re
≈
130. At this point the flow becomes unstable and the ring vor-
tex that makes up the wake begins to oscillate (Taneda 1956). However, it
continues to be attached to the sphere until about
Re
= 500 (Torobin and
Gauvin 1959).
At Reynolds numbers above about 500, vortices begin to be shed and
then convected downstream. The frequency of vortex shedding has not been
studied as extensively as in the case of a circular cylinder and seems to vary
more with Reynolds number. In terms of the conventional Strouhal number,
Str
, defined as
Str
=2
fR/W
(2.22)
the vortex shedding frequencies,
f
, that Moller (1938) observed correspond
to a range of
Str
varying from 0
.
3at
Re
= 1000 to about 1
.
8at
Re
= 5000.
Furthermore, as
Re
increases above 500 the flow develops a fairly steady
near-wake
behind which vortex shedding forms an unsteady and increasingly
turbulent
far-wake.
This process continues until, at a value of
Re
of the order
of 1000, the flow around the sphere and in the near-wake again becomes
quite steady. A recognizable boundary layer has developed on the front of
the sphere and separation settles down to a position about 84
◦
from the
front stagnation point. Transition to turbulence occurs on the free shear
layer (which defines the boundary of the near-wake) and moves progressively
forward as the Reynolds number increases. The flow is similar to that of the
top picture in figure 2.2. Then the events described in the previous section
occur with further increase in the Reynolds number.
Since the Reynolds number range between 0
.
5 and several hundred can
often pertain in multiphase flows, one must resort to an empirical formula
for the drag force in this regime. A number of empirical results are available;
for example, Klyachko (1934) recommends
F
1
=
−
6
πρ
C
ν
C
WR
1+
Re
2
3
6
(2.23)
which fits the data fairly well up to
Re
≈
1000. At
Re
= 1 the factor in the
59
Re
=9
.
15
Re
=37
.
7
Re
=17
.
9
Re
=73
.
6
Re
=25
.
5
Re
= 118
Re
=26
.
8
Re
= 133
Figure 2.4.
Streamlines of steady flow (from left to right) past a sphere at
various Reynolds numbers (from Taneda 1956, reproduced by permission
of the author).
60
square brackets is 1
.
167, whereas the same factor in equation 2.20 is 1
.
187.
On the other hand, at
Re
= 1000, the two factors are respectively 17
.
7and
188
.
5.
2.2.3 Molecular effects
When the mean free path of the molecules in the surrounding fluid,
λ
,be-
comes comparable with the size of the particles, the flow will clearly deviate
from the continuum models, that are only relevant when
λ
R
. The Knud-
sen number,
Kn
=
λ/
2
R
, is used to characterize these circumstances, and
Cunningham (
1910) showed that the first-order correction for small but finite
Knudsen number leads to an additional factor, (1 + 2
AKn
), in the Stokes
drag for a spherical particle. The numerical factor, A, is roughly a constant
of order unity (see, for example, Green and Lane 1964).
When the impulse generated by the collision of a single fluid molecule
with the particle is large enough to cause significant change in the particle
velocity, the resulting random motions of the particle are called
Brownian
motion
(Einstein 1956). This leads to diffusion of solid particles suspended
in a fluid. Einstein showed that the diffusivity,
D
, of this process is given by
D
=
kT/
6
πμ
C
R
(2.24)
where
k
is Boltzmann’s constant. It follows that the typical
rms
displace-
ment of the particle in a time,
t
,isgivenby(
kT t/
3
πμ
C
R
)
1
2
.Brownian
motion is usually only significant for micron- and sub-micron-sized parti-
cles. The example quoted by Einstein is that of a 1
μm
diameter particle
in water at 17
◦
C
for which the typical displacement during one second is
0
.
8
μm
.
A third, related phenomenon is the response of a particle to the collisions
of molecules when there is a significant temperature gradient in the fluid.
Then the impulses imparted to the particle by molecular collisions on the
hot side of the particle will be larger than the impulses on the cold side. The
particle will therefore experience a net force driving it in the direction of the
colder fluid. This phenomenon is known as
thermophoresis
(see, for example,
Davies 1966). A similar phenomenon known as
photophoresis
occurs when
a particle is subjected to nonuniform radiation. One could include in this
list the Bjerknes forces described in the section 3.4 since they constitute
sonophoresis
, namely forces acting on a particle in a sound field.
61
2.3 UNSTEADY EFFECTS
2.3.1 Unsteady particle motions
Having reviewed the steady motion of a particle relative to a fluid, we must
now consider the consequences of unsteady relative motion in which either
the particle or the fluid or both are accelerating. The complexities of fluid
acceleration are delayed until the next section. First we shall consider the
simpler circumstance in which the fluid is either at rest or has a steady
uniform streaming motion (
U
= constant) far from the particle. Clearly the
second case is readily reduced to the first by a simple Galilean transformation
and it will be assumed that this has been accomplished.
In the ideal case of unsteady inviscid potential flow, it can then be shown
by using the concept of the total kinetic energy of the fluid that the force
on a rigid particle in an incompressible flow is given by
F
i
,where
F
i
=
−
M
ij
dV
j
dt
(2.25)
where
M
ij
is called the
added mass matrix
(or tensor) though the name
induced inertia tensor
used by Batchelor (1967) is, perhaps, more descrip-
tive. The reader is referred to Sarpkaya and Isaacson (1981), Yih (1969), or
Batchelor (1967) for detailed descriptions of such analyses. The above men-
tioned methods also show that
M
ij
for any finite particle can be obtained
from knowledge of several
steady
potential flows. In fact,
M
ij
=
ρ
C
2
volume
of f luid
u
ik
u
jk
d
(
volume
)
(2.26)
where the integration is performed over the entire volume of the fluid. The
velocity field,
u
ij
, is the fluid velocity in the
i
direction caused by the
steady
translation of the particle with unit velocity in the
j
direction. Note that
this means that
M
ij
is necessarily a symmetric matrix. Furthermore, it is
clear that particles with planes of symmetry will not experience a force
perpendicular to that plane when the direction of acceleration is parallel to
that plane. Hence if there is a plane of symmetry perpendicular to the
k
direction, then for
i
=
k
,
M
ki
=
M
ik
= 0, and the only off-diagonal matrix
elements that can be nonzero are
M
ij
,
j
=
k
,
i
=
k
. In the special case of
the sphere
all
the off-diagonal terms will be zero.
Tables of some available values of the diagonal components of
M
ij
are
given by Sarpkaya and Isaacson (1981) who also summarize the experi-
mental results, particularly for planar flows past cylinders. Other compila-
tions of added mass results can be found in Kennard (1967), Patton (1965),
62
Table 2.1. Added masses (diagonal terms in
M
ij
) for some three-dimensional
bodies (particles): (T) Potential flow calculations, (E) Experimental data
from Patton (1965).
and Brennen (1982). Some typical values for three-dimensional particles are
listed in Table 2.1. The uniform diagonal value for a sphere (often referred to
simply as the added mass of a sphere) is 2
ρ
C
πR
3
/
3 or one-half the displaced
mass of fluid. This value can readily be obtained from equation 2.26 using
the steady flow results given in equations 2.7 to 2.10. In general, of course,
there is
no
special relation between the added mass and the displaced mass.
63
Consider, for example, the case of the infinitely thin plate or disc with zero
displaced mass which has a finite added mass in the direction normal to the
surface. Finally, it should be noted that the literature contains little, if any,
information on off-diagonal components of added mass matrices.
Now consider the application of these potential flow results to real viscous
flows at high Reynolds numbers (the case of low Reynolds number flows will
be discussed in section 2.3.4). Significant doubts about the applicability of
the added masses calculated from potential flow analysis would be justified
because of the experience of D’Alembert’s paradox for steady potential flows
and the substantial difference between the streamlines of the potential and
actual flows. Furthermore, analyses of experimental results will require the
separation of the
added mass
forces from the viscous drag forces. Usually
this is accomplished by heuristic summation of the two forces so that
F
i
=
−
M
ij
dV
j
dt
−
1
2
ρ
C
AC
ij
|
V
j
|
V
j
(2.27)
where
C
ij
is a lift and drag coefficient matrix and
A
is a typical cross-
sectional area for the body. This is known as Morison’s equation (see Morison
et al.
1950).
Actual unsteady high Reynolds number flows are more complicated and
not necessarily compatible with such simple superposition. This is reflected
in the fact that the coefficients,
M
ij
and
C
ij
, appear from the experimental
results to be not only functions of
Re
but also functions of the reduced
time or frequency of the unsteady motion. Typically experiments involve
either oscillation of a body in a fluid or a
cceleration from rest. The most
extensively studied case involves planar flow past a cylinder (for example,
Keulegan and Carpenter 1958), and a detailed review of this data is included
in Sarpkaya and Isaacson (1981). For oscillatory motion of the cylinder with
velocity amplitude,
U
M
,andperiod,
t
∗
, the coefficients are functions of both
the Reynolds number,
Re
=2
U
M
R/ν
C
, and the reduced period or Keulegan-
Carpenter number,
Kc
=
U
M
t
∗
/
2
R
. When the amplitude,
U
M
t
∗
,islessthan
about 10
R
(
Kc <
5), the inertial effects dominate and
M
ii
is only a little less
than its potential flow value over a wide range of Reynolds numbers (10
4
<
Re <
10
6
). However, for larger values of
Kc
,
M
ii
can be substantially smaller
than this and, in some range of
Re
and
Kc
, may actually be negative. The
values of
C
ii
(the drag coefficient) that are deduced from experiments are
also a complicated function of
Re
and
Kc
. The behavior of the coefficients
is particularly pathological when the reduced period,
Kc
, is close to that of
vortex shedding (
Kc
of the order of 10). Large transverse or
lift
forces can be
generated under these circumstances. To the author’s knowledge, detailed
64
investigations of this kind have not been made for a spherical body, but one
might expect the same qualitative phenomena to occur.
2.3.2 Effect of concentration on added mass
Though most multiphase flow effects are delayed until later chapters it is
convenient at this point to address the issue of the effect on the added mass
of the particles in the surrounding mixture. It is to be expected that the
added mass coefficient for an individual particle would depend on the void
fraction of the surrounding medium. Zuber (
1964) first addressed this issue
using a
cell method
and found that the added mass,
M
ii
, for spherical bubbles
increased with volume fraction,
α
, like
M
ii
(
α
)
M
ii
(0)
=
(1 + 2
α
)
(1
−
α
)
=1+3
α
+
O
(
α
2
)
(2.28)
The simplistic geometry assumed in the cell method (a concentric spheri-
cal shell of fluid surrounding each spherical particle) caused later researchers
to attempt improvements to Zuber’s analysis; for example, van Wijngaarden
(1976) used an improved geometry (and the assumption of potential flow)
to study the
O
(
α
) term and found that
M
ii
(
α
)
M
ii
(0)
=1+2
.
76
α
+
O
(
α
2
)
(2.29)
which is close to Zuber’s result. However, even more accurate and more
recent analyses by Sangani
et al.
(1991) have shown that Zuber’s original
result is, in fact, remarkably accurate even up to volume fractions as large
as 50% (see also Zhang and Prosperetti 1994).
2.3.3 Unsteady potential flow
In general, a particle moving in any flow other than a steady uniform stream
will experience fluid accelerations, and it is therefore necessary to consider
the structure of the equation governing the particle motion under these
circumstances. Of course, this will include the special case of acceleration of
a particle in a fluid at rest (or with a steady streaming motion). As in the
earlier sections we shall confine the detailed solutions to those for a spherical
particle or bubble. Furthermore, we consider only those circumstances in
which both the particle and fluid acceleration are in one direction, chosen
for convenience to be the
x
1
direction. The effect of an external force field
65
such as gravity will be omitted; it can readily be inserted into any of the
solutions that follow by the addition of the conventional buoyancy force.
All the solutions discussed are obtained in an accelerating frame of refer-
ence
fixed
in the center of the fluid particle. Therefore, if the velocity of the
particle in some original, noninertial coordinate system,
x
∗
i
,was
V
(
t
)inthe
x
∗
1
direction, the Navier-Stokes equations in the new frame,
x
i
,fixedinthe
particle center are
∂u
i
∂t
+
u
j
∂u
i
∂x
j
=
−
1
ρ
C
∂P
∂x
i
+
ν
C
∂
2
u
i
∂x
j
∂x
j
(2.30)
where the pseudo-pressure,
P
, is related to the actual pressure,
p
,by
P
=
p
+
ρ
C
x
1
dV
dt
(2.31)
Here the conventional time derivative of
V
(
t
) is denoted by
d/dt
, but it
should be noted that in the original
x
∗
i
frame it implies a Lagrangian deriva-
tive following the particle. As before, the fluid is assumed incompressible
(so that continuity requires
∂u
i
/∂x
i
= 0) and Newtonian. The velocity that
the fluid would have at the
x
i
origin in the absence of the particle is then
W
(
t
)inthe
x
1
direction. It is also convenient to define the quantities
r
,
θ
,
u
r
,
u
θ
as shown in figure 2.1 and the Stokes streamfunction as in equations
2.6. In some cases we shall also be able to consider the unsteady effects due
to growth of the bubble so the radius is denoted by
R
(
t
).
First consider inviscid potential flow for which equations 2.30 may be
integrated to obtain the Bernoulli equation
∂φ
∂t
+
P
ρ
C
+
1
2
(
u
2
θ
+
u
2
r
)=
constant
(2.32)
where
φ
is a velocity potential (
u
i
=
∂φ/∂x
i
)and
ψ
must satisfy the equation
Lψ
=0 where
L
≡
∂
2
∂r
2
+
sin
θ
r
2
∂
∂θ
1
sin
θ
∂
∂θ
(2.33)
This is of course the same equation as in steady flow and has harmonic
solutions, only five of which are necessary for present purposes:
ψ
=sin
2
θ
−
Wr
2
2
+
D
r
+cos
θ
sin
2
θ
2
Ar
3
3
−
B
r
2
+
E
cos
θ
(2.34)
φ
=cos
θ
−
Wr
+
D
r
2
+(cos
2
θ
−
1
3
)
Ar
2
+
B
r
3
+
E
r
(2.35)
66
u
r
=cos
θ
−
W
−
2
D
r
3
+(cos
2
θ
−
1
3
)
2
Ar
−
3
B
r
4
−
E
r
2
(2.36)
u
θ
=
−
sin
θ
−
W
+
D
r
3
−
2cos
θ
sin
θ
Ar
+
B
r
4
(2.37)
The first part, which involves
W
and
D
, is identical to that for steady
translation. The second, involving
A
and
B
, will provide the fluid velocity
gradient in the
x
1
direction, and the third, involving
E
, permits a time-
dependent particle (bubble) radius. The
W
and
A
terms represent the fluid
flow in the absence of the particle, and the
D, B,
and
E
terms allow the
boundary condition
(
u
r
)
r
=
R
=
dR
dt
(2.38)
to be satisfied provided
D
=
−
WR
3
2
,B
=
2
AR
5
3
,E
=
−
R
2
dR
dt
(2.39)
In the absence of the particle the velocity of the fluid at the origin,
r
=0,is
simply
−
W
in the
x
1
direction and the gradient of the velocity
∂u
1
/∂x
1
=
4
A/
3. Hence
A
is determined from the fluid velocity gradient in the original
frame as
A
=
3
4
∂U
∂x
∗
1
(2.40)
Now the force,
F
1
, on the bubble in the
x
1
direction is given by
F
1
=
−
2
πR
2
π
0
p
sin
θ
cos
θdθ
(2.41)
which upon using equations 2.31, 2.32, and 2.35 to 2.37 can be integrated
to yield
F
1
2
πR
2
ρ
C
=
−
D
Dt
(
WR
)
−
4
3
RW A
+
2
3
R
dV
dt
(2.42)
Reverting to the original coordinate system and using
v
as the sphere volume
for convenience (
v
=4
πR
3
/
3), one obtains
F
1
=
−
1
2
ρ
C
v
dV
dt
∗
+
3
2
ρ
C
v
DU
Dt
∗
+
1
2
ρ
C
(
U
−
V
)
dv
dt
∗
(2.43)
67
where the two Lagrangian time derivatives are defined by
D
Dt
∗
≡
∂
∂t
∗
+
U
∂
∂x
∗
1
(2.44)
d
dt
∗
≡
∂
∂t
∗
+
V
∂
∂x
∗
1
(2.45)
Equation 2.43 is an important result, and care must be taken not to confuse
the different time derivatives contained in it. Note that in the absence of
bubble growth, of viscous drag, and of body forces, the equation of motion
that results from setting
F
1
=
m
p
dV /dt
∗
is
1+
2
m
p
ρ
C
v
dV
dt
∗
=3
DU
Dt
∗
(2.46)
where
m
p
is the mass of the
particle.
Thus for a massless bubble the accel-
eration of the bubble is three times the fluid acceleration.
In a more comprehensive study of unsteady potential flows Symington
(1978) has shown that the result for more general (i.e., noncolinear) a
ccel-
erations of the fluid and particle is merely the vector equivalent of equation
2.43:
F
i
=
−
1
2
ρ
C
v
dV
i
dt
∗
+
3
2
ρ
C
v
DU
i
Dt
∗
+
1
2
ρ
C
(
U
i
−
V
i
)
dv
dt
∗
(2.47)
where
d
dt
∗
=
∂
∂t
∗
+
V
j
∂
∂x
∗
j
;
D
Dt
∗
=
∂
∂t
∗
+
U
j
∂
∂x
∗
j
(2.48)
The first term in equation 2.47 represents the conventional added mass effect
due to the particle acceleration. The factor 3
/
2 in the second term due to
the fluid acceleration may initially seem surprising. However, it is made up
of two components:
1.
1
2
ρ
C
dV
i
/dt
∗
, which is the added mass effect of the fluid acceleration
2.
ρ
C
vDU
i
/Dt
∗
,whichisa
buoyancy
-like force due to the pressure gradient asso-
ciated with the fluid acceleration.
The last term in equation 2.47 is caused by particle (bubble) volumetric
growth,
dv/dt
∗
, and is similar in form to the force on a source in a uniform
stream.
Now it is necessary to ask how this force given by equation 2.47 should
be used in the practical construction of an equation of motion for a particle.
Frequently, a viscous drag force
F
D
i
, is quite arbitrarily added to
F
i
to
68
obtain some total
effective
force on the particle. Drag forces,
F
D
i
,withthe
conventional forms
F
D
i
=
C
D
2
ρ
C
|
U
i
−
V
i
|
(
U
i
−
V
i
)
πR
2
(
Re
1)
(2.49)
F
D
i
=6
πμ
C
(
U
i
−
V
i
)
R
(
Re
1)
(2.50)
have both been employed in the literature. It is, however, important to
recognize that there is no fundamental analytical justification for such su-
perposition of these forces. At high Reynolds numbers, we noted in the
last section that experimentally observed added masses are indeed quite
close to those predicted by potential flow within certain parametric regimes,
and hence the superposition has some experimental justification. At low
Reynolds numbers, it is improper to use the results of the potential flow
analysis. The appropriate analysis under these circumstances is examined in
the next section.
2.3.4 Unsteady Stokes flow
In order to elucidate some of the issues raised in the last section, it is instruc-
tive to examine solutions for the unsteady flow past a sphere in low Reynolds
number Stokes flow. In the asymptotic case of zero Reynolds number, the so-
lution of section 2.2.2 is unchanged by unsteadiness, and hence the solution
at any instant in time is identical to the steady-flow solution for the same
particle velocity. In other words, since the fluid has no inertia, it is always
in static equilibrium. Thus the instantaneous force is identical to that for
the steady flow with the same
V
i
(
t
).
The next step is therefore to investigate the effects of small but nonzero
inertial contributions. The Oseen solution provides some indication of the
effect of the
convective
inertial terms,
u
j
∂u
i
/∂x
j
,insteadyflow.Herewe
investigate the effects of the
unsteady
inertial term,
∂u
i
/∂t
. Ideally it would
be best to include
both
the
∂u
i
/∂t
term
and
the Oseen approximation to
the convective term,
U∂u
i
/∂x
. However, the resulting
unsteady
Oseen flow
is sufficiently difficult that only small-time expansions for the impulsively
started motions of droplets and bubbles exist in the literature (Pearcey and
Hill 1956).
Consider, therefore the
unsteady
Stokes equations in the absence of the
convective inertial terms:
ρ
C
∂u
i
∂t
=
−
∂P
∂x
i
+
μ
C
∂
2
u
i
∂x
j
∂x
j
(2.51)
69
Since both the equations and the boundary conditions used below are linear
in
u
i
, we need only consider colinear particle and fluid velocities in one
direction, say
x
1
. The solution to the general case of noncolinear particle and
fluid velocities and accelerations may then be obtained by superposition. As
in section 2.3.3 the colinear problem is solved by first transforming to an
accelerating coordinate frame,
x
i
, fixed in the center of the particle so that
P
=
p
+
ρ
C
x
1
dV /dt
. Elimination of
P
by taking the curl of equation 2.51
leads to
(
L
−
1
ν
C
∂
∂t
)
Lψ
= 0
(2.52)
where
L
is the same operator as defined in equation 2.33. Guided by both
the steady Stokes flow and the unsteady potential flow solution, one can
anticipate a solution of the form
ψ
=sin
2
θf
(
r, t
)+cos
θ
sin
2
θg
(
r, t
)+cos
θh
(
t
)
(2.53)
plus other spherical harmonic functions. The first term has the form of the
steady Stokes flow solution; the last term would be required if the parti-
cle were a growing spherical bubble. After substituting equation 2.53 into
equation 2.52, the equations for
f, g, h
are
(
L
1
−
1
ν
C
∂
∂t
)
L
1
f
=0
where
L
1
≡
∂
2
∂r
2
−
2
r
2
(2.54)
(
L
2
−
1
ν
C
∂
∂t
)
L
2
g
=0
where
L
2
≡
∂
2
∂r
2
−
6
r
2
(2.55)
(
L
0
−
1
ν
C
∂
∂t
)
L
0
h
=0
where
L
0
≡
∂
2
∂r
2
(2.56)
Moreover, the form of the expression for the force,
F
1
, on the spherical
particle (or bubble) obtained by evaluating the stresses on the surface and
integrating is
F
1
4
3
ρ
C
πR
3
=
dV
dt
+
1
r
∂
2
f
∂r∂t
+
ν
C
r
2
r
2
∂f
∂r
+
2
r
∂
2
f
∂r
2
−
∂
3
f
∂r
3
r
=
R
(2.57)
It transpires that this is
independent
of
g
or
h
. Hence only the solution
to equation 2.54 for
f
(
r, t
) need be sought in order to find the force on a
spherical particle, and the other spherical harmonics that might have been
included in equation 2.53 are now seen to be unnecessary.
Fourier or Laplace transform methods may be used to solve equation 2.54
for
f
(
r, t
), and we choose Laplace transforms. The Laplace transforms for
70
the relative velocity
W
(
t
), and the function
f
(
r, t
) are denoted by
ˆ
W
(
s
)and
ˆ
f
(
r, s
):
ˆ
W
(
s
)=
∞
0
e
−
st
W
(
t
)
dt
;
ˆ
f
(
r, s
)=
∞
0
e
−
st
f
(
r, t
)
dt
(2.58)
Then equation 2.54 becomes
(
L
1
−
ξ
2
)
L
1
ˆ
f
= 0
(2.59)
where
ξ
=(
s/ν
C
)
1
2
, and the solution after application of the condition that
ˆ
u
1
(
s, t
) far from the particle be equal to
ˆ
W
(
s
)is
ˆ
f
=
−
ˆ
Wr
2
2
+
A
(
s
)
r
+
B
(
s
)(
1
r
+
ξ
)
e
−
ξr
(2.60)
where
A
and
B
are functions of
s
whose determination requires application
of the boundary conditions on
r
=
R
.Intermsof
A
and
B
the Laplace
transform of the force
ˆ
F
1
(
s
)is
ˆ
F
1
4
3
ρ
C
πR
3
=
ˆ
dV
dt
+
s
r
∂
ˆ
f
∂r
+
ν
C
R
−
4
ˆ
W
r
+
8
A
r
4
+
CBe
−
ξr
r
=
R
(2.61)
where
C
=
ξ
4
+
3
ξ
3
r
+
3
ξ
2
r
2
+
8
ξ
r
3
+
8
r
4
(2.62)
The classical solution (see Landau and Lifshitz 1959) is for a solid sphere
(i.e., constant R) using the no-slip (Stokes) boundary condition for which
f
(
R, t
)=
∂f
∂r
r
=
R
= 0
(2.63)
and hence
A
=+
ˆ
WR
3
2
+
3
ˆ
WRν
C
2
s
{
1+
ξR
}
;
B
=
−
3
ˆ
WRν
C
2
s
e
ξR
(2.64)
so that
ˆ
F
1
4
3
ρ
C
πR
3
=
ˆ
dV
dt
−
3
2
s
ˆ
W
−
9
ν
C
ˆ
W
2
R
2
−
9
ν
1
2
C
2
R
s
1
2
ˆ
W
(2.65)
For a motion starting at rest at
t
=0theinverseLaplacetransformofthis
71
yields
F
1
4
3
ρ
C
πR
3
=
dV
dt
−
3
2
dW
dt
−
9
ν
C
2
R
2
W
−
9
2
R
(
ν
C
π
)
1
2
t
0
dW
(
̃
t
)
d
̃
t
d
̃
t
(
t
−
̃
t
)
1
2
(2.66)
where
̃
t
is a dummy time variable. This result must then be written in the
original coordinate framework with
W
=
V
−
U
and can be generalized to
the noncolinear case by superposition so that
F
i
=
−
1
2
vρ
C
dV
i
dt
∗
+
3
2
vρ
C
dU
i
dt
∗
+
9
vμ
C
2
R
2
(
U
i
−
V
i
)
+
9
vρ
C
2
R
(
ν
C
π
)
1
2
t
∗
0
d
(
U
i
−
V
i
)
d
̃
t
d
̃
t
(
t
∗
−
̃
t
)
1
2
(2.67)
where
d/dt
∗
is the Lagrangian time derivative following the particle. This
is then the general force on the particle or bubble in unsteady Stokes flow
when the Stokes boundary conditions are applied.
Compare this result with that obtained from the potential flow analysis,
equation 2.47 with
v
taken as constant. It is striking to observe that the coef-
ficients of the added mass terms involving
dV
i
/dt
∗
and
dU
i
/dt
∗
are identical
to those of the potential flow solution. On superficial examination it might
be noted that
dU
i
/dt
∗
appears in equation 2.67 whereas
DU
i
/Dt
∗
appears
in equation 2.47; the difference is, however, of order
W
j
∂U
i
/dx
j
and terms
of this order have already been dropped from the equation of motion on the
basis that they were negligible compared with the temporal derivatives like
∂W
i
/∂t
. Hence it is inconsistent with the initial assumption to distinguish
between
d/dt
∗
and
D/Dt
∗
in the present unsteady Stokes flow solution.
The term 9
ν
C
W/
2
R
2
in equation 2.67 is, of course, the steady Stokes drag.
The new phenomenon introduced by this analysis is contained in the last
term of equation 2.67. This is a fading memory term that is often named the
Basset term after one of its identifiers (Basset 1888). It results from the fact
that additional vorticity created at the solid particle surface due to relative
acceleration diffuses into the flow and creates a temporary perturbation in
the flow field. Like all diffusive effects it produces an
ω
1
2
term in the equation
for oscillatory motion.
Before we conclude this section, comment should be included on three
other analytical results. Morrison and Stewart (1976) considered the case of
a spherical bubble for which the Hadamard-Rybczynski boundary conditions
rather than the Stokes conditions are applied. Then, instead of the conditions
of equation 2.63, the conditions for zero normal velocity and zero shear stress
72
on the surface require that
f
(
R, t
)=
∂
2
f
∂r
2
−
2
r
∂f
∂r
r
=
R
= 0
(2.68)
and hence in this case (see Morrison and Stewart 1976)
A
(
s
)=+
ˆ
WR
3
2
+
3
ˆ
WR
(1 +
ξR
)
ξ
2
(3 +
ξR
)
;
B
(
s
)=
−
3
ˆ
WRe
+
ξR
ξ
2
(3 +
ξR
)
(2.69)
so that
ˆ
F
1
4
3
πρ
C
R
3
=
ˆ
dV
dt
−
9
ˆ
Wν
C
R
2
−
3
2
ˆ
Ws
+
6
ν
C
ˆ
W
R
2
1+
s
1
2
R/
3
ν
1
2
C
(2.70)
The inverse Laplace transform of this for motion starting at rest at
t
=0is
F
1
4
3
ρ
C
πR
3
=
dV
dt
−
3
2
dW
dt
−
3
ν
C
W
R
2
(2.71)
−
6
ν
C
R
2
t
0
dW
(
̃
t
)
d
̃
t
exp
9
ν
C
(
t
−
̃
t
)
R
2
erfc
9
ν
C
(
t
−
̃
t
)
R
2
1
2
d
̃
t
Comparing this with the solution for the Stokes conditions, we note that the
first two terms are unchanged and the third term is the expected Hadamard-
Rybczynski steady drag term (see equation 2.16). The last term is signifi-
cantly different from the Basset term in equation 2.67 but still represents a
fading memory.
More recently, Magnaudet and Legendre (
1998) have extended these re-
sults further by obtaining an expression for the force on a particle (bubble)
whose radius is changing with time.
Another interesting case is that for unsteady Oseen flow, which essentially
consists of attempting to solve the Navier-Stokes equations with the convec-
tive inertial terms approximated by
U
j
∂u
i
/∂x
j
. Pearcey and Hill (1956) have
examined the small-time behavior of droplets and bubbles started from rest
when this term is included in the equations.
2.4 PARTICLE EQUATION OF MOTION
2.4.1 Equations of motion
In a multiphase flow with a very dilute discrete phase the fluid forces dis-
cussed in sections 2.1 to 2.3.4 will determine the motion of the particles that
73
constitute that discrete phase. In this section we discuss the implications of
some of the fluid force terms. The equation that determines the particle
velocity,
V
i
, is generated by equating the total force,
F
T
i
, on the particle
to
m
p
dV
i
/dt
∗
. Consider the motion of a spherical particle (or bubble) of
mass
m
p
and volume
v
(radius
R
)ina
uniformly
accelerating fluid. The
simplest example of this is the vertical motion of a particle under gravity,
g
, in a pool of otherwise quiescent fluid. Thus the results will be written
in terms of the buoyancy force. However, the same results apply to mo-
tion generated by any uniform acceleration of the fluid, and hence
g
can be
interpreted as a general uniform fluid acceleration (
dU/dt
). This will also
allow some tentative conclusions to be drawn concerning the relative mo-
tion of a particle in the nonuniformly accelerating fluid situations that can
occur in general multiphase flow. For the motion of a sphere at small rela-
tive Reynolds number,
Re
1(where
Re
=2
WR/ν
C
and
W
is the typical
magnitude of the relative velocity), only the forces due to buoyancy and the
weight of the particle need be added to
F
i
as given by equations 2.67 or 2.71
in order to obtain
F
T
i
. This addition is simply given by (
ρ
C
v
−
m
p
)
g
i
where
g
is a vector in the vertically upward direction with magnitude equal to the
acceleration due to gravity. On the other hand, at high relative Reynolds
numbers,
Re
1, one must resort to a more heuristic approach in which
the fluid forces given by equation 2.47 are supplemented by drag (and lift)
forces given by
1
2
ρ
C
AC
ij
|
W
j
|
W
j
as in equation 2.27. In either case it is useful
to nondimensionalize the resulting equation of motion so that the pertinent
nondimensional parameters can be identified.
Examine first the case in which the relative velocity,
W
(defined as positive
in the direction of the acceleration,
g
, and therefore positive in the vertically
upward direction of the rising bubble or sedimenting particle), is sufficiently
small so that the relative Reynolds number is much less than unity. Then,
using the Stokes boundary conditions, the equation governing
W
may be
obtained from equation 2.66 as
w
+
dw
dt
∗
+
9
π
(1 + 2
m
p
/ρ
C
v
)
1
2
t
∗
0
dw
d
̃
t
d
̃
t
(
t
∗
−
̃
t
)
1
2
= 1
(2.72)
where the dimensionless time,
t
∗
=
t/t
u
and the relaxation time,
t
u
,isgiven
by
t
u
=
R
2
(1 + 2
m
p
/ρ
C
v
)
/
9
ν
C
(2.73)
74