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 Erlangen-Nimrberg, Haberstr. 2, 91058 Erlangen, Germany
prignitz @ am.uni-erlangen.de and baensch @ am.uni-erlangen.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 finite-element discretisation
in space and an operator-splitting 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 3d-case 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 Navier-Stokes-equations. 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 non-dimensional
form reads


tu + u Vu + V pi D[u]) f in Q (t),
(Pff 1 Re
(1)
V-u 0 in (t), (2)
u 0 on FD, (3)
u = U + x r on OP(t), (4)


MU =F+j a-nds,


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+u-Vu-v acc+ Vpv- 1 V-D[u]-vbc
n(t) M.n(t) Re
+ MU V + Id; ( + w x (Iiw) -


12(t)


f .vcb+ F V +T


+ o. n-V+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 Vp-v- V-D[u]-vc-

pV v 1 D[u] : D[v] c
J-n(t) 2Re

SDpn-eD[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)}


.-n-V+rxa-n-. 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, )
S-2 (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)+F-V, (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 Navier-Stokes
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 Navier-Stokes 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 non-standard 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 pseudo-code 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):181-191, 1991. doi:
10.1016/0899-8248(91)90006-G.

Eberhard Binsch. Simulation of instationary, incompressible
flows. Acta Math. Univ. Comenianae, 67(1):101-114, 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:1434-1445, 2007.

A. Fage. Experiments on a sphere at critical reynolds numbers.
Aeronautical Research Committee Reports and Memoranda,
(1766):1-20, 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:95-134, 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:755-794, 1999. doi:
10.1016/S0301-9322(98)00048-2.

J.L. Guermond and Jie Shen. On the error estimates for the
rotational pressure-correction projection methods. Math. of
Comp., 73:1719-1737, 2004.

H.H. Hu. Direct simulation of flows of solid-liquid mix-
tures. International Journal of Multiphase Flow, 22(2):335-
352, 1996.

N.S. Martys and R.D. Mountain. Velocity Verlet algorithm for
dissipative-particle-dynamics-based models of suspensions.
Physical Review E, 59(3):3733-3736, 1999.

John McNown and John Newlin. Drag of spheres within cylin-
drical boundaries. US National Congress ofApplied Mechan-
ics, pages 801-806, 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:28-56,
2007.