A Coupled Lattice Boltzmann Method and Discrete Element Method for
Discrete Particle Simulations of Particulate Flows
C. Rettingera,∗, U. Rüdea,b

arXiv:1711.00336v1 [] 1 Nov 2017

a Chair

for System Simulation, Friedrich–Alexander–Universität Erlangen–Nürnberg, Cauerstraße 11, 91058 Erlangen,
Germany
b CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse Cedex 1, France

Abstract
Discrete particle simulations are widely used to study large–scale particulate flows in complex geometries
where particle–particle and particle–fluid interactions require an adequate representation but the computational cost has to be kept low. In this work, we present a novel coupling approach for such simulations.
A lattice Boltzmann formulation of the generalized Navier–Stokes equations is used to describe the fluid
motion. This promises efficient simulations suitable for high performance computing and, since volume displacement effects by the solid phase are considered, our approach is also applicable to non–dilute particulate
systems. The discrete element method is combined with an explicit evaluation of interparticle lubrication
forces to simulate the motion of individual submerged particles. Drag, pressure and added mass forces
determine the momentum transfer by fluid–particle interactions. A stable coupling algorithm is presented
and discussed in detail. We demonstrate the validity of our approach for dilute as well as dense systems by
predicting the settling velocity of spheres over a broad range of solid volume fractions in good agreement
with semi–empirical correlations. Additionally, the accuracy of particle–wall interactions in a viscous fluid
is thoroughly tested and established. Our approach can thus be readily used for various particulate systems
and can be extended straightforward to e.g. non–spherical particles.
Keywords: particulate flow, discrete particle simulation, fluid–particle coupling, lattice Boltzmann
method, discrete element method

1. Introduction
In the broad field of particulate flows, simulations are becoming increasingly important to study their
complex behavior and to enhance our understanding of the underlying processes. Ideally, one would like to
apply so–called direct numerical simulations (DNS) that are characterized by geometrically fully resolved
particles and a fine enough numerical resolution to capture all scales of fluid motion. They have been applied
successfully to e.g. fluidized beds [1–3] and sediment beds in rivers [4–7]. As a main benefit, they enable a
detailed analysis of flow properties inside closely packed particle clusters and naturally allow for an accurate
evaluation of trajectories and interaction forces of individual particles. Such quantities are usually difficult
to obtain in laboratory experiments. However, the large number of particles typically encountered in such
systems poses a major challenge for DNS and restricts their applicability to relatively small system sizes
with currently at most O(106 ) particles [1, 5, 6].
Alternative simulation approaches exist that introduce certain degrees of modeling to account for then
unresolved physical processes. This decreases the computational cost significantly and enables simulations
of larger systems over a longer period of time. Those approaches can be roughly classified as either two–fluid
models (TFM) or discrete particle simulation models (DPS) [8]. The TFM is an Eulerian–Eulerian model
∗ E-mail:

christoph.rettinger@fau.de

Preprint submitted to Elsevier

November 2, 2017

where the solid phase, like the fluid phase, is regarded as a continuum and its evolution is described by a
set of additional model equations. This approach requires that several modeling assumptions are made and
naturally can not represent the details of particle–particle or particle–fluid interactions adequately [8]. In
DPS models, the solid phase is represented by the actual particles, using a Lagrangian description of their
motion. Particle collisions are usually treated with the discrete element method (DEM) [9]. Different from
a DNS, however, the particle size is typically of the order or smaller than the numerical grid spacing used
for the representation of the fluid phase. The fluid–particle interaction, as well as possible turbulent flow
structures, are thus no longer resolved by the numerical simulation and have to be modeled. The overall idea
of these approaches is to establish the coupling of the fluid and the solid phase via interaction forces. For each
particle, such a fluid–particle interaction force is evaluated and applied, and a corresponding reaction force
is employed on the surrounding fluid. Therefore, closure correlations for most of the different contributions
to this interaction force are required. In most applications, the drag force is the most important contribution
and several correlations have been proposed here [10, 11].
DPS is a well–established method in classical CFD, where it is often referred to as CFD–DEM. It is also
the basis of several commercial tools and frequently used in various application fields, see the reviews in
[12–14]. In recent years, the lattice Boltzmann method (LBM) has become a prominent alternative to those
CFD methods [15]. Since it essentially requires only local operations to simulate fluid flow, it is well–suited
for highly parallel and scalable simulations to be executed an large supercomputers [16]. This feature makes
the LBM a viable choice for DNS [1, 2, 5, 7]. Furthermore, different approaches have been developed to
carry out DPS by coupling the LBM with the DEM [17–20]. We will briefly review these LBM–DEM DPS
approaches and summarize their properties and areas of applicability, see Sec. 2. As we will see, there is
currently no single approach available that has demonstrated to yield satisfactory results for the whole range
of dilute to dense particulate systems in terms of accuracy and stability.
The aim of this work is thus to develop a new coupled LBM–DEM scheme for DPS that is carefully
validated and can then be readily applied to different scenarios. It should retain the unique features of the
LBM that are important for the efficient parallel execution and allow for a flexible extension of the applied
models if required by a specific simulation. Additionally, measures for a stable momentum transfer between
the phases have to be considered. We will mainly focus on fluid–solid systems with spherical particles
which will guide our choice of interaction force contributions, namely lubrication and added mass forces,
and the empirical correlations that we include in our model. Nevertheless, the resulting approach will also
be applicable to gas–solid flows and can be extended to non–spherical particles.
The remainder of this work is structured as follows: After the review of existing LBM–DEM DPS
approaches in Sec. 2, we present our numerical model in detail in Sec. 3. This is split into discussions of
the treatment of the solid phase in Sec. 3.1 and the fluid phase in Sec. 3.2. The complete algorithm is
summarized in Sec. 3.3. We demonstrate the validity of our approach for the prediction of settling velocities
over a large range of solid volume fractions in Sec. 4. The accuracy of particle–wall interactions is analyzed in
Sec. 5. In Sec. 6, we summarize the approach and the results, and give an outlook to further developments.
2. Review of existing LBM–DEM DPS approaches
Several DPS approaches that couple the LBM and the DEM have been proposed. We briefly review
their main features and discuss their applicability. We subdivide the DPS approaches into two variants: the
two–way coupling and the four–way coupling approaches. They all have in common that interaction forces
between the particles and the fluid phase are evaluated, typically via interpolations to transfer information
from the Eulerian (fluid) to the Lagrangian (particle) description and vice versa. This force is then applied
onto the particles. A corresponding reaction force is employed on the surrounding fluid cells. In this
overview, we do not consider approaches that neglect the effect of the particles onto the fluid, referred to as
one–way coupling.
2.1. Two–way coupling approaches
In these approaches, the fluid flow is described via the incompressible Navier–Stokes equations and the
interaction between the two phases is established solely via interaction forces. In the context of the LBM,
2

this was explored in a pioneering paper [17] for the coupling with molecular dynamics simulations. Here,
the interaction force F f →p = −ζ(up − uf (xp , t)), with a friction coefficient ζ, is evaluated at the polymer
position xp and requires the interpolated fluid velocity uf as well as the polymer velocity up . An overview
of this method as well as further developments can be found in [21].
Another method is the so–called Particulate Immersed Boundary Method (PIBM) from [18] that can be
interpreted as a limiting case of the immersed boundary method of [22]. The interaction force is calculated
by interpolating the particle distribution functions available in the LBM, see Sec. 3.2.2, to the particle
position and applying the momentum exchange formulation in LBM. It can be shown that this interaction
force density takes the form f f →p,PIBM = 2ρf Ap (uf (xp , t) − up ), see e.g. [23], with Ap the cross sectional
area of the particle. Thus, the PIBM is very similar to the aforementioned coupling approach of [17]. It
has certain inconsistencies in terms of the physical units of f f →p,PIBM , and it lacks the contribution from
the fluid viscosity in comparison to e.g. Stokes’ drag force. This might explain the undesired viscosity
dependency of the single sphere settling velocity, as observed in [18]. Recently, the PIBM approach has
been enhanced by modifying the PIBM force with empirical drag correlations [24], which has been further
investigated by [25]. Both studies confirm that the enhanced PIBM is capable of predicting the settling
velocity of a single sphere correctly, but a thorough investigation of its accuracy for non–dilute systems is
still missing.
2.2. Four–way coupling approaches
As has been investigated in [26], it is crucial to account for fluid volume displacement effects when
considering non–dilute systems, or particles with sizes comparable to the size of the control volumes, i.e. grid
cells, to obtain accurate results. Such approaches are often referred to as four–way or volumetric coupling
approaches since they typically include the local fluid volume fraction εf into the Navier Stokes equations.
A commonly used variant are the volume–averaged Navier Stokes (VANS) equations [27], especially in the
area of gas–fluidized beds [8, 28]. A LBM formulation of the VANS equations exists [29, 30] and has been
combined with empirical drag correlations to carry out simulations of gas–solid fluidized beds [19]. However,
it has been shown that this formulation yields spurious velocities in the presence of larger solid volume
fraction gradients which could render the simulations inaccurate and unstable in heterogeneous systems
[31]. The proposed new scheme in [31] overcomes these problems but introduces several correction terms
that in turn render the LBM inefficient and tedious. An alternative LBM formulation of the VANS equations
exists [32] which, however, necessitates the inclusion of mass sources and a local solution of matrix systems,
which again adds substantial complexity to the LBM.
Another variant of such a coupling approach was presented in [20], which uses the partially saturated cells
method (PSM) from [33] to include the effect of fluid volume displacement. This method is usually applied
in fully resolved simulations, as e.g. in [34] where particles span several grid cells. Since the particles in
DPS are smaller than a grid cell, the computed forces from PSM are expected to be inaccurate and are thus
rescaled in [20] with the help of empirical drag laws. This approach has been applied successfully to gas–solid
fluidization simulations. However, a thorough analysis of the PSM is not yet available and it thus remains
unclear which macroscopic equations are solved with the approach. Additionally, the applicability of this
approach for dilute systems is questionable since the proposed rescaling relies on averaging the interaction
force over all particles inside a computational cell and thus the accuracy might deteriorate if only a few
particles are present in a cell.
3. Numerical Method
In this section, we present the equations describing the particle motion, Sec. 3.1, and the fluid dynamics,
Sec. 3.2. We discuss the applied numerical methods and the included models to account for various physical effects. All these components are implemented in the open–source LBM framework waLBerla1 . A
summary of the complete algorithm is given in Sec. 3.3.
1 http://walberla.net

3

3.1. Solid phase
The motion of a particle, given by the temporal evolution of its position xp , linear velocity up and
angular velocity ω p , is described by the following set of equations
dxp
= up ,
dt
dup
ρp Vp
= F c + F f →p + F lub + F g ,
dt
dω p
Ip
= T c,
dt

(1a)
(1b)
(1c)

with the particle density ρp , volume Vp and moment of inertia Ip . The net force acting on a particle consists
of different contributions. These are: collision forces F c , fluid–particle interaction forces F f →p , lubrication
forces F lub , and gravity F g . The fluid–particle interaction force is the sum of all interphase momentum
transfer contributions which are discussed in Sec. 3.1.2. The time integration of above equations is carried
out by an explicit time–stepping scheme [35] with a time step ∆tDEM which is efficient and accurate enough
for the regarded systems [8].
3.1.1. Discrete element method
To compute the forces and torques arising from particle–particle or particle–wall collisions, the discrete
element method (DEM) from [9] is used. Collisions are modeled based on a linear mass–spring–damper
system and thus result in a so–called soft–sphere model. In this work, we will restrict ourselves to spherical
particles for simplicity, noting however that an extension to more complex shapes is possible [36]. When
considering two particles, denoted by indices i and j, inside a system of Np particles, we can define an
overlap length
(2)
δij = 21 (dp,i + dp,j ) − |xp,j − xp,i |,
with the particle diameter dp and its center position xp , see the sketch in Fig. 1. The contact normal unit
vector of the collision is defined as
xp,j − xp,i
,
(3)
nij =
|xp,j − xp,i |
and the relative velocity at the contact point C is
ur,ij = up,i − up,j + ( 12 dp,i ω p,i + 21 dp,j ω p,j ) × nij .

(4)

This velocity can be decomposed into the normal and the tangential part, given as
unr,ij = (ur,ij · nij )nij ,

(5a)

utr,ij

(5b)

= ur,ij −

unr,ij ,

respectively.
Following the spring–damper collision model, with a stiffness constant kn and a normal damping constant
ξn , the normal collision force on particle i due to this contact is
(
−kn δij nij − ξn unr,ij , δij > 0,
n
F c,ij =
(6)
0,
else.
The collision force in tangential direction distinguishes between a sliding and a sticking regime, and is given
as
(
−min(µc |F nc,ij |, ξt |utr,ij |)tij , if |utr,ij | > 0,
t
F c,ij =
(7)
0,
else,

4

ω p,j

δij
up,j
ω p,i

dp,j

dp,i
C
xp,i

xp,j
nij
tij
up,i

Figure 1: Schematic representation of two colliding spheres, i and j, together with reference quantities required by the DEM
collision model.

with the tangential unit vector tij = utr,ij /|utr,ij |, the friction coefficient µc , and the tangential damping
constant ξt . Consequently, the total force on particle i due to possible collisions with all other particles in
the system is then given as
Np
X
(F nc,ij + F tc,ij ).
(8)
Fc =
j=1,i6=j

The tangential collision force also generates a torque on the particle and thus the total collisional torque
on particle i is
Np
X
dp,i t
Tc =
F
× nij .
(9)
2 c,ij
j=1,i6=j

In total, this DEM model introduces four parameters, namely kn , ξn , ξt and µc . In our simulations, we
always choose ξn = ξt . Additionally, we use
p
−2 ln en mij kn
ξn = p
,
(10)
π 2 + ln2 en
which is an analytical solution of the linear spring–damper problem [37], and relates ξn to kn . It features
the dry coefficient of restitution, en , that is often available in literature for different materials, and the
reduced mass mij = mp,i mp,j /(mp,i + mp,j ) for a particle pair and mij = mp,i for a particle–wall contact.
Furthermore, the duration of a collision event can be estimated as [37]:
Tc = p

2πmij
.
4mij kn − ξn2

(11)

By choosing a specific Tc , and together with Eq. (10), the stiffness and damping coefficients can be computed.
3.1.2. Fluid–particle interaction forces
The hydrodynamic force acting on the individual particles contains several contributions that account
for different physical effects. Specifically, we here include drag F d , pressure F pr , lift F l , and added mass
F am effects:
F f →p = F d + F pr + F l + F am
(12)
The choice for or against the inclusion of certain effects depends on the targeted setup. In simulations of
gas–solid systems, such as gas–fluidized beds, often only the first two components are modeled. In fluid–solid
systems, like sediment transport in river beds, the lift and added mass contributions cannot be neglected
5

force

formula

drag force

F d = 3πdp µ0 (1 − εp|p )Cd (Rep|p , εp|p )(uf |p − up ),


Cd (Rep ,0)
Cd (Rep , εp ) = (1 − εp ) (1−ε
+
A(ε
)
+
B(Re
,
ε
)
3
p
p
p
p)
Cd (Rep , 0) = 1 + 0.15Re0.687
p
5.81εp
(1−εp )3

ε1/3

p
+ 0.48 (1−ε
4
p)


0.61ε3
B(Rep , εp ) = ε3p Rep 0.95 + (1−εp p)2

A(εp ) =
pressure gradient force
lift force
added (virtual) mass force

F pr = −Vp (∇P )|p
q µρ

0 f
(uf |p − up ) × (∇ × uf )|p
F l = 1.61d2p |∇×u
f ||p



Duf
dup
F am = Cam ρf Vp
−
, Cam = 0.5
Dt
dt
|p

Table 1: Considered contributions to the fluid–particle interaction force and the applied (closure) relations.

xp

∆x
∆x
Figure 2: Schematic representation of the interpolation operator in Eq. (14). The gray circle depicts a particle at position xp
on an Eulerian grid, with cell centers x marked by black squares. The blue circle around the particle is of radius 1.5∆x and
thus shows the extent of the filter kernel. All cells with cell centers inside this circle contribute to the interpolation, displayed
by the black lines.

and are thus typically included [38, 39]. The models for the different contributions are displayed in Tab. 1
and are briefly discussed next.
In all cases, different Eulerian field variables must be evaluated at the particle position, which will be
denoted as φ|p . This requires an interpolation from the cell centers x of the Eulerian grid cells to the particle
position xp and is given as an interpolation operator I (x, xp ), such that:
X
φ|p =
I (x, xp )φ
(13)
x

In this work, we use the three point discrete delta functions from [40] that we will also need for the opposite
operation later in Sec. 3.2.1:

√
1
2

 3 (1

 + −3rp+ 1,
 r < 0.5,
1
I (x, xp ) = 6 5 − 3r − −3(1 − r)2 + 1 , 0.5 ≤ r < 1.5,
(14)


0,
1.5 ≤ r,
with r = |x − xp |/∆x and the Eulerian grid spacing ∆x, see also Fig. 2. Other possibilities are trilinear
interpolation [38], Lagrange polynomial interpolation [18], or other filter kernels [41]. Note, that apart from
the fluid quantities like velocity and pressure, we also interpolate the solid volume fraction εp .
6

The drag force acting on a single sphere in the Stokes flow regime is given by Stokes’ law as F d,St =
3πdp µ0 (uf |p − up ). In order to extend the applicability to higher Reynolds number flows, it is scaled by
a drag coefficient Cd (Rep , εp ). Various such drag correlations are available [10, 11]. The particle Reynolds
number based on the superficial fluid velocity is
Rep =

εf |p ρf dp |up − uf |p |
,
µ0

(15)

with the local fluid volume fraction εf |p = 1 − εp|p , the fluid density ρf and the dynamic viscosity µ0 . Here,
we make use of the drag correlation from [42] since it is derived as an extension to the single sphere case,
Cd (Rep , 0), and is valid for the regarded Reynolds numbers.
The pressure force, which features the gradient of the total fluid pressure, allows the particle to respond
to local changes in the dynamic pressure. It also includes the buoyancy term if gravity is acting on the
fluid in the simulation. Otherwise, the buoyancy force F buoy = −ρf Vp g must be added to the particles in
addition to the gravitational force F g = ρp Vp g, with the gravitational acceleration g.
The added, or virtual, mass force models the observation that an accelerating particle experiences a force
against its acceleration since the fluid surrounding the particle has to be accelerated as well [37, 43]. It thus
appears as if the particle has gained, or added, mass. Since in our simulations, the densities of the solid and
the fluid phase are of the same order, we expect that the added mass force may have a non–negligible effect.
At last, also shear forces act on the particle and we account for them by including a shear–induced lift
force [44].
3.1.3. Lubrication force
When two immersed objects approach each other, the fluid between these objects gets squeezed out
of the gap and this flow can produce a significant resistance to such a (wet) collision. As a result, the
rebound velocity after the collision is suppressed and a lower effective coefficient of restitution is observed.
Simulation models, that do not accurately resolve the flow inside these gaps, thus have to take care of the
unresolved lubrication forces to obtain the correct rebound behavior. In this work, we explicitly calculate
the lubrication force F lub and add it to the particle, Eq. (1). The lubrication force acting on a sphere i due
to another sphere j can be computed as [45]:

F lub,i = −asq (upr,ij · nij )nij − ash

2
hij + hdi

2
[upr,ij − (upr,ij · nij )nij ],

with nij from Eq. (3), the relative particle velocity upr,ij = up,i − up,j , a mean diameter hdi =
and the gap size hij , calculated as
hij = |xp,j − xp,i | − 21 (dp,i + dp,j ).

(16)
2dp,i dp,j
dp,i +dp,j ,

(17)

The coefficients are given as
asq
ash





hdi
18
hdi
9 hij
hdi
+
ln
+
ln
,
4hij
40
2hij
84 hdi
2hij


1
hdi (hdi + hij )2
= πµ0 hdi ln
.
2
2hij
4

3
= πµ0 hdi
2



(18a)
(18b)

We currently do not include contributions from angular particle velocities and lubrication effects on the
torque, which are also given in [45]. For a sphere–wall interaction, we assume dp,j → ∞ and obtain
hdi = 2dp,i . We consider those forces only if the gap size is below a certain cut off value hc . Such a
procedure is also commonly applied in fully resolved simulations [46]. Additionally, we pose a lower limit
on the gap size, such that hij := max(hij , 10−5 ∆x), to avoid too large forces [47]. The accuracy of this
approach is demonstrated in Sec. 5.
7

3.2. Fluid phase
Next, we present and discuss the models of the fluid flow where we account for the fluid volume displaced
by the solid phase. The corresponding LBM formulation is presented, as well as the evaluation of the fluid–
particle interaction forces.
3.2.1. Generalized Navier–Stokes equations
The fluid flow inside a porous medium is commonly described via the so–called generalized Navier–Stokes
(GNS) equations [48, 49]:

ρf

∂ ūf
+ (ūf · ∇)
∂t



∇ · ūf = 0,

ūf
= −∇(εf p) + µe ∇2 ūf + f p→f + εf f b ,
εf

(19a)
(19b)

where εf is the fluid volume fraction that can vary locally, ρf is the fluid density, p the fluid pressure,
ūf = εf uf is the volume–averaged fluid velocity, and µe is the effective fluid viscosity. For εf = 1,
the incompressible Navier–Stokes equations are recovered. The forcing acting on the fluid comprises two
parts: The body force density f b is induced by an external force. The other part, f p→f , accounts for the
interphase momentum transfer due to the presence of the porous medium. It usually consists of the Darcy
and Forchheimer drag correlations which are the linear and nonlinear parts of Ergun’s experimental drag
correlation [50]. In order to account for the velocity of the particle phase and to incorporate other important
contributions, we enhance and modify this interphase momentum transfer term.
To define those quantities, a filter operation D must be introduced first that allows to distribute a particle
property φp to the surrounding fluid cells in order to obtain the continuous field
φf (x) =

Np
X

D(x, xp,i )φp,i .

(20)

i=1

In this work, we choose D to be the same as the previously defined interpolation filter I from Eq. (14).
With the help of Eq. (20), the fluid volume fraction in a cell can now be calculated from the particle
positions via


Np
1 X
D(x, xp,i )Vp,i  .
(21)
εf (x) = 1 − εp (x) = 1 −
(∆x)3 i=1
In the same way, and with the modifications from Sec. 3.1.2, the force density accounting for the interphase
momentum transfer is given as
f p→f = −

Np
1 X
D(x, xp )(F d + F l + F am )i .
(∆x)3 i=1

(22)

Note, that in contrast to the hydrodynamic force on the particle, Eq. (12), we purposely do no include the
pressure force F pr here since this contribution is already present in the pressure term of the momentum
GNS equation, Eq. (19b).
Following [38], a modified version of Eilers equation is used to include the effect of the solid phase on
the dynamic viscosity µ0 [51]:
2

0.5[µ]εp
,
(23)
µ∗ = µ0 1 +
(1 − εp )/εcp
with the intrinsic viscosity [µ] = 2.5, that models the effect of the particle shape on the rheology of the
fluid–particle system, and the solid volume fraction for a close sphere packing, εcp = 0.64. Additionally,
a large–eddy (LES) approach is applied to model the unresolved subgrid–scale turbulence which, for the
simplest case of a Smagorinsky type turbulence model, adds a turbulent contribution µt to Eq. (23) and will
be discussed in more detail in the next section. Thus, the effective viscosity µe is given as
µe = µ∗ + µt .
8

(24)

3.2.2. Lattice Boltzmann formulation
For the simulation of the incompressible Navier–Stokes equations, we employ the lattice Boltzmann
method [15]. It computes the evolution of particle distribution functions (PDFs) on a Cartesian lattice by
solving the discretized version of the Boltzmann equation. Each of these PDFs fq corresponds to a discrete
lattice velocity cq . A common choice is the D3Q19 lattice model [52] that features a set of 19 lattice
velocities, i.e. q ∈ {0, . . . , 18}. A lattice Boltzmann model for the generalized Navier–Stokes equations,
Eqs. (19), has been proposed by [49] and is given as

eq
fq (x + cq ∆t, t + ∆t) = fq (x, t) + ∆t
(25)
τe fq (ρf , ūf , εf ) − fq (x, t) + Fq (x, ūf , εf , t)∆t.
It uses a single relaxation time parameter τe ∈ ( 12 , ∞) to linearly relax the PDFs towards their equilibrium
values fqeq , which can be computed via


(cq · ūf )2
ūf · ūf
cq · ūf
eq
+
−
.
(26)
fq (ρf , ūf , εf ) = wq ρf 1 +
c2s
2εf c4s
2εf c2s
The lattice weights wq are as given e.g. in [52] and cs is the lattice speed of sound. The forcing operator in
Eq. (25) is used to incorporate external forces and can be written as [49, 53]:



∆t
ūf
cq · ūf
cq
Fq (x, ūf , εf , t) = wq 1 −
−
+
cq · f ext (x, t),
(27)
2τe
c2s
εf c2s
εf c4s
with a force density f ext = f p→f + εf f b . The fluid density ρf and the volume–averaged fluid velocity ūf
are cell local quantities and calculated via moments of the PDFs:
X
ρf (x, t) =
fq (x, t),
(28)
q

∆t ext
1 X
fq (x, t)cq +
f (x, t).
ūf (x, t) =
ρf q
2ρf

(29)

The pressure is given as p = c2s ρf /εf . The effective dynamic viscosity µe is related to the relaxation time τe
via
2
µe /ρf = (τe − ∆t
(30)
2 )cs .
Analogously to Eq. (24), this relaxation time includes two contributions,
τe = τ∗ + τt ,

(31)

where τ∗ and τt are the relaxation times that correspond to the modified fluid viscosity µ∗ from Eq. (23)
and the turbulence viscosity µt , respectively [54, 55]:
µ∗ /ρf = ν∗ = (τ∗ −
µt /ρf = νt =

∆t 2
2 )cs ,

(32a)

τt c2s

(32b)

For the Smagorinsky type LES turbulence model, the turbulence viscosity is obtained from the strain rate
tensor S = Sα,β = 21 ((∇ūf ) + (∇ūf )> ) as [54]:
s X
Q̄
2
νt = (Cs ∆xLES ) S̄, S̄ =
Q̄α,β Q̄α,β ,
(33)
, Q̄ = 2
2ρf c2s τe
α,β

with the Smagorinsky constant Cs = 0.1 [55], a filter length ∆xLES = ∆x, and the (filtered) mean momentum
flux Q̄. The latter one is computed from the momentum fluxes Qα,β , obtained as the second–order moments
of the nonequilibrium parts of the distribution functions [49, 54]:
X
Qα,β =
cq,α cq,β [fq − fqeq ].
(34)
q

9

This leads to a quadratic equation for τt = νt /c2s , with the solution given as [54]
q

√
1
2
2
4
−1
τt =
τ∗ + 2 2(CS ∆xLES ) (ρf cs ) Q̄ − τ∗ .
2

(35)

In the context of the LBM, all quantities are commonly expressed in terms of lattice units which results
in ∆t = 1, ∆x = 1, cs = √13 , and ρf ≈ 1.
3.2.3. Fluid–particle interaction forces
The components of the fluid–particle interaction force density f p→f are given in Eq. (22) and are included via f ext in Eq. (27). Their definitions are as discussed in Sec. 3.1.2. Their evaluation requires
macroscopic fluid quantities and their derivatives instead of the readily available LBM quantities, i.e. the
particle distribution functions. We therefore briefly discuss the LBM specific aspects of the force evaluation
here.
In general, we compute and store all required quantities explicitly for each computational cell in the
simulation domain. Those can then be interpolated to the particle position via Eq. (13). The velocity used
to compute the interaction forces is the fluid–phase velocity uf rather than volume–averaged fluid velocity.
To evaluate the appearing derivatives, we use the lattice differential operators from [56]. Consequently, the
gradient of a scalar quantity φ, as e.g. P = εf p in the pressure gradient force, can be approximated as
∇φ ≈

1 X
wq cq φ(x + cq ∆t).
c2s q

(36)

The rotation of a vector field ψ, as e.g. the velocity field in the lift force, is computed as
∇×ψ ≈

1 X
wq cq × ψ(x + cq ∆t).
c2s q

(37)

A factor that significantly influences the stability of this specific fluid–particle coupling approach originates from the definition of the macroscopic velocity in Eq. (29). When the external force is changed, e.g.
by adding the interaction force, the macroscopic fluid velocity changes also. As a consequence, a reevaluation of the interaction force components, that depend on the velocity, with the updated velocity values
will yield distinct force values. If not handled correctly, this can easily lead to severe oscillations of the
force term which then affect the overall stability of the simulation. A cell–local implicit treatment of this
dependency for each component of the interaction force would be desirable, as e.g. done in [49] for Ergun’s
drag correlation. However, this is not applicable here since in general the applied (drag) correlations are
more complex than Ergun’s correlation and the influence of all present particles would have to be treated
implicitly at the same time. We therefore rely on an iterative reevaluation of the required fluid quantities
with the updated force contributions. This procedure works well for e.g. the drag force which is usually the
dominant contribution. In our simulations, we observe that small oscillations in the drag force in combination with the added mass force, which depends on the temporal change of the fluid velocity, are often the
source of growing instabilities, especially for large solid volume fractions. In order to damp this feedback
Du
du
mechanism, we refrain from reevaluating Dtf , see Tab. 1, during those internal iterations. Since dtp is
updated in each iteration, the added mass force still changes. For additional stability, especially in cases of
large solid volume fractions, we reevaluate the fluid velocity after the contributions from the added mass
and lift forces have been computed and before the drag force is evaluated.
3.3. Complete DPS algorithm for coupled LBM–DEM simulations
A detailed overview of our proposed algorithm to carry out discrete particle simulations by coupling the
lattice Boltzmann method and the discrete element method is presented in Alg. 1. It features an outer time
loop with time step size ∆t which is also the time step used in the LBM simulations. However, the typical
time scales of particle collisions are usually significantly smaller than the time scales of the fluid solver.
10

Initially, evaluate and store fluid volume fraction εf in each cell, Eq. (21);
Apply external body force, εf f b , to cells;
for each time step do
Du
Evaluate and store Dtf ;
for each interaction subcycle do
Evaluate and store uf , ∇P , ∇ × uf ;
for each particle do
Evaluate F l and distribute to cells;
Evaluate F am and distribute to cells;
end
Reevaluate uf ;
for each particle do
Evaluate F pr ;
Evaluate F d and distribute to cells;
Apply gravitational force F g ;
end
for each DEM substep do
Compute lubrication forces F lub on all particles, Eq. (16);
Compute collision forces F c , Eq. (8), and torques, T c , Eq. (9);
Perform time integration of particles, Eqs. (1), with time step size ∆tDEM ;
end
Reevaluate fluid volume fraction εf ;
end
Compute τ∗ , Eq. (32a), from modified viscosity µ∗ , Eq. (23);
Compute τt , Eq. (35), and then τe , Eq. (31);
Perform GNS–LBM step, Eq. (25), with time step size ∆t;
end
Algorithm 1: Overview of the LBM–DEM DPS algorithm.
Therefore, the time loop is subdivided twice, such that during one time step nint interaction subcycles are
performed which include nDEM DEM substeps. The ratio between the LBM and the DEM time step is then
∆t/∆tDEM = nint nDEM and easily reaches O(100), see Sec. 5.
The interaction subcycles are introduced to improve the accuracy of the interaction force evaluations
since the particles’ position, needed for the interpolation of fluid quantities, as well as their velocity are
updated in each subcycle. Additionally, this reevaluation of the interaction forces helps to avoid oscillations
in the force components effectively, see the discussion in Sec. 3.2.3.
The additional DEM substeps are used to deal efficiently with simulation scenarios that require accurate
resolution of the collision and lubrication interactions between particles, as e.g. in Sec. 5, or flow through
densely packed particle beds. During those substeps, the currently acting forces on the particles, F f →p and
F g , are kept constant, and only collision and lubrication contributions are updated. Since the computationally expensive evaluation of the interaction forces is thus omitted, the usage of DEM substeps is more
efficient than increasing the number of interaction subcycles.
Finally, we note that no global data exchange is required in our algorithm at any time since all parts
operate only on locally available data. This key feature ensures its efficient applicability to large–scale
simulations on HPC systems.
4. Hindered settling behavior of spheres
In order to reliably capture the fluid–solid interaction in various applications, the DPS approach must
predict the correct behavior over a wide range of solid volume fractions, i.e. from single particles to densely
11

packed scenarios. A suitable test setup to validate these cases is the evaluation of the hindered settling
behavior of spheres at various solid volume fractions. It is well known that the average relative settling
velocity of the spheres hupS i in an unbounded domain is affected by the average solid volume fraction hεp i
as follows [57]:
hupS i
= (1 − hεp i)κ ,
(38)
upT
where upT is the terminal settling velocity of a single sphere in an infinite fluid, which can be measured in
experiments or found by direct numerical simulations. The coefficient κ is a function of the Reynolds number
ReT = ρf upT dp /µ0 and its value for unbounded flows can be estimated by the widely used semi–empirical
correlations from Richardson and Zaki [57]:

4.65,



4.35 Re−0.03 ,
T
κ=
−0.1

4.45
Re
,

T


2.39,

ReT < 0.2
0.2 ≤ ReT < 1
1 ≤ ReT < 500
500 ≤ ReT

(39)

4.1. Setup
We use a fully periodic domain of size [Lx × Ly × Lz ]/dp = 32 × 32 × 32. Initially, spheres with a uniform
diameter dp are randomly generated inside this domain until the desired average solid volume fraction is
obtained. In order to have non–overlapping spheres before starting the coupled simulations, several DEM
time steps are carried out to resolve possible overlaps. Then, their velocity is set to zero. These spheres at
rest then begin to settle under gravity in the viscous fluid. Regarding the physical parameters, we follow
the experimental setup from [58] which is also used in simulations of [38]. Thus, the fluid phase is water
with ρf = 1000 kg/m3 and µ0 = 0.001 Pa s. The glass spheres have the properties dp = 0.35 · 10−3 m,
ρp = 2500 kg/m3 , en = 0.88, and µc = 0.25. Furthermore, we use Tc = ∆t/2, hc = dp , nint = 10, and
nDEM = 50. The gravitational acceleration is g = −9.81 m/s2 and acts only on the spheres. A uniform body
force density f b = (0, 0, −hεp ig(ρp − ρf ))> is applied on the fluid to balance the spheres’ weight. In the
experiments of [58], the terminal settling velocity of a single sphere upT = −0.048 m/s was measured, which
results in ReT = 16.8. The average relative sphere velocity in gravitational direction is evaluated as
hupS i =

Np
(∆x)3 X
1 X
up,i −
uf (x)
Np i=1
Lx Ly Lz x

(40)

The simulations are run for 30tSt time steps, with tSt = ρp d2p /18µ0 , which is long enough to obtain a
stationary average settling velocity [38]. The number of spheres ranges from 1 to 37549, corresponding to
hεp i ≈ 0 up to 0.6. We use three different sphere diameters, dp /∆x = 14 , 12 , and 1, to investigate resolution
effects.
4.2. Results
A visualization of the outcome of this test case is shown in Fig. 3. As can be seen from the spheres’
velocities, clusters of similar velocity form during the simulation. The results of the study are plotted in
Fig. 4 with the average relative settling velocity over the average solid volume fraction. Apparently, the
simulated settling velocities match the predictions from the Richardson–Zaki formula very well over the
whole range of solid volume fractions. In particular, the settling behavior of a single sphere is captured well.
For the densely packed scenarios slightly larger settling velocities are obtained. Altering the ratio between
sphere diameter and cell size does not change the settling velocity.

12

Figure 3: Visualization of the setup for hεp i = 0.3 and dp = 0.5. The color indicates the normalized relative settling velocity
of the individual spheres, upS /upT . Gravity acts in vertical direction.

case
A0
A1
A2
A3
A4
A5

lift force
X
X
X
X
X

lubrication force
X
X
X
X
X

turbulence model
X
X
X
X
X

mod. viscosity
X
X
X
X
X

nint
10
10
10
10
10
1

hupS i/upT
0.288
0.297
0.288
0.288
0.289
0.288

Table 2: Simulation parameters for the sensitivity study and their influence on the settling velocity for hεp i = 0.3 and
dp /∆x = 12 . A0 is the baseline case.

4.3. Sensitivity analysis
Next, we focus on a specific scenario and investigate the influence of different model parameters on the
obtained settling velocity. With such a sensitivity analysis, the most relevant parameters can be identified.
This is essential to develop computationally efficient numerical approaches since contributions with small
influence can possibly be neglected. As a setup, we chose dp /∆x = 21 and hεp i = 0.3, as visualized in Fig. 3.
The parameters of the cases under comparison are as summarized in Tab. 2. We evaluate the contributions of
the lift force, lubrication forces, the turbulence model, the modified viscosity of Eq. (23), and the number of
interaction subcycles. Without the lift force (A1), the settling velocity increases by around 3% in comparison
to the baseline case A0. Lubrication forces are seemingly not important in this rather homogeneous setup
since the settling velocity remains the same when neglected (A2). Also the influence of the turbulence
model is negligible here as the Reynolds number is not large enough such that turbulent flow structures
would become dominant (A3). The modifications of the fluid viscosity introduced by Eq. (23) also do not
alter the settling velocity significantly (A4). Finally, decreasing the number of interaction subcycles yields
the same settling velocity (A5), but the overlaps between the spheres become around four times larger since
also ∆tDEM is reduced.
4.4. Discussion
The results of this test case show that our approach is capable of correctly predicting the fluid–particle
interactions over a wide range of solid volume fractions, from a single sphere up to systems that are close
to the maximum packing fraction. The drag force, and thus the applied drag correlation, is the most
important factor in this setup but also the volume exclusion effects on the fluid phase become relevant for
larger solid volume fractions. Even though the sensitivity study showed no changes for a smaller number
13

hupS i/upT

100

10−1
Richardson–Zaki
dp /∆x = 1/4
dp /∆x = 1/2
dp /∆x = 1
10−2

0.0

0.1

0.2

0.3
hεp i

0.4

0.5

0.6

Figure 4: Average relative settling velocity normalized by the single sphere terminal settling velocity as a function of the average
solid volume fraction hεp i. For comparison, the Richardson–Zaki correlation [57], Eqs. (38) & (39), is included.

of interaction subcycles, scenarios with large solid volume fractions require more than one internal iteration
to yield stable results for the reasons mentioned in Sec. 3.2.3. We expect that the influence from using the
modified viscosity or the turbulence model will become more important in other setups like a shear flow
over a sediment bed at high Reynolds numbers which also features large velocity gradients. In such a case,
also lift forces will play a more prominent role on the sediment motion. The observation that the simulated
settling velocities are independent of the applied sphere diameter without tuning other numerical parameters
like the relaxation time, necessary e.g. in the PIBM approach from [18], again demonstrates the validity of
our approach. This property offers a great amount of flexibility as the ratio between diameter and cell size
can thus be chosen with respect to the requirements of the application setup. It also enables the accurate
simulation of polydisperse systems, i.e. mixtures of small and large particles. Furthermore, the approach
can be combined with grid refinement techniques where the particles move across meshes with different grid
sizes ∆x.
5. Sphere–wall collision in a viscous fluid
In this scenario, a sphere, heavier than the surrounding viscous fluid, is dropped in a box, and its
trajectory and velocity are tracked. When it hits the bottom plane, a rebound will happen which allows to
define an effective, or wet, coefficient of restitution
ewet
=−
n

upR
,
upT

(41)

where upR is the rebound velocity of the sphere and upT the terminal settling velocity. From experiments
[59, 60], it is known that the sphere motion, and thus the ratio between the wet and dry coefficient of
restitution, is mainly a function of the Stokes number,
St =

ρp ReT
,
ρf 9
14

(42)

1.0

1

80
60
40

0.5
0
−1 0

up /upT

(zp − dp /2)/dp

100

1

20
0
−100

−50
0
50
(t − tI )upT /dp

1
0

0.0 −1
−1
−0.5

0

1

−1.0
−30 −20 −10 0 10 20
(t − tI )upT /dp

100

30

Figure 5: Left: Wall–normal distance of the sphere surface to the wall over time, shifted such that the times of the first impact
coincide, for all simulated St and dp /∆x = 21 . Right: Normalized sphere velocity of St = 128 over time, again shifted. Inset
with marker of the rebound velocity evaluated at tI + tR .

with ReT = ρf upT dp /µ0 . With this setup, the accuracy and correctness of the collision and lubrication
forces as well as the substepping procedure will be tested.
5.1. Setup
In our simulations, we use a rectangular domain of size [Lx × Ly × Lz ]/dp = 32 × 32 × 512 and initialize a
resting sphere horizontally centered and close to the upper boundary. The domain is large enough such that
the terminal velocity is always reached before the sphere impacts on the bottom wall. For the fluid phase,
no–slip boundary conditions are applied in all directions. The gravity, only acting on the sphere, is given by
a gravitational acceleration g = −0.0001 in lattice units. The dry coefficient of restitution is en = 0.97 and
the coefficient of friction is µc = 0.1, which corresponds to a steel sphere and is similar to the one used in
the experiments of [60]. In order to keep the maximum penetration depth below 5% of the sphere diameter,
dp
we use Tc = ∆x
∆t, see Eq. (11). We keep the number of interaction subcycles at nint = 10 and choose
to
maintain a constant number of DEM steps per Tc . The cut–off gap size for the lubrication
nDEM = 25 ∆x
dp
forces is hc = dp . We define the impact time tI as the instance in time where the sphere loses contact with
the wall, i.e. the collision force becomes zero after the first collision. As in [61], various combinations of
density ratios and fluid viscosities are employed in order to obtain different Stokes numbers St ∈ [10, 2000].
Additionally, three different sphere diameters are tested, dp /∆x = 41 , 12 , and 1.
5.2. Results
In the left plot of Fig. 5, the trajectories of all simulations with dp /∆x = 21 are shown as wall–normal
distances to the bottom wall as a function of time, which is shifted by the impact time and normalized
with the reference time dp /upT . As can be seen from the overlapping, straight lines before the impact, all
spheres obtain their terminal settling velocity before any interaction with the wall happens. The zoomed
inset shows that in some cases the distance gets negative, indicating that overlaps between the sphere and
the wall happen which, however, are very small. The maximum rebound height ranges from negligible small
values for St below 17, which can be quantified as no rebound, up to around 48dp for the largest simulated
St. The right plot depicts the temporal evolution of the sphere velocity for St = 128. The inset of this plot
shows that already a few time steps before the actual collision with the wall, the sphere gets decelerated due
to lubrication forces. As marked in the inset, the rebound velocity is not taken as the maximum velocity
after the collision but evaluated some time tR = 0.1dp /upT after the impact. This has been proposed by
[61], to account for the temporal resolution of experimental measurement facilities used by [59] and [60], to
which the results are compared next.
The normalized wet coefficient of restitution as a function of St is evaluated according to Eq. (41)
and compared to experiments, all shown in Fig. 6. It can be seen that the simulation results fit to the
15

1.0

ewet
n /en

0.8
0.6
Joseph et al. (2001)
Gondret et al. (2002)

0.4

dp /∆x = 1/4
dp /∆x = 1/2

0.2

dp /∆x = 1
0.0

101

102
St

103

Figure 6: Ratio between the wet and dry coefficient of restitution as a function of the Stokes number St. For comparison,
experimental data of [59] and [60] are included.

case
B0
B1
B2
B3
B4
B5

nDEM
50
5
500
50
50
50

added mass
X
X
X
X
X

hc
dp
dp
dp
dp
0
dp /2

ewet
n /en
0.792
0.899
0.803
0.831
0.976
0.797

Table 3: Simulation parameters for the sensitivity study and their influence on the ratio between wet and dry coefficient of
restitution for St = 128 and dp /∆x = 12 . B0 is the baseline case.

experimental data very well over the whole range of St and also show an exponential increase. In particular,
the critical Stokes number, below which no rebound can be seen, agrees with the experimentally reported
ones. For St > 600, the wet coefficient of restitution approaches the dry one which is also in accordance
with the reported values of this so–called elastic limit. The overall behavior is independent of the particle
diameter used. Some perturbations can be observed, however, especially for dp /∆x = 1, which can be
explained by the low temporal resolution applied for evaluating the rebound velocity since tR has typically
a magnitude of a few simulation time steps.
5.3. Sensitivity analysis
We again carry out a sensitivity analysis for this setup to quantify the influence of different numerical and
model parameters. We here investigate the effect on the relative coefficient of restitution and the rebound
trajectory for the scenario with St = 128 and dp /∆x = 12 . These parameters are nDEM , and consequently
∆tDEM , the added mass force contribution, and the cut off distance hc in the lubrication forces, with values
given in Tab. 3. The resulting trajectories are shown in Fig. 7. When decreasing nDEM , and thus increasing
∆tDEM , the maximum rebound height as well as ewet
increases since the collision and lubrication interaction
n
is less accurately resolved. Increasing nDEM hardly changes the outcome in comparison to the baseline case,
indicating that nDEM = 50 is a good choice in this setup in terms of accuracy and computational efficiency.
Neglecting the added mass contributions to the interaction force (B3) slightly increases the rebound height
16

6
B0
B1
B2

(zp − dp /2)/dp

5

B3
B4
B5

4
3
2
1
0

−5

0

5

10
15
20
(t − tI )upT /dp

25

30

35

Figure 7: Wall–normal distance of the sphere surface to the wall over time for the setups of the sensitivity study, given in
Tab. 3, for St = 128.

and the coefficient of restitution. The largest influence has the lubrication force since the rebound behavior
is completely changed when no lubrication forces are considered (B4). The value of the wet coefficient of
restitution approaches the dry one and the rebound is significantly higher. As (B5) shows, the results vary
only slightly when choosing a smaller cut–off distance than in the baseline case.
5.4. Discussion
With the accuracy of the simulated settling velocity demonstrated in Sec. 4, this benchmark demonstrates
that also single collision events in a submerged setup are well predicted with our approach. In such a scenario,
a fine enough temporal resolution of the collision dynamics has to be ensured to prevent large overlaps
during the collision. Furthermore, depending on fluid flow properties, the inclusion of lubrication forces can
be crucial for a correct behavior. In our approach, where we explicitly compute and set the lubrication forces
for each sphere–sphere or sphere–wall pair, this evaluation can become the computationally most expensive
part of the simulation in systems with a large number of particles. Therefore, an other approach could be
used where the coefficient of restitution is changed dynamically for each collision, like in [38]. This requires
to compute a local impact Stokes number and to explicitly prescribe how the wet coefficient of restitution
changes depending on this number. However, such an approach introduces further modeling assumptions
and is unable to capture physical effects like the deceleration of the sphere already before the collision
happens. Besides the lubrication and collision forces, also added mass forces can influence the dynamics of
the system and should thus be included.
6. Conclusion
We have presented a novel algorithm for discrete particle simulations that couples the lattice Boltzmann method with the discrete element method. The main feature of such simulations is that the solid
phase consists of Lagrangian particles that are generally smaller than a computational cell used in the fluid
simulation. The coupling was realized by evaluating interphase momentum transfer forces which include
drag, pressure gradient, lift, and added mass forces. Here, special care and measures had to be taken inside
the algorithm to prevent destabilizing oscillations. The DEM was used to model the collision between the
17

particles, accounting for normal and tangential collision forces. For the description of the fluid flow, a LBM
formulation of the generalized Navier–Stokes equations was chosen to include volume displacement effects
due to the presence of the particulate phase. Furthermore, we used a modified fluid viscosity based on the
local fluid volume fraction and a Smagorinsky–type LES turbulence model. The simulation results of our
approach agreed well with semi–empirical predictions for the settling velocity of spheres over a large range
of solid volume fractions, from a single sphere up to densely packed scenarios. Additionally, we explicitly
computed interparticle lubrication forces to obtain the correct collision behavior for submerged particles.
The major benefit of such a DPS approach in comparison to fully resolved simulations is that significantly
more particles can be simulated in less time. This enables efficient simulations of larger systems like fluidized
beds or riverbeds as they are relevant in industrial and engineering applications. In order to obtain predictive
results with the DPS technique, a careful validation of the method is essential. For our approach, we
demonstrated accurate simulation results for different configurations, demonstrating their suitability for
more complex scenarios. This approach also offers the flexibility to be further extended by exchanging the
applied models and adding further ones. For example, even though we focused on spherical particles in
this work, the presented approach can be extended to other shapes, which then allows for a more realistic
representation of e.g. river sediments, as done in [38]. Also, this approach can be used in combination with
a fully resolved simulation approach such that polydisperse systems with large particle size ratios can be
simulated. Then, the large particles could be geometrically fully resolved and the smaller particle are treated
by DPS.
References
[1] J. Götz, K. Iglberger, M. Stürmer, U. Rüde, Direct numerical simulation of particulate flows on 294912 processor cores,
in: Proceedings of the 2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage
and Analysis, IEEE Computer Society, 2010, pp. 1–11. doi:10.1109/sc.2010.20.
[2] Q. Xiong, B. Li, G. Zhou, X. Fang, J. Xu, J. Wang, X. He, X. Wang, L. Wang, W. Ge, J. Li, Large–scale DNS of gas–solid
flows on Mole–8.5, Chemical Engineering Science 71 (2012) 422–430. doi:10.1016/j.ces.2011.10.059.
[3] A. Esteghamatian, A. Hammouti, M. Lance, A. Wachs, Particle resolved simulations of liquid/solid and gas/solid fluidized
beds, Physics of Fluids 29 (3) (2017) 033302. doi:10.1063/1.4979137.
[4] A. G. Kidanemariam, M. Uhlmann, Interface–resolved direct numerical simulation of the erosion of a sediment bed sheared
by laminar channel flow, International Journal of Multiphase Flow 67 (2014) 174–188. doi:10.1016/j.ijmultiphaseflow.
2014.08.008.
[5] C. Rettinger, C. Godenschwager, S. Eibl, T. Preclik, T. Schruff, R. Frings, U. Rüde, Fully resolved simulations of dune
formation in riverbeds, in: International Conference on High Performance Computing, Springer International Publishing,
2017, pp. 3–21. doi:10.1007/978-3-319-58667-0_1.
[6] A. G. Kidanemariam, M. Uhlmann, Formation of sediment patterns in channel flow: minimal unstable systems and their
temporal evolution, Journal of Fluid Mechanics 818 (2017) 716–743. doi:10.1017/jfm.2017.147.
[7] P. Seil, S. Pirker, T. Lichtenegger, Onset of sediment transport in mono– and bidisperse beds under turbulent shear flow,
Computational Particle Mechanicsdoi:10.1007/s40571-017-0163-6.
[8] M. van der Hoef, M. van Sint Annaland, N. Deen, J. Kuipers, Numerical simulation of dense gas–solid fluidized beds: A
multiscale modeling strategy, Annual Review of Fluid Mechanics 40 (1) (2008) 47–70. doi:10.1146/annurev.fluid.40.
111406.102130.
[9] P. A. Cundall, O. D. Strack, A discrete numerical model for granular assemblies, geotechnique 29 (1) (1979) 47–65.
[10] S. Bogner, S. Mohanty, U. Rüde, Drag correlation for dilute and moderately dense fluid–particle systems using the
lattice Boltzmann method, International Journal of Multiphase Flow 68 (Supplement C) (2015) 71–79. doi:10.1016/j.
ijmultiphaseflow.2014.10.001.
[11] Y. (Yali) Tang, E. A. J. F. (Frank) Peters, J. A. M. (Hans) Kuipers, S. H. L. (Sebastian) Kriebitzsch, M. A. (Martin)
van der Hoef, A new drag correlation from fully resolved simulations of flow past monodisperse static arrays of spheres,
AIChE Journal 61 (2) (2015) 688–698. doi:10.1002/aic.14645.
[12] S. Subramaniam, Lagrangian–Eulerian methods for multiphase flows, Progress in Energy and Combustion Science 39 (2)
(2013) 215–245. doi:10.1016/j.pecs.2012.10.003.
[13] W. Zhong, A. Yu, G. Zhou, J. Xie, H. Zhang, CFD simulation of dense particulate reaction system: Approaches, recent
advances and applications, Chemical Engineering Science 140 (Supplement C) (2016) 16–43. doi:10.1016/j.ces.2015.
09.035.
[14] W. Ge, L. Wang, J. Xu, F. Chen, G. Zhou, L. Lu, Q. Chang, J. Li, Discrete simulation of granular and particlefluid flows: from fundamental study to engineering application, Reviews in Chemical Engineering 0 (0). doi:10.1515/
revce-2015-0079.
[15] S. Chen, G. D. Doolen, Lattice Boltzmann Method for Fluid Flows, Annual Review of Fluid Mechanics 30 (1) (1998)
329–364. doi:10.1146/annurev.fluid.30.1.329.

18

[16] C. Godenschwager, F. Schornbaum, M. Bauer, H. Köstler, U. Rüde, A framework for hybrid parallel flow simulations with
a trillion cells in complex geometries, in: Proceedings of the International Conference on High Performance Computing,
Networking, Storage and Analysis, ACM Press, 2013, pp. 1–12. doi:10.1145/2503210.2503273.
[17] P. Ahlrichs, B. Dünweg, Simulation of a single polymer chain in solution by combining lattice Boltzmann and molecular
dynamics, The Journal of Chemical Physics 111 (17) (1999) 8225–8239. doi:10.1063/1.480156.
[18] H. Zhang, F. X. Trias, A. Oliva, D. Yang, Y. Tan, S. Shu, Y. Sheng, PIBM: Particulate immersed boundary method for
fluid–particle interaction problems, Powder Technology 272 (2015) 1–13. doi:10.1016/j.powtec.2014.11.025.
[19] Q. Xiong, E. Madadi-Kandjani, G. Lorenzini, A LBM–DEM solver for fast discrete particle simulation of particle–fluid
flows, Continuum Mechanics and Thermodynamics 26 (6) (2014) 907–917. doi:10.1007/s00161-014-0351-z.
[20] L. Wang, B. Zhang, X. Wang, W. Ge, J. Li, Lattice Boltzmann based discrete simulation for gas–solid fluidization,
Chemical Engineering Science 101 (2013) 228–239. doi:10.1016/j.ces.2013.06.019.
[21] B. Dünweg, A. J. C. Ladd, Lattice Boltzmann simulations of soft matter systems, in: Advanced Computer Simulation
Approaches for Soft Matter Sciences III, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009, pp. 89–166. doi:10.1007/
978-3-540-87706-6_2.
[22] X. Niu, C. Shu, Y. Chew, Y. Peng, A momentum exchange–based immersed boundary–lattice Boltzmann method for
simulating incompressible viscous flows, Physics Letters A 354 (3) (2006) 173–182. doi:10.1016/j.physleta.2006.01.060.
[23] Y. Hu, H. Yuan, S. Shu, X. Niu, M. Li, An improved momentum exchanged–based immersed boundary–lattice Boltzmann
method by using an iterative technique, Computers & Mathematics with Applications 68 (3) (2014) 140–155. doi:
10.1016/j.camwa.2014.05.013.
[24] M. A. Habte, C. Wu, Particle sedimentation using hybrid lattice Boltzmann–immersed boundary method scheme, Powder
Technology 315 (2017) 486–498. doi:10.1016/j.powtec.2017.04.032.
[25] M.-L. Maier, T. Henn, G. Thaeter, H. Nirschl, M. J. Krause, Multiscale simulation with a two–way coupled lattice
Boltzmann method and discrete element method, Chemical Engineering & Technology 40 (9) (2017) 1591–1598. doi:
10.1002/ceat.201600547.
[26] A. J. Cihonski, J. R. Finn, S. V. Apte, Volume displacement effects during bubble entrainment in a travelling vortex ring,
Journal of Fluid Mechanics 721 (2013) 225–267. doi:10.1017/jfm.2013.32.
[27] D. Joseph, T. Lundgren, R. Jackson, D. Saville, Ensemble averaged and mixture theory equations for incompressible fluid–
particle suspensions, International Journal of Multiphase Flow 16 (1) (1990) 35–42. doi:10.1016/0301-9322(90)90035-H.
[28] J. Kuipers, K. V. Duin, F. V. Beckum, W. V. Swaaij, A numerical model of gas–fluidized beds, Chemical Engineering
Science 47 (8) (1992) 1913–1924. doi:10.1016/0009-2509(92)80309-Z.
[29] T. Wang, J. Wang, Two–fluid model based on the lattice Boltzmann equation, Phys. Rev. E 71 (2005) 045301. doi:
10.1103/PhysRevE.71.045301.
[30] J. Zhang, L. Wang, J. Ouyang, Lattice Boltzmann model for the volume–averaged Navier–Stokes equations, EPL (Europhysics Letters) 107 (2) (2014) 20001. doi:10.1209/0295-5075/107/20001.
[31] B. Blais, J.-M. Tucny, D. Vidal, F. Bertrand, A conservative lattice Boltzmann model for the volume–averaged Navier–
Stokes equations based on a novel collision operator, Journal of Computational Physics 294 (2015) 258–273. doi:10.1016/
j.jcp.2015.03.036.
[32] F. Song, W. Wang, J. Li, A lattice Boltzmann method for particle–fluid two–phase flow, Chemical Engineering Science
102 (2013) 442–450. doi:10.1016/j.ces.2013.08.037.
[33] D. R. Noble, J. R. Torczynski, A Lattice–Boltzmann Method for Partially Saturated Computational Cells, International
Journal of Modern Physics C 09 (08) (1998) 1189–1201. doi:10.1142/S0129183198001084.
[34] C. Rettinger, U. Rüde, A comparative study of fluid–particle coupling methods for fully resolved lattice Boltzmann
simulations, Computers & Fluids 154 (2017) 74–89, iCCFD8. doi:10.1016/j.compfluid.2017.05.033.
[35] T. Preclik, U. Rüde, Ultrascale simulations of non–smooth granular dynamics, Computational Particle Mechanics 2 (2)
(2015) 173–196. doi:10.1007/s40571-015-0047-6.
[36] K. Iglberger, U. Rüde, Massively parallel granular flow simulations with non–spherical particles, Computer Science Research and Development 25 (1) (2010) 105–113. doi:10.1007/s00450-010-0114-4.
[37] C. T. Crowe, J. D. Schwarzkopf, M. Sommerfeld, Y. Tsuji, Multiphase flows with droplets and particles, CRC press, 2011.
[38] J. R. Finn, M. Li, S. V. Apte, Particle based modelling and simulation of natural sand dynamics in the wave bottom
boundary layer, Journal of Fluid Mechanics 796 (2016) 340–385. doi:10.1017/jfm.2016.246.
[39] R. Sun, H. Xiao, SediFoam: A general–purpose, open–source CFD-–DEM solver for particle–laden flow with emphasis on
sediment transport, Computers & Geosciences 89 (2016) 207–219. doi:10.1016/j.cageo.2016.01.011.
[40] A. M. Roma, C. S. Peskin, M. J. Berger, An adaptive version of the immersed boundary method, Journal of Computational
Physics 153 (2) (1999) 509–534. doi:10.1006/jcph.1999.6293.
[41] X. Yang, X. Zhang, Z. Li, G.-W. He, A smoothing technique for discrete delta functions with application to immersed
boundary method in moving boundary simulations, Journal of Computational Physics 228 (20) (2009) 7821–7836. doi:
10.1016/j.jcp.2009.07.023.
[42] S. Tenneti, R. Garg, S. Subramaniam, Drag law for monodisperse gas–solid systems using particle–resolved direct numerical
simulation of flow past fixed assemblies of spheres, International Journal of Multiphase Flow 37 (9) (2011) 1072–1092.
doi:10.1016/j.ijmultiphaseflow.2011.05.010.
[43] T. R. Auton, J. C. R. Hunt, M. Prud’Homme, The force exerted on a body in inviscid unsteady non–uniform rotational
flow, Journal of Fluid Mechanics 197 (1988) 241–257. doi:10.1017/S0022112088003246.
[44] P. G. Saffman, The lift on a small sphere in a slow shear flow, Journal of Fluid Mechanics 22 (2) (1965) 385–400.
doi:10.1017/S0022112065000824.
[45] R. Ball, J. Melrose, A simulation technique for many spheres in quasi–static motion under frame–invariant pair drag

19

[46]
[47]
[48]

[49]
[50]
[51]
[52]
[53]
[54]

[55]
[56]
[57]
[58]
[59]
[60]
[61]

and Brownian forces, Physica A: Statistical Mechanics and its Applications 247 (1) (1997) 444–472. doi:10.1016/
S0378-4371(97)00412-3.
N.-Q. Nguyen, A. J. C. Ladd, Lubrication corrections for lattice–Boltzmann simulations of particle suspensions, Phys.
Rev. E 66 (2002) 046708. doi:10.1103/PhysRevE.66.046708.
D. Bartuschat, U. Rüde, Parallel multiphysics simulations of charged particles in microfluidic flows, Journal of Computational Science 8 (Supplement C) (2015) 1–19. doi:10.1016/j.jocs.2015.02.006.
P. Nithiarasu, K. Seetharamu, T. Sundararajan, Natural convective heat transfer in a fluid saturated variable porosity
medium, International Journal of Heat and Mass Transfer 40 (16) (1997) 3955–3967. doi:10.1016/S0017-9310(97)
00008-2.
Z. Guo, T. S. Zhao, Lattice Boltzmann model for incompressible flows through porous media, Phys. Rev. E 66 (2002)
036304. doi:10.1103/PhysRevE.66.036304.
S. Ergun, Fluid flow through packed columns, Chem Eng Prog 48 (1952) 89.
F. Ferrini, D. Ercolani, B. de Cindio, L. Nicodemo, L. Nicolais, S. Ranaudo, Shear viscosity of settling suspensions,
Rheologica Acta 18 (2) (1979) 289–296. doi:10.1007/BF01542776.
Y. H. Qian, D. D’Humières, P. Lallemand, Lattice BGK Models for Navier–Stokes Equation, EPL (Europhysics Letters)
17 (6) (1992) 479. doi:10.1209/0295-5075/17/6/001.
Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Physical Review
E 65 (4) (2002) 046308. doi:10.1103/PhysRevE.65.046308.
S. Hou, J. Sterling, S. Chen, G. Doolen, A lattice Boltzmann subgrid model for high reynolds number flows, in: A. T.
Lawniczak, R. Kapral (Eds.), Pattern formation and lattice gas automata, Vol. 6 of Fields Institute Communications,
American Mathematical Society, Providence, RI, 1996, p. 151–166.
H. Yu, S. S. Girimaji, L.-S. Luo, DNS and LES of decaying isotropic turbulence with and without frame rotation using
lattice Boltzmann method, Journal of Computational Physics 209 (2) (2005) 599–616. doi:10.1016/j.jcp.2005.03.022.
R. Ramadugu, S. P. Thampi, R. Adhikari, S. Succi, S. Ansumali, Lattice differential operators for computational physics,
EPL (Europhysics Letters) 101 (5) (2013) 50006. doi:10.1209/0295-5075/101/50006.
J. Richardson, W. Zaki, The sedimentation of a suspension of uniform spheres under conditions of viscous flow, Chemical
Engineering Science 3 (2) (1954) 65–73. doi:10.1016/0009-2509(54)85015-9.
T. Baldock, M. Tomkins, P. Nielsen, M. Hughes, Settling velocity of sediments at high concentrations, Coastal Engineering
51 (1) (2004) 91–100. doi:10.1016/j.coastaleng.2003.12.004.
G. G. Joseph, R. Zenit, M. L. Hunt, A. M. Rosenwinkel, Particle–wall collisions in a viscous fluid, Journal of Fluid
Mechanics 433 (2001) 329–346. doi:10.1017/S0022112001003470.
P. Gondret, M. Lance, L. Petit, Bouncing motion of spherical particles in fluids, Physics of Fluids 14 (2) (2002) 643–652.
doi:10.1063/1.1427920.
A. Kidanemariam, M. Uhlmann, Direct numerical simulation of pattern formation in subaqueous sediment, Journal of
Fluid Mechanics 750. doi:10.1017/jfm.2014.284.

20

