7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Particulate Flows
Rodolphe Prignitz and Eberhard Binsch
Chair of Applied Mathematics 3, University of ErlangenNimrberg, Haberstr. 2, 91058 Erlangen, Germany
prignitz @ am.unierlangen.de and baensch @ am.unierlangen.de
Keywords: CFD, multiphase flows, particulate flow, finite elements, subspace projection, rheology
Abstract
The flow of particles in a fluid can be found in numerous industrial settings which use sedimentation, fluidisation
and lubricated transport such as food processing, catalytic processing, slurries, coating, paper manufacturing, particle
injection molding and filter operation. The ability to understand rheology effects of particulate flows is important for
the design, operation and efficiency of the underlying processes. Despite the fact that particle technology is widely
used, still it is an enormous experimental challenge to determine the correct parameters for the process used.
For the simulation of such processes we consider a Newtonian fluid with rigid particles in our numerical model.
The particles can interact through molecular or electrical forces. The method uses a finiteelement discretisation
in space and an operatorsplitting technique for discretisation in time and has its basis in work by Glowinski
et al. (1999). However, rather than using a Lagrange multiplier, a subspace projection is used to couple the
particle motion with the fluid motion. The advantage of this approach is the implicit treatment of the particle
and the fluid velocity. Combined with local mesh refinement, the method results in a fast and accurate algorithm
for the simulation of a huge number of particles in a flow field. Validation is achieved using the sedimentation of
one particle and comparing the resulting drag coefficient with theoretical and experimental results in three dimensions.
Introduction
Mathematical Model
The aim of the present paper is to analyze and develop
an algorithm for the simulation of a huge amount of par
ticles suspended in a viscous liquid. Related work can be
found for instance in Glowinski et al. (1999), Feng et al.
(1994), BOnisch and Heuveline (2007), Hu (1996) and
Wan and Turek (2007). We adapted the work of Glowin
ski et al. (1999), such that it fits our requirements to ob
tain a fast and accurate code for the simulation of a huge
number of particles in a viscous fluid. Our own con
tribution is the representation of the particles utilizing
adaptive finite element techniques, the time discretiza
tion as well as the incorporation of a subspace projection
method. In this paper we report on 3 dimensional simu
lations. This is an extension to 2 dimensional simulation
which we published in Prignitz and Bansch (2010) The
remainder of the paper is organized as follows. In the
first sections we present a short overview of the mathe
matical model and numerical method used, while in the
last section we present some numerical results.
In this section we introduce a model for particulate
flows. For the ease of presentation we restrict our
selves to the 3dcase with one particle. The exten
sion to more particles is straightforward. We denote
by Q(t) C R3 the area occupied by the fluid with ho
mogeneous Dirichlet boundary condition on its outer
boundary FD. P(t) C R3 is the particle and its center
Figure 1: Fluid domain with arbitrary particle inside.
of mass is denoted by X = 1 j(t) x cb, while
r x X is its relative coordinate. It should be empha
sized that the fluid area and the particle area don't inter
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
sect, Q(t)nP(t) = 0. The unknowns in the fluid area are
the velocity u and the pressure p, which are described by
the NavierStokesequations. The particle, being a rigid
body, is described by Newton's law. The describing val
ues are the translatorial and angular velocities U, w, re
spectively, the position X and the orientation in space
given by the angle e. The system in nondimensional
form reads
tu + u Vu + V pi D[u]) f in Q (t),
(Pff 1 Re
(1)
Vu 0 in (t), (2)
u 0 on FD, (3)
u = U + x r on OP(t), (4)
MU =F+j ands,
I + uc x (I ) P r x (o nds,
P(t)
1 Dtu v + u Vu v
1
pVv 1D[u] : D[v] c
2Re
+ (1 a) MU V+
(1 a) [Ij + w x (Iw)] .
f vcb+ F.V,
SV uqak =0,
2
6= e .
Hereby a is the density fraction defined by a = p/pp.
(6) In order to obtain this weak formulation one has to per
form the symbolic calculation
see for example Feng et al. (1994). The two integrals on
the right hand side represent the force and the torque, re
spectively, exerted by the fluid on the particle. Here, a is
the stress tensor, M denotes the mass of the particle and
I its inertia. F is an external force acting on the particle,
Re denotes the Reynolds number and the deformation
tensor is defined by D [u] j .= ., + ) .,
Following an idea of Glowinski et al. (1999), we de
rive a weak formulation of the problem in the form of
the one domain approach. To this end, the space of com
bined velocities is defined by
Hc() { (v, V, ) v ve (H1())3,
VER3,eRv E onFD,
v V+ xr inP(t)}. (9)
Note that the combined fluid/particle domain 2 =
S() U P (t) does not depend on time t. The space He (Q)
is called the space of combined velocities, since it con
tains all velocities posed in this problem, the fluid ve
locity, the translational particle velocity and the angular
velocity represented by v, V and respectively. We em
phasize that the fluid velocity is defined on the whole
domain 2 and is restricted to the particle velocity inside
the particle by the above definition. With test functions
(v, V, >) e HK(Q) and q e L2 (Q) we can state the weak
formulation of our problem (1)(8) as
resulting in an intermediate weak formulation of the
problem
S tu.v+uVuv acc+ Vpv 1 VD[u]vbc
n(t) M.n(t) Re
+ MU V + Id; ( + w x (Iiw) 
12(t)
f .vcb+ F V +T
+ o. nV+r x o n ds. (14)
ap(t)
Integration the stress tensor term by parts and the use of
v V + x r on OP(t) leads to
j Vpv VD[u]vc
pV v 1 D[u] : D[v] c
Jn(t) 2Re
SDpneD[u] n) (V + x r) (15)
.n
With a (b x c) = (c x a) b we receive for the boundary
term in (15)
.n.(V + x r) ds = j
Jap(t)
IP(t)
1 (1) .v + (5) V + (7) ,
0(nt)}
.nV+rxan. ds,
which is exactly the last term in the weak formulation
(14), but with different sign and hence vanishes. This is
the first advance of this method, that one does not need
to evaluate the stress tensor. This is caused by the special
choice of the functional space He.
Our set of formulas now reads
S tu v+u Vu v c
n(t)
S pV. 1 D[u] : D[v] cb
J(t) D[2Re
+ MU V + Id. + C x (j1).
Sf vd+ F V+ T (16)
J V.uqck 0, (17)
Q (t)
X U, (18)
( w. (19)
This integrals are defined on the time depended domain
Q(t). For our numerical treatment (one domain ap
proach) we want to avoid this and extend this set of for
mulas to the whole constant domain Q. For the time
derivative of the velocity we obtain
u = dtu + u Vu.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
calculate
/P t u v + u Vu v d
JP(t) Jp(t)
JP( (V + x r) (+
Px()
S(W x r + w x (w x r)) (V + x r) &d =
P(t)
JP(t)
U. V c+
JP(t) (cW x r + u x (u x r)) ( x r) cd ,
Jp(t)
SJP(t)
(cW x r + c x (c x r)) ( x r) cd
c JP(t)
aMU V + a(Id + w(xlw)) (24)
In addition we set f 0 in P(t). This will result in the
weak formulation presented in (10)(13)
For a shorter notation we define the linear forms
m(u, v) / u. v d, (25)
s(u, v) 2 D[u] : D[v] dc, (26)
2Re Qfi
k(u, v,w) / w Vu. vc, (27)
b(p, v) 2 pV v c, (28)
Inside the particle we demand u = U + w x r which
leads to
d
du U+j x r +wu x (w x r).
dt
Its easy to conduct V v = D[v] = 0 in P(t) Vv E
H,(Q). Let pp be the constant particle density and a
" the density ratio. For the dimensionless mass and
PP
inertia we have
M f= 1 (b
a P(t)
1i = l (11rl126,
a JP (t)
For a test function v E He and from (20) and (21) we
rirj) dx
L(U,V)
Si(w, )
S2 (Cll, CI2, )
a)MU V,
a) (i X (12))
a) (Cil X (JCI))) ,
yielding in the following formulation of the problem
m(dtu, v) + k(u, v, u) + s(u, v) b(p, v)
+ L(U, V) + S(, () + S,(W, W, )
=m(f,v)+FV, (33)
b(q, u) = 0,
S=U,
Q=e .
Numerical method
Our numerical scheme to solve the problem (33)(36) is
based on the following four steps
Splitting: in order to reduce the complexity of the
problem we separate the calculation of the parti
cle's time dependent location from the calculation
of the flow field.
Time discretization: a BDF2 projection scheme is
used for the efficient solution of the NavierStokes
equations.
Subspace projection: a novel method to take into
account the restriction of the function space He (().
Adaptivity: local (time dependent) grid refinement
is crucial to represent the particle's geometry.
Splitting
For an efficient numerical treatment we initially split the
position of the particle from the remaining variables.
1. Predictor: Solve for X, U, e, c using equations
(35), (36), MU F and Ij = c x (Iw).
2. Velocities:
m(dtu, v) + k(u, v, u) + s(u, v) b(p, v)
+L(U,V) +S(, E)=m(f,v), (37)
b(q, u) 0.
In the first step the positions (X, 0) and velocities (U,
w) of the particle are calculated, while in the second step
the velocities of the particle are corrected and the fluid
velocity and pressure is computed. The effect of the ex
ternal force F is only taken into account in the predictor
step. The force the fluid cause on the particle is created
implicitly by the terms L(U, V), S(9, I) and the usage
of the space of combined velocities He ().
Time discretization
For the time discretization of the predictor a velocity ver
let method with a = for the variables X and U is
used:
2
X+1 = X + TUk + a, (39)
2
2
a+1 = evaluate forces at new positions Xk+1, (41)
Uk+1 = Ukfc + Ta kc+1. (42)
For a more detailed description of the algorithm see for
example Martys and Mountain (1999). The same is valid
for the variables ( and cw.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The flow field governed by the NavierStokes equa
tions is solved by a BDF2 based projection method in
rotational form, see Guermond and Shen (2004). Intro
ducing the time step T, and 7 = this scheme is based
on three steps to solve equations (38):
1. Burgers problem
m(u+1, v) + yk(u+ 1, v, u+ 1) + _s(uk +, v)
2L(Uk+1, V) +2S1(w 1)
3 3
Sb(i.', v) + 7m(f(tk+l), v)
+ m(u, _ u v)
3 3
+ b( X _X v)
+ L(U, V) + 2S(w, ) (43)
3 3
2. Poisson problem
m(Vk+1, V') = b(t,, uk+1) (44)
"/
3. Update
m/' +1,q) =m(/' +Xk+ q)
b(q, u2 )
Re
(45)
In the calculations presented later, the form k is lin
earized by k(u +l,v, 2uk u 1).
Spatial discretization
The crucial point in the spatial discretization is to define
a discrete counterpart of He (() and, moreover, the con
crete realization of the this nonstandard finite element
space. A brief description of how to solve this problem
is given in the sequel. A more detailed presentation will
be published in near future.
Let T be a triangulation of L. Define the usual finite
element space by
X() { (v, V, ) v E (C())3
ve (P"(T))3VT E T,
VER3, I eR 3,V 0 onFD}.
A discrete subspace of Hc() is now given by
X,(Q) { (v, V, ,) e X()
vc V+ xr in P(t)}
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
For an arbitrary time step k the linearized equation (43)
may be rewritten with the bilinear form a, the corre
sponding operator A, and the cumulative right hand side
g: find u E Xc(Q) such that for all v E Xc(Q) it holds
a(u, v) =: (Au, v) (g, v). (46)
To circumvent the explicit representation of H6,(), a
subspace projection P : X Xc is used. With this
operator (46) may be formulated in terms of the standard
finite element space X(Q): find u E X(Q) such that for
all v E X(() it holds
(A'Pi, Pv)= (g, Pv). (47)
Note that the solution u is now easily found by taking
u = Pu, where u is a solution of equation (47). The
above system now leads to the linear system of equations
for the nodal vector U of the form
PTAPU = PTG, (48)
where A is the system matrix corresponding to opera
tor A and P is a matrix representation of P. Note that,
when using iterative solvers, one can bypass to explicitly
compute the modified system matrix PTAP, but rather
just slightly needs to modify the matrix vector product,
because one only has to take into account the action of
PTAP on a vector. We call this method subspace pro
jection method. Because the matrix P is quite simple, its
not necessary to store it explicitly. Instead, a short rou
tine can perform the multiplication of P and PT with a
vector v. This pseudocode shows that computation.
! Multiplication (u,U,omega)=P*(v,V,xi)
subroutine Pmul(v,V,xi,u,U,omega)
U, omega
do ii=l,npart Number of particles
U(:,ii) = V(:,ii)
omega(ii) = xi(ii)
end do
! u = rigid body motion in th<
do i=l,nk Number of DOFs
if( isparticle(i) ) then
ii= numpart(i)
r(:)= x(:,i) xpart(:,
u(l,i) = V(l,ii)
+ xi(2,ii)*r(3)
xi (3, ii) *r (2)
u(2,i) = V(2,ii)
+ xi(3,ii) *r (1)
xi (1, ii) *r (3)
u(3,i) = V(3,ii)
+ xi (1, ii) *r (2)
xi(2,ii) *r (1)
particle
ii)
end if
end do
nd subroutin
Adaptivity
Another important point is how to represent the parti
cle's geometry. In Hu (1996) a remeshing technique was
used to explicitly follow the geometry in time, Wan and
Turek (2007) introduced a mesh deformation technique
and Glowinski et al. (1999) used Lagrange multipliers.
All these techniques have some pros and cons, in partic
ular in 3d.
In contrast to the above mentioned methods, we use
time dependent adaptively refined grids based on the bi
section method by Bansch (1991) to sufficiently resolve
the region around the particle.
Figure 2: Adaptive refined mesh around a particle. For
an accurate representation its useful to refine the mesh
on the particle boundary.
The above algorithm was implemented in the finite el
ement flow solver NAVIER, for more details see Bansch
(1998).
Application to sedimentation
The test case presented here is the sedimentation of a
single particle in a lighter fluid. Though conceptually
rather simple, this test case might be already quite re
vealing. Here, the quantity of interest is the terminal
particle velocity, denoted by U and the corresponding
drag CD defined by
1
F PCDdU2.
2
We compared the drag coefficient with experimental
data by Fage (1936) and McNown and Newlin (1951),
resulting in a pretty good agreement, see Figure 3. It is
worth noting that care must be taken regarding wall ef
fects. This can be clearly seen in Figure 3, where curves
for different values of the ratio A d/1 are plotted.
Here, d is the particle's diameter and I the diameter of
u(:,i)= v(:,i)
the fluid domain. For low Reynolds numbers Re = ,
where large wall friction is present, a large ratio A pro
duces results that are far off the desired ones.
250
Numerical results x=0.1
200 Experimental data =0.0
Experimental data =0.2
SNumerical results x=0.2
150
100
Q 0
Figure 3: Comparison of drag curves.
Conclusion
In this paper a finite element method to compute par
ticular flow was presented. The methods is based on a
weak formulation in the form of a one domain approach,
a splitting scheme in time, adaptive grids and a subspace
projection method.
A validation with a sedimenting particle results in
good agreement to experimental results.
Acknowledgement
This work was supported by the Bayerische Forschungss
., ni,. which is gratefully acknowledged.
References
Eberhard Binsch. Local mesh refinement in 2 and 3 dimen
sions. IMPACT Comput. Sci. Engrg, 3(3):181191, 1991. doi:
10.1016/08998248(91)90006G.
Eberhard Binsch. Simulation of instationary, incompressible
flows. Acta Math. Univ. Comenianae, 67(1):101114, 1998.
S. Bonisch and V. Heuveline. On the numerical simulation of
the instationary free fall of a solid in a fluid. i. the newtonian
case. Computers and Fluids, 36:14341445, 2007.
A. Fage. Experiments on a sphere at critical reynolds numbers.
Aeronautical Research Committee Reports and Memoranda,
(1766):120, 1936.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
J. Feng, H.H. Hu, and D.D. Joseph. Direct simulation of initial
value problems for the motion of solid bodies in a newtonian
fluid part 1. sedimentation. J. Fluid Mech., 261:95134, 1994.
R. Glowinski, T.W. Pan, T.I. Hesla, and D.D. Joseph. A dis
tributed lagrange multiplier/fictitious domain method for par
ticulate flows. Int. J. Multiphase Flow, 25:755794, 1999. doi:
10.1016/S03019322(98)000482.
J.L. Guermond and Jie Shen. On the error estimates for the
rotational pressurecorrection projection methods. Math. of
Comp., 73:17191737, 2004.
H.H. Hu. Direct simulation of flows of solidliquid mix
tures. International Journal of Multiphase Flow, 22(2):335
352, 1996.
N.S. Martys and R.D. Mountain. Velocity Verlet algorithm for
dissipativeparticledynamicsbased models of suspensions.
Physical Review E, 59(3):37333736, 1999.
John McNown and John Newlin. Drag of spheres within cylin
drical boundaries. US National Congress ofApplied Mechan
ics, pages 801806, 1951.
Rodolphe Prignitz and Eberhard Binsch. Simulation of sus
pension induced rheology. To be published in Kybernetika,
2010.
Decheng Wan and Stefan Turek. Fictitious boundary and mov
ing mesh methods for the numerical simulation of rigid par
ticulate flows. Journal of Computational Physics, 222:2856,
2007.
