05 book
2007/5/15
page 1
✐
✐
✐
✐
✐
✐
✐
✐
Chapter 1
Fundamentals
When the dimensions of a body are insignificant to the description of its motion or the action
of forces on it, the body may be idealized as a particle, i.e., a piece of material occupying
a point in space and perhaps moving as time passes. In the next few sections, we briefly
review some essential concepts that will be needed later in the analysis of particles.
1.1 Notation
In this work, boldface symbols imply vectors or tensors. A fixed Cartesian coordinate
system will be used throughout. The unit vectors for such a system are given by the mutually
orthogonal triad (e
1
, e
2
, e
3
). For the inner product of two vectors u and v, we have in three
dimensions
u · v =
3


(1.2)
represents the Euclidean norm in R
3
and θ is the angle between the two vectors. We recall
that a norm has three main characteristics for any two bounded vectors u and v (||u|| < ∞
and ||v|| < ∞):
• ||u|| > 0, and ||u|| = 0 if and only if u = 0,
• ||u + v||≤||u||+||v||, and
• ||γ u||≤|γ |||u||, where γ is a scalar.
Two vectors are said to be orthogonal if u ·v = 0. The cross (vector) product of two vectors
is
u × v =−v × u =






e
1
e
2
e
3
u
1
u
2
u
3

u(t) =
du
1
(t)
dt
e
1
+
du
2
(t)
dt
e
2
+
du
3
(t)
dt
e
3
=˙u
1
e
1
+˙u
2
e
2
+˙u

. The divergence of a vector (a contraction to a scalar) is defined by
∇·u =

e
1
∂
∂x
1
+ e
2
∂
∂x
2
+ e
3
∂
∂x
3

·
(
u
1
e
1
+ u
2
e
2
+ u


e
1
e
2
e
3
∂
∂x
1
∂
∂x
2
∂
∂x
3
u
1
u
2
u
3






. (1.7)
1.2 Kinematics of a single particle

+ r
2
e
2
+ r
3
e
3
, (1.10)
v =
˙
r =˙r
1
e
1
+˙r
2
e
2
+˙r
3
e
3
, (1.11)
and
a =
¨
r =¨r
1
e

, and a
i−j
= a
i
− a
j
.
05 book
2007/5/15
page 3
✐
✐
✐
✐
✐
✐
✐
✐
1.3. Kinetics of a single particle 3
1.3 Kinetics of a single particle
Throughout this monograph, the fundamental relation between force and acceleration is
given by Newton’s second law of motion, in vector form:
 = ma, (1.13)
where  is the sum (resultant) of all the applied forces acting on mass m.
1.3.1 Work, energy, and power
A closely related concept is that of work and energy. The differential amount of work done
by a force acting through a differential displacement is
dW =  · dr. (1.14)
Therefore, the total amount of work performed by a force over a displacement history is
W

·v
2
−v
1
·v
1
)
def
= T
2
−T
1
,
(1.15)
where T
def
=
1
2
mv ·v is known as the kinetic energy.
7
Therefore, we may write
T
1
+ W
1→2
= T
2
. (1.16)
If the forces can be written in the form

)
dV = V(r(t
1
)) − V(r(t
2
)) =

r(t
1
)
r(t
2
)
dV ⇒

r(t
2
)
r(t
1
)
dV +

r(t
1
)
r(t
2
)
dV = 0. (1.21)

4 Chapter 1. Fundamentals
1.3.2 Properties of a potential
As we have indicated, a force field  is said to be conservative if and only if there exists a
continuously differentiable scalar field V such that  =−∇V . Therefore, a necessary and
sufficient condition for a particle to be in equilibrium is that
 =−∇V = 0. (1.24)
In other words,
∂V
∂x
1
= 0,
∂V
∂x
2
= 0, and
∂V
∂x
3
= 0. (1.25)
Forces acting on a particle (1) that are always directed toward or away from another point
and (2) whose magnitude depends only on the distance between the particle and the point
of attraction/repulsion are called central forces. They have the form
 =−C(||r −r
o
||)
r −r
o
||r −r
o
||

||
−β
1
+1
−β
1
+ 1

 
attraction
−
α
2
||r −r
o
||
−β
2
+1
−β
2
+ 1

 
repulsion
, (1.30)
where all of the parameters, the α’s and β’s, are nonnegative. The gradient yields
−∇V =  =

α

∂
2
V
∂x
1
∂x
1
∂
2
V
∂x
1
∂x
2
∂
2
V
∂x
1
∂x
3
∂
2
V
∂x
2
∂x
1
∂
2

3
∂x
3





, (1.32)
05 book
2007/5/15
page 5
✐
✐
✐
✐
✐
✐
✐
✐
1.3. Kinetics of a single particle 5
around an equilibrium point. A sufficient condition for V to attain a minimum at an equilib-
rium point is for the Hessian to be positive definite (which implies that V is locally convex).
For more details, see Hale and Kocak [88].
Remark. Provided that the α’s and β’s are selected appropriately, the chosen central
force potential form is stable for motion in the normal direction, i.e., the line connecting the
centers of particles in particle-particle interaction.
8
In order to determine stable parameter
combinations, consider a potential function for a single particle, in one-dimensional motion,

 =−
∂V
∂r
=

α
1
|r −r
o
|
−β
1
− α
2
|r −r
o
|
−β
2

n, (1.34)
where n =
r
o
−r
|r−r
o
|
. For stability, we require
∂

o
|) =−α
1
|r
e
−r
o
|
−β
1
+
α
2
|r
e
− r
o
|
−β
2
= 0, which implies
|r
e
− r
o
|=

α
2
α

t
1
 dt = G(t
2
), (1.38)
8
For disturbances in directions orthogonal to the normal direction, the potential is neutrally stable, i.e., the
Hessian’s determinant is zero, thus indicating that the potential does not change for such perturbations. The
motion analysis in the normal direction is relevant for central forces of the type under consideration.
05 book
2007/5/15
page 6
✐
✐
✐
✐
✐
✐
✐
✐
6 Chapter 1. Fundamentals
where
G(t
1
) = (mv)|
t=t
1
(1.39)
is the linear momentum. Clearly, if
 = 0, (1.40)

dt = H
o
(t
2
). (1.43)
Thus, if
M = 0, (1.44)
then
H
o
(t
1
) = H
o
(t
2
), (1.45)
and angular momentum is said to be conserved.
1.4 Systems of particles
We now discuss the dynamics of a system of N
p
particles. Let r
i
, i = 1, 2, 3, ,N
p
,be
the position vectors of a system of particles.
1.4.1 Linear momentum
The position vector of the center of mass of the system is given by
r

r
i
= r
cm
+ r
i−cm
. (1.47)
The linear momentum of a system of particles is given by
N
p

i=1
m
i
˙
r
i

G
i
=
N
p

i=1
m
i
(
˙
r

2007/5/15
page 7
✐
✐
✐
✐
✐
✐
✐
✐
1.4. Systems of particles 7
since
N
p

i=1
m
i
˙
r
i−cm
= 0. (1.49)
Thus, the linear momentum of any system with constant mass is the product of the mass
and the velocity of its center of mass; furthermore,
˙
G
cm
= M
¨
r

Now sum over all the particles in the system to obtain
N
p

i=1
m
i
¨
r
i
= M
¨
r
cm
=
N
p

i=1


EXT
i
+ 
INT
i

=
N
p

= M
¨
r
cm
=
N
p

i=1

EXT
i
. (1.54)
Thus, the impulse-momentum relation reads
G
cm
(t
1
) +
N
p

i=1

t
2
t
1

EXT

where

W
i,1→2
represents all of the work done by the external and internal forces. It is
advantageous to decompose the kinetic energy into the translation of the center of mass and
the motion relative to the center of mass. This is achieved by writing
v
i
= v
cm
+
˙
r
i−cm
, (1.57)
05 book
2007/5/15
page 8
✐
✐
✐
✐
✐
✐
✐
✐
8 Chapter 1. Fundamentals
which yields
N


i=1
1
2
m
i
v
cm
· v
cm
+
N
p

i=1
1
2
m
i
˙
r
i−cm
·
˙
r
i−cm
.
(1.58)
If the entire system is rigid, the second term takes on the meaning of rotation around the
center of mass.

¨
r
i
=
N
p

j=i

α
1ij
||r
i
− r
j
||
−β
1
− α
2ij
||r
i
− r
j
||
−β
2

n
ij

following dimensionless parameters:
• r
∗
def
=
r
L
,
• t
∗
def
=
t
T
.
The quantities that appear in Equation (1.59) become
• m
i
¨
r
i
= m
i
L
T
2
d
2
r
∗

2ij
||r
i
− r
j
||
−β
2
= α
2ij
L
−β
2
||r
∗
i
− r
∗
j
||
−β
2
,
where n
ij
remains unchanged. Substituting these relations into Equation (1.59) yields
d
2
r
∗

α
2ij
m
i
T
2
L
−(β
2
+1)
||r
∗
i
− r
∗
j
||
−β
2

n
ij
.
(1.61)
Thus, two dimensionless parameters, which must be the same for two systems to exhibit
similitude between one another, are
•
α
1ij
m

1
+1)

system 1
=

α
1ij
m
i
T
2
L
−(β
1
+1)

system 2
(1.62)
and

α
2ij
m
i
T
2
L
−(β
2

05 book
2007/5/15
page 11
✐
✐
✐
✐
✐
✐
✐
✐
Chapter 2
Modeling of particulate
flows
As indicated in the preface, in this introductory monograph the objects in the flow are
assumed to be small enough to be considered (idealized) as particles, spherical in shape,
and the effects of their rotation with respect to their mass center are assumed unimportant
to their overall motion.
2.1 Particulate flow in the presence of near-fields
We consider a group of nonintersecting particles (N
p
in total).
9
The equation of motion for
the ith particle in a flow is
m
i
¨
r
i

(2.2)
represents the sum of forces due to near-field interaction (
nf
), normal contact forces
(
con
), and friction (
fric
). We consider the following relatively general central-force
attraction-repulsion form for the near-field forces induced by all particles on particle i:

nf
i
=
N
p

j=i



α
1ij
||r
i
− r
j
||
−β
1

r
j
− r
i
||r
i
− r
j
||
. (2.4)
9
The approach in this chapter draws from methods developed in Zohdi [212] and [217].
11
05 book
2007/5/15
page 12
✐
✐
✐
✐
✐
✐
✐
✐
12 Chapter 2. Modeling of particulate flows
RECOVERY
COMPRESSION
CONTACT
INITIAL
Figure 2.1. Compression and recovery of two impacting particles (Zohdi [212]).

m
i
v
in
(t) +m
j
v
jn
(t) +

t+δt
t
E
i
·n
ij
dt +

t+δt
t
E
j
·n
ij
dt = m
i
v
in
(t +δt)+m
j

(t +δt), (2.6)
where

t+δt
t
I
n
dt isthe totalnormal impulse dueto impact. For a pair ofparticles undergoing
impact, let usconsider a decomposition of the collision event (Figure 2.1) into a compression
(δt
1
) and a recovery (δt
2
) phase, i.e., δt = δt
1
+δt
2
. Between the compression and recovery
10
Alternatively, if the near-fields are related to the amount of surface area, this scaling could be done per unit
area.
05 book
2007/5/15
page 13
✐
✐
✐
✐
✐
✐

ij
dt = m
i
v
cn
,
(2.7)
and, in the recovery phase,
m
i
v
cn
+

t+δt
t+δt
1
I
n
dt +

t+δt
t+δt
1
E
i
· n
ij
dt = m
i

(2.9)
and, in the recovery phase,
m
j
v
cn
−

t+δt
t+δt
1
I
n
dt +

t+δt
t+δt
1
E
j
· n
ij
dt = m
j
v
jn
(t +δt).
(2.10)
This leads to an expression for the coefficient of restitution:
e

(v
cn
− v
in
(t)) − E
in
(t, t +δt
1
)
=
−m
j
(v
jn
(t +δt) − v
cn
) + E
jn
(t +δt
1
,t + δt)
−m
j
(v
cn
− v
jn
(t)) + E
jn
(t, t +δt

t+δt
1
E
j
· n
ij
dt,
E
in
(t, t +δt
1
)
def
=

t+δt
1
t
E
i
· n
ij
dt,
E
jn
(t, t +δt
1
)
def
=

1
)
,
(2.13)
11
A common normal velocity for particles should be interpreted as indicating that the relative velocity in the
normal direction between particle centers is zero.
05 book
2007/5/15
page 14
✐
✐
✐
✐
✐
✐
✐
✐
14 Chapter 2. Modeling of particulate flows
where
12

ij
(t +δt
1
,t + δt)
def
=
1
m

1
m
j
E
jn
(t, t +δt
1
).
(2.15)
Thus, we may rewrite Equation (2.13) as
v
jn
(t +δt) = v
in
(t +δt) − 
ij
(t +δt
1
,t + δt) + e

v
in
(t) − v
jn
(t) + 
ij
(t, t +δt
1
)


in
(t) − v
jn
(t)))
m
i
+ m
j
+
(
E
in
+ E
jn
)δt −m
j
(e
ij
(t, t +δt
1
) − 
ij
(t +δt
1
,t + δt))
m
i
+ m
j
,

E
in
. (2.19)
In particular, as will be done later in the analysis, when we discretize the equations of
motion with a discrete (finite difference) time step of t, where δt  t, we shall define
the impulsive normal contact contribution to the total force acting on a particle, 
tot
i
=

nf
i
+ 
con
i
+ 
fric
i
(Equation (2.2)), to be

con
=
I
n
δt
t
n
ij
. (2.20)
Furthermore, at the implementation level, we choose δt = γt, where 0 <γ  1 and t

jn
(t +δt) −v
in
(t +δt)
v
in
(t) − v
jn
(t)
.
13
A typical choice is 0 <γ ≤ 0.01. Typically, the system is insensitive to γ below 0.01.
05 book
2007/5/15
page 15
✐
✐
✐
✐
✐
✐
✐
✐
2.3. Kinetic energy dissipation 15
V(0)V(0)
t
n
Figure 2.2. Two identical particles approaching one another (Zohdi [212]).
These results are consistent with the fact that the recovery time vanishes (all compression
and no recovery) for e → 0 (asymptotically “plastic”) and, as e → 1, the recovery time

M
N
p

i=1
m
i
v
i
(2.23)
and M =

N
p
i=1
m
i
, leading to
v
i
(t) = v
cm
(t) + δv
i
(t), (2.24)
05 book
2007/5/15
page 16
✐
✐

(v
cm
(t) + δv
i
(t)) · (v
cm
(t) + δv
i
(t))
= Mv
cm
(t) · v
cm
(t) + 2v
cm
(t) ·
N
p

i=1
m
i
δv
i
(t)

 
=0
+
N

p

i=1
m
i
δv
i
(t +δt) · δv
i
(t +δt).
(2.26)
Subtracting Equation (2.25) from Equation (2.26) yields
N
p

i=1
m
i
v
i
(t +δt) · v
i
(t +δt) −
N
p

i=1
m
i
v

p

i=1
m
i
δv
i
(t) · δv
i
(t) −
N
p

i=1
m
i
δv
i
(t) · δv
i
(t)
= (e
2
− 1)
N
p

i=1
m
i

Remark. In order to help characterize the overall behavior of the motion, it is advan-
tageous to decompose the kinetic energy per unit mass into the bulk motion of the center of
mass and the motion relative to the center of mass:
T(t)=
T(t)
M
=
1
2
v
cm
(t) · v
cm
(t)

 
def
=T
b
= bulk motion energy
+
1
2M
N
p

i=1
m
i
δv

Remark. Sometimes expressions of the form
N
p

i=1
m
i
v
i
· v
i
− Mv
cm
· v
cm
=
N
p

i=1
m
i
δv
i
· δv
i
(2.30)
are termed “granular gas temperatures.”
2.4 Incorporating friction
To incorporate frictional stick-slip phenomena during impact, for a general particle pair (i

δt + E
it
δt = m
i
v
ct
, (2.32)
where the friction contribution is
I
f
=
1
δt

t+δt
t
I
f
dt, (2.33)
the total contribution from all other particles in the tangential direction (τ
ij
)is
E
it
=
1
δt

t+δt
t

f
, which can be
solved for as
I
f
=

E
it
m
i
−
E
jt
m
j

δt + v
it
(t) − v
jt
(t)

1
m
i
+
1
m
j

s
|I
n
|, (2.37)
where
µ
s
≥ µ
d
(2.38)
is the coefficient of static friction, then slip must occur and a dynamic sliding friction model
is used. If sliding occurs, the friction force is assumed to be proportional to the normal force
and opposite to the direction of relative tangential motion, i.e.,

fric
i
def
= µ
d
||
con
||
v
jt
− v
it
||v
jt
− v
it

τ
(t))
(2.41)
and
T(t + δt) =
1
2
m(v
2
n
(t +δt) + v
2
τ
(t +δt)).
(2.42)
Assuming sliding takes place, for either particle, the impulse-momentum relation can be
written as
mv
n
(t) +

t+δt
t
I
n
dt = mv
n
(t +δt) (2.43)
05 book
2007/5/15

dt = m(v
n
(t +δt) − v
n
(t)) =−(1 + e)mv
n
(t). (2.45)
Substituting this relation into the conservation of momentum relation in the tangential di-
rection, we have
v
τ
(t +δt) = v
τ
(t) − (1 + e)v
n
(t)µ
d
. (2.46)
Now consider the restriction that the friction forces cannot be so large that they reverse the
initial tangential motion. Mathematically, this restriction can be written as
v
τ
(t +δt) = v
τ
(t) − (1 + e)v
n
(t)µ
d
≥ 0, (2.47)
which leads to the expression

v
n
(t)(1 + e)
, (2.49)
which is the maximum value of µ
d
dictated by Equation (2.48).
16
2.4.2 Velocity-dependent coefficients of restitution
It is important to realize that, in reality, the phenomenological parameter e depends on the
severity of the impact velocity. For extensive experimental data, see Goldsmith [79], or
see Johnson [111] for a more detailed analytical treatment. Qualitatively, the coefficient of
restitution has behavior as shown in Figure 2.4. A mathematical idealization of the behavior
can be constructed as
e
def
= max

e
o

1 −
v
n
v
∗

,e
−
