05 book
2007/5/15
page 20
✐
✐
✐
✐
✐
✐
✐
✐
20 Chapter 2. Modeling of particulate flows
IMPACT VELOCITY
e
−
EMPIRICALLY
OBSERVED
e
e
o
IDEALIZATION
V*
Figure 2.4. Qualitative behavior of the coefficient of restitution with impact ve-
locity (Zohdi [212]).
where v
∗
is a critical threshold velocity (normalization) parameter, the relative velocity of
approach is defined by
v
n
def

≈ G(r, t). (3.2)
A primary question is, at which time should we evaluate the equation? If we use time = t,
then
˙r|
t
=
r(t +t) −r(t)
t
= G(r(t), t) ⇒ r(t +t) = r(t) + tG(r(t), t), (3.3)
which yields an explicit expression for r(t + t). This is often referred to as a forward
Euler scheme. If we use time = t + t , then
˙r|
t+t
=
r(t +t) −r(t)
t
= G(r(t + t), t + t), (3.4)
and therefore
r(t +t) = r(t) +tG(r(t + t ), t + t), (3.5)
which yields an implicit expression, which can be nonlinear in r(t +t), depending on G.
This is often referred to as a backward Euler scheme. These two techniques illustrate the
most basic time-stepping schemes used in the scientific community, which form the founda-
tion for the majority of more sophisticated methods. Two main observations can be made:
• The implicit method usually requires one to solve a (nonlinear) equation in r(t +t).
• The explicit method has the major drawback that the step size t may have to be very
small to achieve acceptable numerical results. Therefore, an explicit simulation will
usually require many more time steps than an implicit simulation.
21
05 book
2007/5/15

o
(1 − ct)
L
, (3.8)
where L indicatesthe time step counter, t = Lt for uniformtime steps (as inthis example),
and r
L
def
= r(t), etc. It is stable if |1 − ct| < 1. For the implicit method,
˙r ≈
r(t +t) −r(t)
t
=−cr(t +t), (3.9)
which leads to the time-stepping scheme
r(Lt) =
r
o
(1 + ct)
L
. (3.10)
Since
1
1+ct
< 1, it is always stable. Note that the approximation in Equation (3.8) oscillates
in an artificial, nonphysical manner when
t >
2
c
. (3.11)
If c  1, then Equation (3.6) is a so-called stiff equation, and t may have to be very small

✐
✐
✐
✐
✐
✐
3.3. Application to particulate flows 23
one arrives at the following abstract form, for the entire system of particles:
A(r
L+1
) = F. (3.13)
It is convenient to write
A(r
L+1
) − F = G(r
L+1
) − r
L+1
+ R = 0, (3.14)
where R is a remainder term that does not depend on the solution, i.e.,
R = R(r
L+1
). (3.15)
A straightforward iterative scheme can be written as
r
L+1,K
= G(r
L+1,K−1
) + R, (3.16)
where K = 1, 2, 3, is the index of iteration within time step L +1. The convergence of

L+1,K−1
) − G(r
L+1
)|| ≤ η
L+1,K
||r
L+1,K−1
− r
L+1
||,
(3.19)
where, if
0 ≤ η
L+1,K
< 1 (3.20)
for each iteration K, then

L+1,K
→ 0 (3.21)
for any arbitrary starting value r
L+1,K=0
,asK →∞. This type of contraction condition is
sufficient, but not necessary, for convergence. In order to control convergence, we modify
the discretization of the acceleration term:
17
¨
r
L+1
≈
˙

. (3.22)
Inserting this into
m
¨
r = 
tot
(r) (3.23)
17
This collapses to a stencil of
¨
r
L+1
=
r
L+1
−2r
L
+r
L−1
(t)
2
when the time step size is uniform.
n05 book
2007/5/15
page 24
✐
✐
✐
✐
✐

L


 
R
, (3.24)
whose convergence is restricted by
η ∝ EIG(G) ∝
t
2
m
. (3.25)
Therefore, we see that the eigenvalues of G are (1) directly dependent on the strength of
the interaction forces, (2) inversely proportional to the mass, and (3) directly proportional
to (t )
2
(at time = t). Therefore, if convergence is slow within a time step, the time step
size, which is adjustable, can be reduced by an appropriate amount to increase the rate of
convergence. Thus, decreasing the time step size improves the convergence; however, we
want to simultaneously maximize the time step sizes to decrease overall computing time
while still meeting an error tolerance on the numerical solution’s accuracy. In order to
achieve this goal, we follow an approach found in Zohdi [208], [209], originally developed
for continuum thermochemical multifield problems in which (1) one approximates
η
L+1,K
≈ S(t)
p
(3.26)
(S is a constant) and (2) one assumes that the error within an iteration behaves according to
(S(t)

is not met in the desired number of iterations, the contraction constant η
L+1,K
is too large.
Accordingly, one can solve for a new smaller step size under the assumption that S is
constant:
t
tol
= t




TOL
||
L+1,0
||

1
pK
d

||
L+1,K
||
||
L+1,0
||

1
pK

✐
✐
3.3. Application to particulate flows 25
Remark. Classical solution methods require O(N
3
) operations, whereas iterative
schemes, such as the one presented, typically require order N
q
, where 1 ≤ q ≤ 2. For
details, see Axelsson [11]. Also, such solvers are highly advantageous, since solutions to
previous time steps can be used as the first guess to accelerate the solution procedure.
Remark. A recursive iterative scheme of Jacobi type, where the updates are made
only after one complete system iteration, was illustrated here only for algebraic simplicity.
The Jacobi method is easier to address theoretically, while the Gauss–Seidel method, which
involves immediately using the most current values, when they become available, is usually
used at the implementation level. As is well known, under relatively general conditions, if
the Jacobi method converges, the Gauss–Seidel method converges at a faster rate, while if
the Jacobi method diverges, the Gauss–Seidel method diverges at a faster rate (for example,
see Ames [5] or Axelsson [11]). The iterative approach presented can also be considered
as a type of staggering scheme. Staggering schemes have a long history in the computa-
tional mechanics community. For example, see Park and Felippa [161], Zienkiewicz [206],
Schrefler [173], Lewis et al. [133], Doltsinis [52], [53], Piperno [162], Lewis and Schrefler
[132], Armero and Simo [7]–[9], Armero [10], Le Tallec and Mouro [131], Zohdi [208],
[209], and the extensive works of Farhat and coworkers (Piperno et al. [163], Farhat et al.
[65], Lesoinne and Farhat [130], Farhat and Lesoinne [66], Piperno and Farhat [164], and
Farhat et al. [67]).
Remark. It is important to realize that the Jacobi method is perfectly parallelizable.
In other words, the calculations for each particle are uncoupled, with the updates only
coming afterward. Gauss–Seidel, since it requires the most current updates, couples the
particle calculations immediately. However, these methods can be combined to create

i
(r) = 
EXT
i
(r) + 
INT
i
(r), (3.31)
to obtain
N
p

i=1
m
i
¨
r
i
=
N
p

i=1
(
EXT
i
(r) + 
INT
i
(r)) =

✐
26 Chapter 3. Iterative solution schemes
Thus, a consistency check can be made by tracking the condition






N
p

i=1

INT
i
(r)






= 0. (3.33)
This condition is usually satisfied, to an extremely high level of accuracy, by the previously
presented temporally adaptive scheme. However, clearly, this is only a necessary, but not
sufficient, condition for zero error.
Remark. An alternative solution scheme would be to attempt to compute the solution
by applying a gradient-based method like Newton’s method. However, for the class of
systems under consideration, there are difficulties with such an approach.

− (A
TAN ,K−1
)
−1
(r
K−1
), (3.37)
where
A
TAN ,K
=
(
∇
r
A(r)
)
|
r
K
=
(
∇
r
(r)
)
|
r
K
(3.38)
is the tangent. Therefore, in the fixed-point form, one has the operator

✐
3.4. Algorithmic implementation 27
(1) GLOBAL FIXED-POINT ITERATION (SET i = 1 AND K = 0):
(2) IF i>N
p
, THEN GO TO (4);
(3) IF i ≤ N
p
, THEN
(a) COMPUTE POSITION: r
L+1,K
i
≈
t
2
m
i


tot
i
(r
L+1,K−1
)

+ r
L
i
+ t
˙

i
||
(normalized);
(b) Z
K
def
=

K
TOL
r
;
(c) 
K
def
=


(
TOL

0
)
1
pK
d
(

K


error measures wereused. As with theunnormalized case, oneapproximates the error within
an iteration to behave according to
(S(t)
p
)
K
||r
L+1,1
− r
L+1,0
||
||r
L+1,0
− r
L
||

 

0
=
||r
L+1,K
− r
L+1,K−1
||
||r
L+1,K
− r
L

✐
✐
28 Chapter 3. Iterative solution schemes
γ , into the iterations:
r
L+1,K
= γ(G(r
L+1,K−1
) + R) + (1 − γ)r
L+1,K−1
. (3.41)
Since the scheme must reproduce the exact solution, we have
r
L+1
= γ(G(r
L+1
) + R) + (1 − γ)r
L+1
. (3.42)
Subtracting Equation (3.42) from Equation (3.41) yields
r
L+1,K
− r
L+1
= γ

G(r
L+1,K−1
) − G(r
L+1

2
, ,a
K
, ,a implies
a − a
K+1
a − a
K
= <1, (3.45)
where  is a constant and a is the limit. Now consider the following sequence of terms:
a ≈ a
K
+ C
K
⇒ a − a
K
≈ C
K
,
⇒ a − a
K+1
≈ C
K+1
= (a − a
K
),
⇒ a − a
K+2
≈ C
K+2

K
− (a
K
)
2
a
K+2
+ a
K
− 2a
K+1
. (3.48)
We then repeat the procedure on the newly generated sequence:
a
K,2
=
a
K+2,1
a
K,1
− (a
K,1
)
2
a
K+2,1
+ a
K,1
− 2a
K+1,1

✐
✐
05 book
2007/5/15
page 31
✐
✐
✐
✐
✐
✐
✐
✐
Chapter 4
Representative numerical
simulations
In order to illustrate how to simulate a particulate flow, we consider a group of N
p
randomly
positioned particles ina cubical domain withdimensions D×D×D. During the simulation,
if a particle escapes from the control volume, the position component is reversed and the
velocity component is retained (now incoming). Thus, for example, if the x
1
component
of the position vector for the ith particle exceeds the boundary of the control volume, then
r
ix
1
=−r
ix

3
3
. (4.3)
Thus, the total volume occupied by the particles, denoted by , can be written as
ν = v
f
N
p
V, (4.4)
and the total mass is
M =
N
p

i=1
m
i
= ρν, (4.5)
20
There are many variants of this procedure.
21
D is normalized to unity in these simulations.
31
05 book
2007/5/15
page 32
✐
✐
✐
✐

however, such methods appear to offer distinct computational advantages if extremely high
volume fractions are desired.
4.1 Simulation parameters
The relevant simulation parameters were
• number of particles = 100,
• (normalized) box dimension D = 1m,
• initial mean velocity field = (1.0, 0.1, 0.1) m/s,
• initial random perturbations around mean velocity = (±1.0, ±0.1, ±0.1) m/s,
• (normalized) length scale of the particles, L = 0.25, with corresponding volume
fraction v
f
=
4πL
3
3
= 0.0655 and radius b = 0.0539 m,
• mass density of the particles = 2000 kg/m
3
,
• simulation duration = 1s,
• initial time step size = 0.001 s,
• time step upper bound = 0.01 s,
• tolerance for the fixed-point iteration = 10
−6
.
The parameters
α
1
and α
2

0.2
0.4
0.6
0.8
Z
0.2
0.4
0.6
0.8
Figure 4.1. A typical starting configuration for the types of particulate systems
under consideration.
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ENERGY FRACTION
TIME
RELATIVE MOTION
CENTER OF MASS MOTION
0.25
0.3
0.35
0.4
0.45

0.4
0.45
0.5
0.55
0.6
0.65
0.7
0.75
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ENERGY FRACTION
TIME
RELATIVE MOTION
CENTER OF MASS MOTION
Figure 4.2. The proportions of the kinetic energy that are bulk and relative motion.
Top to bottom and left to right, for e
o
= 0.5, µ
s
= 0.2, µ
d
= 0.1: (1) no near-field
interaction, (2)
α
1
= 0.1 and α
2
= 0.05, (3) α
1
= 0.25 and α
2

0.59
0.6
0.61
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ENERGY (N-m)
TIME
TOTAL KINETIC ENERGY
0.59
0.6
0.61
0.62
0.63
0.64
0.65
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ENERGY (N-m)
TIME
TOTAL KINETIC ENERGY
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ENERGY (N-m)
TIME
TOTAL KINETIC ENERGY

1
= 0.5 and α
2
= 0.25 (Zohdi
[212]).
tive kinetic energy in the system, is dramatically different with increasing severity of the
near-field forces.
22
Notice that the kinetic energy per unit mass is nonmonotone when the
near-field interactions are taken into account (Figure 4.3). One may observe that, from
Figure 4.2, as the near-field strength is increased, the component of the kinetic energy cor-
responding to the relative motion does not decay and actually becomes dominant with time.
Essentially, the near-field interaction becomes strong enough that the flowing system expe-
riences a transition to a vibrating ensemble. This transition can be qualitatively examined
by recognizing that the governing equations are formally similar to classical, normalized,
linear (or linearized) second-order equations governing a one degree of freedom harmonic
oscillator of the form
¨r +2ζω
n
˙r +ω
2
n
r =
f(t)
m
, (4.7)
where
ω
n
=

d
def
=
2π
ω
d
, (4.10)
where
ω
d
def
= ω
n

1 − ζ
2
(4.11)
is the “dampednatural frequency.” Using standard procedures, one decomposes thesolution
into homogeneous and particular parts:
r = r
H
+ r
P
. (4.12)
The homogeneous part must satisfy
¨r
H
+ 2ζω
n
˙r

λ
1,2
= ω
n
(−ζ ±

ζ
2
− 1). (4.17)
The general solution is
r = A
1
exp(λ
1
t) +A
2
exp(λ
2
t). (4.18)
Depending on the value of ζ , the solution will have one of three distinct types of behavior:
• ζ>1, overdamped, leading to no oscillation, where the value of r approaches zero
for large values of time. Mathematically, λ
1
and λ
2
are negative numbers, so
r
H
= A
1

= (A
1
+ A
2
t)exp(ω
n
t). (4.20)
05 book
2007/5/15
page 36
✐
✐
✐
✐
✐
✐
✐
✐
36 Chapter 4. Representative numerical simulations
• ζ<1, underdamped, leading to damped oscillation, where the value of r approaches
zero for large values of time, in an oscillatory fashion. Mathematically, ζ
2
− 1 < 0,
so
r
H
= A
1
cos(ω
d

= R sin(t − φ), (4.24)
where
R =
f
o
k


1 −

2
ω
2
n

2
+

2ζ

ω
n

2
(4.25)
and
φ = tan
−1

2ζ

2
r
−β
2
(4.28)
and d is an effective dissipation term that would arise from inelastic impact and friction.
Upon linearization of the nonlinear interaction relation about a point r
∗
,

nf
(r) ≈ 
nf
(r
∗
) +
∂
nf
∂r




r=r
∗
(r − r
∗
) + O(r − r
∗
), (4.29)

4.2. Results and observations 37
where
ω
∗
n
=

−
∂
nf
∂r
|
r=r
∗
m
, (4.31)
ζ
∗
=
d
2mω
∗
n
, (4.32)
and
f
∗
(t) = 
nf
(r

2
β
2
r
−β
2
−1
∗
m
=

−α
1
mβ
1
r
−β
1
−1
∗
+ α
2
mβ
2
r
−β
2
−1
∗
, (4.34)

2
−1
∗
. (4.35)
We note that if the parameters are chosen (as in the preceding simulations) specifically as
(β
1
,β
2
) = (1, 2) and r
∗
is chosen as the static equilibrium point, r
e
, given by Equation
(1.36), then
r
∗
= r
e
=
α
2
α
1
(4.36)
and
ω
∗
n
=

∗
m
, (4.37)
where
k
∗
def
= α
1

α
1
α
2

2
. (4.38)
Thus, in the preceding numerical examples, when we kept the ratio
α
1
α
2
constant, but in-
creased α
1
(while keeping m constant), we were effectively increasing the “stiffness” in the
system and, therefore, the amount of (pre)stored energy available to counteract dissipation.
Clearly, under certain conditions, a particulate flow may “pulse” (oscillate) depending on
the character of the interaction and the contact parameters. Thus, oscillatory behavior is not
unexpected for the multibody system (Figure 4.3). We remark that increasingly smaller ω