05 book
2007/5/15
page 58
✐
✐
✐
✐
✐
✐
✐
✐
58 Chapter 7. Advanced particulate flow models
d
c
Figure 7.3. Introduction of a cutoff function.
Thus, the preceding analysis indicates that, for the three-dimensional case, an interaction
“cutoff” distance (d
c
) should be introduced (Figure 7.3),
||r
i
− r
j
|| = d
c
≤ d
(+)
, (7.14)
to avoid long-range (central-force) instabilities.
Remark. By introducing a cutoff distance, one can circumvent a loss-of-convexity
instability. However, introducing such acutoff can induce another type of instability. Specif-

2
β
1
α
1

1
−β
1
+β
2
. (7.15)
Clearly, since β
2
>β
1
, d
(−)
is a lower bound(dictatedbythe minimum interaction distance),
while d
(+)
is an upper bound (dictated by the (convexity-type) stability).
7.4 A simple model for thermochemical coupling
As indicated earlier, in certain applications, in addition to the near-field and contact effects
introduced thus far, thermal behavior is of interest. For example, applications arise in the
study of interstellar particulate dust flows in the presence of dilute hydrogen-rich gas. In
many cases, the source of heat generated during impact in such flows can be traced to the
reactivity of the particles. This affects the mechanics of impact, for example, due to thermal
softening. For instance, the presence of a reactive substance (gas) adsorbed onto the surface
of interplanetary dust can be a source of intense heat generation, through thermochemical

o

1 −
v
n
v
∗

,e
−

max

1 −
θ
θ
∗

, 0

, (7.16)
where θ
∗
can be considered as a thermal softening temperature. In order to determine the
thermal state of the particles, we shall decompose the heat generation and heat transfer
processes into two stages. Stage I describes the extremely short time interval when impact
occurs, δt  t, and accounts for the effects of chemical reactions, which are relevant in
certain applications, and energy release due to mechanical straining. Stage II accounts for
the postimpact behavior involving convective and radiative effects.
7.4.1 Stage I: An energy balance during impact

particles deform negligibly during impact, we assume that there is an insignificant amount
of mechanically stored energy. The kinetic energy is
K =
1
2
mv ·v. (7.19)
The mechanical power term is due to the forces acting on a particle, namely
P =
dW
dt
= 
tot
· v, (7.20)
and, because
dK
dt
= m
˙
v ·v, (7.21)
and we have a balance of momentum
m
˙
v ·v = 
tot
· v, (7.22)
we have
dK
dt
=
dW

|
I
∗
n
, 1

πb
2
, (7.27)
where I
n
is the normal impact force; κ is a reaction (saturation) constant, energy per unit
area; I
∗
n
is a normalization parameter; and b is the particle radius. For details, see Schmidt
[172], for example. For the grain sizes and material properties of interest, the term in
Equation (7.26),
δH
mC
, indicates that values of approximately κ ≈ 10
6
J/m
2
will generate
05 book
2007/5/15
page 61
✐
✐

a
s
=
b
3
, which indicates that a
uniform temperature distribution is appropriate, since the particles, by definition, are small.
We also assume that the dynamics of the (dilute) gas does not affect the motion of the (much
heavier) particles. The gas only supplies a reactive thin film on the particles’ surfaces. The
first law reads
d(K +U)
dt
= m
˙
v ·v + mC
˙
θ = 
tot
· v

 
mechanical power
− h
c
a
s
(θ − θ
o
)


is the heating due to convection (Newton’s law of cooling) into
the dilute gas; and a
s
is the surface area of a particle. It is assumed that the radiation
exchange between the particles (emission exchange between particles) is negligible.
32
For
the applications considered here, typically, h
c
is quite small and plays a small role in the
heat transfer processes.
33
From a balance of momentum, we have m
˙
v · v = 
tot
· v, and
Equation (7.28) becomes
mC
˙
θ =−h
c
a
s
(θ − θ
o
) − Ba
s
(θ
4

c
a
s
t
¯
θ(t) −
tBa
s

mC + h
c
a
s
t

θ
4
(t +t) −θ
4
s

+
h
c
a
s
tθ
o
mC + h
c

allows the calculation of the amount of radiation emitted
in all directions and over all wavelengths simply from the knowledge of the temperature of
the black-body. We note that Equation (7.30) is of the form
θ(t + t) = G(θ (t + t)) + R, (7.31)
where R = R(θ(t +t)) and G’s behavior is controlled by
tBa
s

mC + h
c
a
s
t
, (7.32)
which is quite small. Thus, a fixed-point iterative scheme such as
θ
K
(t +t) = G(θ
K−1
(t +t)) +R (7.33)
would converge rapidly.
34
For this stage, since δt  t, we assign θ(t) = θ(t + δt) = θ(t) +
δH
mC
and replace θ(t) with it in Equation
(7.30).
35
Radiation is idealized as requiring no medium to transmit energy.
05 book

of computational effort, since within a time step, many multifield system re-solves can take
place. We now develop a staggering scheme by following an approach found in Zohdi
[208]–[210].
7.5.1 A general iterative framework
We consider an abstract setting, whereby one solves for the particle positions, assuming the
thermal fields fixed,
A
1
(r
L+1,K
,θ
L+1,K−1
) = F
1
(r
L+1,K−1
,θ
L+1,K−1
), (7.34)
and then one solves for the thermal fields, assuming the particle positions fixed,
A
2
(r
L+1,K
,θ
L+1,K
) = F
2
(r
L+1,K

||
||θ
L+1,K
− θ
L
||
.
(7.36)
We define the maximum “violation ratio,” i.e., the larger of the ratios of each field variable’s
error to its corresponding tolerance, by Z
K
def
= max(z
rK
,z
θK
), where
z
rK
def
=

rK
TOL
r
and z
θK
def
=




TOL
r

r0

1
pK
d


rK

r0

1
pK


,φ
θK
def
=



TOL
θ


= A(w) − F . (7.40)
Linearization leads to
(w
K
) = (w
K−1
) +∇
w
|
w
K−1
(w
K
− w
K−1
) + O(||w||
2
), (7.41)
and thus the Newton updating scheme can be developed by enforcing
(w
K
) ≈ 0, (7.42)
leading to
w
K
= w
K−1
− (A
TAN ,K−1
)

−1
(w). (7.45)
One immediately sees a fundamental difficulty due to the possibility of a zero or near-zero
tangent when employing a Newton’s method on a nonconvex system, whichcan lose positive
definiteness and which in turn will lead to an indefinite system of algebraic equations.
36
Therefore, while Newton’s method usually converges at a faster rate than a direct fixed-
point iteration, quadratically as opposed to superlinearly, its convergence criteria are less
robust than the presented fixed-point algorithm, due to its dependence on the gradients of
the solution. Furthermore, for the problems considered, the solutions are nonsmooth and
nonconvex, primarily due to the impact events, and thus we opted for the more robust
“gradientless” staggering scheme.
36
Furthermore, the tangent may not exist in some (nonsmooth) cases.
05 book
2007/5/15
page 65
✐
✐
✐
✐
✐
✐
✐
✐
7.5. Staggering schemes 65
(1) GLOBAL FIXED-POINT ITERATION (SET i = 1 AND K = 0):
(2) IF i>N
p
, THEN GO TO (4);

= θ
L
i
+
δH
L+1,K−1
mC
;
θ
L+1,K
i
=
mC
mC + h
c
a
s
t
θ
L+1,K−1
i
−
tBa
s

mC + h
c
a
s
t

N
p
i=1
||r
L+1,K
i
− r
L+1,K−1
i
||

N
p
i=1
||r
L+1,K
i
− r
L
i
||
,
θK
def
=

N
p
i=1
||θ

rK
TOL
r
,z
θK
def
=

θK
TOL
θ
;
(c) 
K
def
= min(φ
rK
,φ
θK
) where φ
rK
def
=




TOL
r


1
pK
d


θK

θ0

1
pK



;
(5) IF TOLERANCE MET (Z
K
≤ 1) AND K<K
d
, THEN
(a) INCREMENT TIME: t = t + t;
(b) CONSTRUCT NEW TIME STEP: t = 
K
t;
(c) SELECT MINIMUM: t = min(t
lim
,t)AND GO TO (1);
(6) IF TOLERANCE NOT MET (Z
K
> 1) AND K = K

2
+ w
3
= 0,
c ¨w
3
+ w
1
= 0.
(7.46)
When this is discretized in time, for example, with a backward Euler scheme, we obtain
˙w
1
L+1
=
w
L+1
1
− w
L
1
t
,
˙w
2
L+1
=
w
L+1
2

01
t
b
(t)
2
c
01








w
L+1
1
w
L+1
2
w
L+1
3





=

010
001








w
L+1,K
1
w
L+1,K
2
w
L+1,K
3





=






t
b
w
L+1,K−1
2
(t)
2
c
w
L+1,K−1
3







. (7.49)
Rewriting this in terms of the standard fixed-point form, G(w
L+1,K−1
) + R = w
L+1,K
,
yields



0
t





 
w
L+1,K−1
+





w
L
1
w
L
2
2w
L
3
− w
L−1
3






ditions, 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] for details. The
Jacobi method is easier to address theoretically, so it is used for proof of convergence, and the Gauss–Seidel
method is used at the implementation level.
05 book
2007/5/15
page 67
✐
✐
✐
✐
✐
✐
✐
✐
7.5. Staggering schemes 67
The eigenvalues of G are given by λ
3
=
(t)
4
abc
and, hence, for convergence we must have
|max λ|=




(t)
4

0
t
a
t
b
0


 
G

w
L+1,K−1
1
w
L+1,K−1
2


 
w
L+1,K−1
+

w
L
1
w
L
2

L+1
1
=
abw
L
1
+ btw
L
2
ab −(t)
2
= w
L
1
−
w
L
2
a
t

 
first staggered iteration
+
w
L
1
ab
(t)
2

L
2
ab
(t)
2
  
second staggered iteration
+···.
(7.56)
As pointed out in Zohdi [208], the time step induced restriction for convergence matches
the radius of analyticity of a Taylor series expansion of the solution around time t, which
converges in a ball of radius from the point of expansion to the nearest singularity in
Equations (7.55) and (7.56). In other words, the limiting step size is given by setting the
denominator to zero,
ab −(t)
2
= 0, (7.57)
which is in agreement with the condition derived from the analysis of the eigenvalues of G.
05 book
2007/5/15
page 68
✐
✐
✐
✐
✐
✐
✐
✐
68 Chapter 7. Advanced particulate flow models

v
f
def
=
4πL
3
3
. Thus, thetotal volume occupied by the particles, denoted by ν, could be written
as ν = v
f
N
p
V , and the total mass could be written as M =

N
p
i=1
m
i
= ρν, while that of
an individual particle, assuming that all are the same size, was m
i
=
ρν
N
p
= ρ
4
3
πb

parameter in the coefficient of restitution relation was θ
∗
= 3000
◦
K. The reaction constant
was varied in the range 10
6
J/m
2
≤ κ ≤ 10
7
J/m
2
, with I
∗
= 10
3
N. The coefficient of
convective heat transfer (h
c
) was set to zero. We introduced the following near-field param-
eters per unit mass
2
: α
1ij
= α
1
m
i
m

−3
s.
The tolerances of both fields (TOL
r
and TOL
θ
) for the fixed-point iterations were set to 10
−6
and the upper limit on the number of fixed-point iterations was set to K
d
= 10
2
.
Two main types of computational tests were conducted:
1. varying κ, for a given field strength,
α
1
= 0.5 and α
2
= 0.25, with a clustering
augmentation of
α
a
= 1.75 (forcing a small gap characterized by d
a
= 1.03(2b)),
β
a
= 1, δ
a

Z
XY
Z
XY
Z
Figure 7.5. Top to bottom and left to right, the dynamics of the particulate flow
with clustering forces: An initially fine cloud of particles that clusters to form structures
within the flow. Blue indicates a temperature of approximately 300
◦
K, while red indicates
a temperature of approximately 400
◦
K (Zohdi [217]).
For each different parameter selection, the initial conditions, i.e., random positions,
velocities, temperatures, etc., were the same. We remark that parameter studies on the near-
field strength, in isolation (without thermochemical coupling), havebeenconductedinZohdi
[209]. The field strength chosen was strong enough to induce vibratory motion and hence
nonmonotone kinetic energy. Frames of the flows for cases 1 and 2, for (typical) values of
κ = 2 ×10
6
J/m
2
, are shown in Figures 7.5 and 7.6. The plots in Figures 7.7–7.10 indicate
the overall energetic and thermal behavior. Typically, the simulations took approximately
between 1 min and 2 min on a standard (Dell, 2.33 GHz) laptop.
39
For the parameter ranges
used in the presented simulations, the degree of adaptivity needed strongly depended on the
presence of the clustering forces, and to a somewhat lesser degree on the thermochemical
parameters. For example, for the 5-s simulation, if the time steps stayed at the starting value

(t = 10
−3
s), then 5000 timestepswouldbe needed if therehadbeenno time step adaptivity
(time step enlargement). Conversely, if the time steps were found to be unnecessarily small
(an overkill) at the starting value (t = 10
−3
s), and, consequently, unrefined to the upper
bound (t
lim
= 10
−2
s), then approximately500 timesteps wouldbe needed. Tables 7.1 and
7.2 indicate that, for the parameter ranges tested, when clustering forces were not present,
the time steps did not need to be refined or unrefined. However, when clustering forces were
present, the time steps could be unrefined for the given tolerances, requiring more internal
fixed-point iterations. This was primarily because cluster structures formed, leading to
fewer collisions between the larger objects, which did not require such small time steps
(Figure 7.11). For the simulations with clustering forces, there was an expected thermal
sensitivity. As the reaction constant κ became stronger, the number of fixed-point iterations
required to achieve convergence increased. These results highlight an essential point of the
adaptive time-stepping process, which is to allow the system to adjust to the physics of the
problem. Some furtherremarks elaborating on this issuecan be found in Zohdi [208]–[210].
05 book
2007/5/15
page 71
✐
✐
✐
✐
✐

0.5
0.55
0.6
0.65
0.7
0.75
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
ENERGY (N-m)
TIME
TOTAL KINETIC ENERGY
0.4
0.45
0.5
0.55
0.6
0.65
0.7
0.75
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
ENERGY (N-m)
TIME
TOTAL KINETIC ENERGY
Figure 7.7. Top to bottom and left to right, with clustering forces: the total kinetic
energy in the system per unit mass with e
o
= 0.5, µ
s
= 0.2, µ
d
= 0.1, α

through the coefficient of restitution, when clustering dominates, the particle dynamics will
be only marginally affected by varying κ (Figure 7.7). However, the temperature of the
particles in the presence of clustering will rise substantially, due to the large compressive
forces between the contacting particles, which activate the chemical reactions. Also, we
remark that the group dynamics, for different κ without clustering forces, deviate much
more from one another than the cases when clustering is present (Figure 7.8). Typically,
when two particles have clustered, since the binding field was strong, the particles rarely
become dislodged. This issue has been been investigated in depth in Zohdi [225].
05 book
2007/5/15
page 72
✐
✐
✐
✐
✐
✐
✐
✐
72 Chapter 7. Advanced particulate flow models
0.6
0.61
0.62
0.63
0.64
0.65
0.66
0.67
0.68
0.69

0.7
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
ENERGY (N-m)
TIME
TOTAL KINETIC ENERGY
0.6
0.61
0.62
0.63
0.64
0.65
0.66
0.67
0.68
0.69
0.7
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
ENERGY (N-m)
TIME
TOTAL KINETIC ENERGY
Figure 7.8. Top to bottom and left to right, without clustering forces: the total
kinetic energy in the system per unit mass with e
o
= 0.5, µ
s
= 0.2, µ
d
= 0.1, α
1
= 0.5,

def
=
v
obn
(t +δt) − v
in
(t +δt)
v
in
(t) − v
obn
(t)
,
(7.58)
where v
obn
remains the same before and after impact. In Figure 7.13, the impact of a cloud
against an obstacle is shown.
40
Let us focus on a particle impacting a massive obstacle
M  m (Figure 7.12). A balance of momentum reads for the particle as
mv(t) −
¯
Iδt ±|
¯
F |δt = mv(t + δt). (7.59)
40
All other parameters are the same as in the previous simulations.
05 book
2007/5/15

320
325
330
335
340
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
TEMPERATURE
TIME
TEMPERATURE
300
400
500
600
700
800
900
1000
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
TEMPERATURE
TIME
TEMPERATURE
200
400
600
800
1000
1200
1400
1600
1800

J/m
2
(Zohdi [217]).
The coefficient of restitution reads as
e
def
=
−v
in
(t +δt)
v
in
(t)
,
(7.60)
so
¯
I =−
m(v(t +δt) −v(t)
δt
±|
¯
F |=−
mv(t)(1 +e)
δt
±|
¯
F |, (7.61)
where ±|
¯

302
303
304
305
306
307
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
TEMPERATURE
TIME
TEMPERATURE
298
300
302
304
306
308
310
312
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
TEMPERATURE
TIME
TEMPERATURE
300
305
310
315
320
325
330
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

κ = 10
6
J/m
2
, (2) κ = 2 × 10
6
J/m
2
, (3) κ = 4 × 10
6
J/m
2
, and (4) κ = 8 × 10
6
J/m
2
(Zohdi [217]).
Table 7.1. The number of time steps and fixed-point iterations, with clustering
forces: the average particle temperature with e
o
= 0.5, µ
s
= 0.2, µ
d
= 0.1, α
1
= 0.5, and
α
2
= 0.25.

✐
7.5. Staggering schemes 75
Table 7.2. The number of time steps and fixed-point iterations, without clustering
forces: the average particle temperature with e
o
= 0.5, µ
s
= 0.2, µ
d
= 0.1, α
1
= 0.5, and
α
2
= 0.25.
κ(J ×10
6
/m
2
)
Time Steps Fixed-Point Iterations
1 5000 5025
2 5000 5024
4
5000 5029
8 5000 5024
XY
Z
Figure 7.11. A zoom on the structures that form with clustering. Blue indicates a
temperature of approximately 300

Y
Z
Figure 7.13. Top to bottom and leftto right, a charged cloud against animmovable
obstacle.