SIMULATION OF IRREGULARLY SHAPED PARTICLES USING COUPLING
METHOD OF LATTICE BOLTZMANN AND DISCRETE ELEMENT MODELLING
by
Mohammadhassan Ahmadian
B.S., Imam Hossein University, 2013
M.S., Ferdowsi University of Mashhad, 2018

THESIS SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
IN
NATURAL RESOURCES AND ENVIRONMENTAL STUDIES

UNIVERSITY OF NORTHERN BRITISH COLUMBIA
April 2024
© Mohammadhassan Ahmadian, 2024

Abstract
In this study, the lattice Boltzmann method (LBM) and the discrete element method (DEM)
are coupled to simulate the interaction between the fluid phase and irregularly shaped solid
particles. For this purpose, the geometry of real particle shapes is represented as clumps of
overlapping spheres and then simulated through a multi-sphere model in DEM and coupled with
LBM using two open-access codes, LIGGGHTS and Palabos. The accuracy of the coupling
method with clumps is demonstrated by simulating several benchmark cases and comparing them
with the results from the literature. The coupled LBM-DEM method is then used to simulate the
collapse and transport of submerged granular particles with four different irregular particle shapes
in addition to sphere particle shapes to highlight the influence of grain morphology in the solidfluid interaction. To the authors’ best knowledge, it is the first time that LIGGGHTS and Palabos
were used to simulate the LBM-DEM coupling of irregular particles. This study demonstrates the
potential and accuracy of LBM-DEM method in simulating particles with irregular shapes. The
research also provides insight into the importance of considering real shape particles in simulating
granular flows in geotechnical engineering.

ii

Table of Contents

Abstract ........................................................................................................................................... ii
Table of Contents ........................................................................................................................... iii
List of Tables .................................................................................................................................. v
List of Figures ................................................................................................................................ vi
Glossary ......................................................................................................................................... ix
Acknowledgement .......................................................................................................................... x
Chapter 1: Introduction ................................................................................................................... 1
1.1 Motivation and background .................................................................................................. 1
1.2 Literature review ................................................................................................................... 2
1.2.1 Numerical approaches for solid-fluid interaction .......................................................... 2
1.2.2 Eulerian-Lagrangian methods ........................................................................................ 5
1.2.3 LBM-DEM coupling method ......................................................................................... 8
1.2.4 Simulation of irregularly shaped particles ................................................................... 10
1.3 Research objective and current study ................................................................................. 11
1.4 Thesis organization ............................................................................................................. 12
Chapter 2: Methodology ............................................................................................................... 14
2.1 Introduction ......................................................................................................................... 14
2.2 Lattice Boltzmann method .................................................................................................. 14
iii

2.3 Discrete element method..................................................................................................... 17
2.4 Immersed moving boundary method .................................................................................. 21
2.5 LBM-DEM coupling scheme .............................................................................................. 24
2.6 Irregular shape simulation................................................................................................... 26
Chapter 3: Validation .................................................................................................................... 31
3.1 Introduction ......................................................................................................................... 31
3.2 Mesh sensitivity, relaxation time, and sub-particle study ................................................... 31
3.3 Benchmark Cases ................................................................................................................ 34
3.3.1 Settling of a single clump ............................................................................................ 34
3.3.2 Settling of two clumps ................................................................................................. 37
3.3.3 Settling of a disc-shaped clump ................................................................................... 40
Chapter 4: Immersed granular collapse ........................................................................................ 43
4.1 Introduction ......................................................................................................................... 43
4.2 Irregularly shaped particles ................................................................................................. 43
4.3 Simulation of column collapse............................................................................................ 47
Chapter 5: Conclusion................................................................................................................... 56
5.1 Contribution ........................................................................................................................ 56
5.2 Limitations .......................................................................................................................... 57
5.3 Future work ......................................................................................................................... 57
References ..................................................................................................................................... 59
iv

List of Tables
Table 2.1: Values of 𝐜𝒊 , 𝑤𝑖 for a D3Q19 lattice structure (G. C. Yang, Jing, Kwok, & Sobral,
2019). ............................................................................................................................................ 16
Table 2.2: Parameters used in CLUMP algorithm for generating clumps .................................... 29
Table 3.1: Non-dimensional number values for two cases of disc-shaped particle benchmark. .. 41
Table 3.2: Characteristics of two cases related to disc-shaped particle benchmark. .................... 41
Table 4.1: 3D morphological descriptors (Zheng et al., 2021) ..................................................... 44
Table 4.2: Morphological descriptors of irregularly shaped particles .......................................... 45
Table 4.3: Characteristics of fluid and solid ................................................................................. 48

v

List of Figures
Figure 1.1: Solid-fluid interaction simulation methods (Ibrahim & Meguid, 2020) ...................... 3
Figure 2.1: D3Q19 lattice structure .............................................................................................. 15
Figure 2.2: Interaction of two particles in DEM a) when two particles are in contact, b) the springdashpot model ............................................................................................................................... 18
Figure 2.3: The process of calculating solid fraction of cells covered by a clump of overlapping
spheres. In this figure, the red box is the bounding box of each sub-particles. ............................ 22
Figure 2.4: The left figure shows how an irregularly shaped particle is modelled by clump. The
right figure depicts how a clump is mapped on an LBM grid. The zoomed-in cell is a cell partially
saturated by solids and is divided into five sub-cells in each direction. ....................................... 22
Figure 2.5: Flowchart of LBM-DEM coupling method (G. C. Yang et al., 2019) ....................... 25
Figure 2.6: CLUMP algorithm flowchart ..................................................................................... 28
Figure 2.7: Illustration of original STL file (left column), corresponding clump (middle column)
and the subsequent shape (right column) of particles after transformation into clump ................ 29
Figure 2.8: Monte Carlo algorithm for calculating the center of mass and moment of inertia in
clumps ........................................................................................................................................... 30
Figure 3.1: Domain of sedimentation of a single clump. All boundaries are considered as no-slip
walls. ............................................................................................................................................. 32
Figure 3.2: Study the effect of gird resolution on the simulation of sedimentation of a spherical
particle in a tank of fluid. .............................................................................................................. 33
Figure 3.3: Study the effect of relaxation time on the simulation of sedimentation of a spherical
particle in a tank of fluid. .............................................................................................................. 33

vi

Figure 3.4: Sedimentation velocity of an irregularly shaped particle with different numbers of subparticles in the clump. The actual shape of the particle is shown at the side................................ 34
Figure 3.5: Comparison of current simulation and experimental results (ten Cate et al., 2002) of a
single clump settling in ambient fluid. .......................................................................................... 35
Figure 3.6: Comparing the simulation results with the experiment (ten Cate et al., 2002) for settling
trajectory of the sphere particle at various Reynolds numbers ..................................................... 36
Figure 3.7: Vertical velocity contour of a clump at different times (Re = 11.6) .......................... 37
Figure 3.8: Domain of the benchmark case settling of two clumps. All boundaries are considered
as no-slip walls. ............................................................................................................................. 38
Figure 3.9: Settling velocity of two clumped particles in comparison with the results from (Dash
& Lee, 2015) ................................................................................................................................. 39
Figure 3.10: Sedimentation of two clumps at different times ....................................................... 39
Figure 3.11: Domain of the benchmark case disc shape clump. All boundaries are considered as
no-slip walls. ................................................................................................................................. 40
Figure 3.12: Settling velocity of disc-shaped particle in two different cases. .............................. 42
Figure 3.13: Z-velocity contour of settling disc-shaped clump at different times (Case 2). ...... 42
Figure 4.1: Zingg diagram. (a) Jordan No. 1, (b) Jordan No. 26, (c) Jordan No. 59, (d) Jordan No.
70................................................................................................................................................... 44
Figure 4.2: Shape of Jordan sands used in this study ................................................................... 45
Figure 4.3: Sedimentation velocity of four different Jordan sand particles .................................. 46
Figure 4.4: Sedimentation of four irregularly shaped particles at t=0.12 s. (a) particle No. 1, (b)
particle No. 26, (c) particle No. 59, (d) particle No. 70. ............................................................... 47

vii

Figure 4.5. Domain dimensions for immersed column collapse. All boundaries are no-slip walls.
....................................................................................................................................................... 48
Figure 4.6: Collapse of column of irregularly shaped particle No. 1 at different times a) t = 0.08 s,
b) t = 0.16 s, c) t = 0.39 s, d) t = 0.52 s ......................................................................................... 49
Figure 4.7: Collapse of column of irregularly shaped particle No. 26 at different times a) t = 0.08
s, b) t = 0.16 s, c) t = 0.33 s, d) t = 0.44 s ..................................................................................... 50
Figure 4.8: Collapse of column of irregularly shaped particle No. 59 at different times a) t = 0.08
s, b) t = 0.16 s, c) t = 0.46 s, d) t = 0.62 s ..................................................................................... 51
Figure 4.9: Collapse of column of irregularly shaped particle No. 70 at different times a) t = 0.08
s, b) t = 0.16 s, c) t = 0.33 s, d) t = 0.44 s ..................................................................................... 52
Figure 4.10: Collapse of column of spherical particle at different times a) t = 0.08 s, b) t = 0.16 s,
c) t = 0.25 s, d) t = 0.33 s .............................................................................................................. 53
Figure 4.11: Final deposition of column collapse of spherical particles. 𝐿𝑓 = 3.45053 𝑐𝑚 ....... 54
Figure 4.12: Final deposition of column collapse of a) irregularly shaped particle No. 1 with 𝐿𝑓 =
3.24631 𝑐𝑚, b) irregularly shaped particle No. 26 with 𝐿𝑓 = 3.45091 𝑐𝑚, c) irregularly shaped
particle No. 59 with 𝐿𝑓 = 3.78061 𝑐𝑚 , d) irregularly shaped particle No. 70 with 𝐿𝑓 =
3.31487 𝑐𝑚 .................................................................................................................................. 55

viii

Glossary
BGK

Bhatnagar-Gross-Krook

CFD

Computational Fluid Dynamics

CFD-DEM

Coupling of Computational Fluid Dynamics and Discrete Element Method

CLUMP

A Code Library to generate Universal Multi-sphere Particles

D3Q19

A 3D lattice structure in LBM

DEM

Discrete Element Method

FVM

Finite Volume Method

IMB

Immersed Moving Boundary

LBM

Lattice Boltzmann Method

LBM-DEM

Coupling of Lattice Boltzmann Method and Discrete Element Method

LIGGGHTS

Open-Source Discrete Element Method Particle Simulation Software

MPM

Material Point Method

MPS

Moving Particle Semi-Implicit

Palabos

Parallel Lattice Boltzmann Solver

SPH

Smoothed Particle Hydrodynamics

TFM

Two-Fluid Model

ix

Acknowledgement
I would like to express my deepest gratitude to my supervisor, Dr. Wenbo Zheng, for his
invaluable guidance, patience, and support throughout the duration of my research over the past
two years. His wisdom and mentorship have been pivotal to my development and success in my
studies.
My sincere thanks also go to Dr. Lu Jing, who was instrumental in helping me compile the
LBM-DEM code. His willingness to answer my questions and provide assistance during my initial
phase of engaging with the current code has been greatly appreciated.
A special note of thanks is reserved for Dr. Nikhil Aravindakshan, the UNBC HPC Senior
Lab Instructor. His expertise and assistance in compiling the code on the UNBC workstation
computer were crucial to my research's progress, and for that, I am truly thankful.
I am also grateful to Natural Sciences and Engineering Research Council of Canada (NSERC)
for providing the funding that supported my master's program. This financial support has been
fundamental to my research and academic endeavors.
My heartfelt appreciation extends to my committee members, Dr. Jueyi Sui, Dr. Fan Jiang,
and Dr. Sumi Siddiqua, in advance, for their time, insightful comments, and constructive criticism
on my thesis. Their feedback has been invaluable in enhancing the quality of my work.

x

Chapter 1: Introduction
1.1 Motivation and background
The fluid-solid interaction is a fundamental concept in geotechnical engineering with practical
importance in various issues involving water and soils, such as scouring phenomena around
offshore mono-pile foundations (Nielsen, Probst, Petersen, & Sumer, 2015) and internal erosion
in embankment dams (Foster, Fell, & Spannagle, 2000; Teng, Johansson, & Hellström, 2023). The
interaction between fluids and sediment particles, characterized by a strong coupling effect, plays
a crucial role in not only shaping particle movements but also influencing the fluid flow patterns
(Jansen & Harting, 2011; Stratford, Adhikari, Pagonabarraga, & Desplat, 2005; T. Wang, Zhang,
& Zheng, 2023). The essence of fluid-solid system behaviour lies at the microscopic scale, where
the interaction between solid particles and fluid depends on parameters like fluid velocity, solid
volume fraction, and particle shape. Such parameters eventually contribute to the interaction in
terms of forces and momentum exchange. As conventional experimental measurements are unable
to capture the intricate interplay of fluid and particle interactions, employing numerical modelling
is more promising in studying the underlying mechanisms involved in a fluid-solid system (Lai,
Zhao, Zhao, & Huang, 2023).
Most research simulating the interaction of solid particles and fluid simplifies the shape of
particles as spheres (Lai et al., 2023). However, the influence of particle irregularity on the overall
characteristics of a granular flow, including its packing density, stiffness, and strength, is widely
acknowledged as significant (Cho, Dodds, & Santamarina, 2007). In addition, the effect of
particle’s shape might get even more apparent as the goal is to eventually simulate real-world
problems such as geohazards where the systems are considered at macroscopic scale. Therefore,
it is crucial to study more on the feasibility of current computational methods in simulating system
1

of real shape particles submerged in fluid and consequently improving the level of accuracy in
granular flows.
1.2 Literature review
Finding solutions for the governing equations of granular flows through analytical methods
becomes inherently complicated because of nonlinear and complex boundary conditions nature of
such equations in real-world situations (Ibrahim & Meguid, 2020). As a result, the conventional
recourse lies in employing numerical techniques to address these equations. Numerical analysis
presents researchers with two fundamental approaches: the Eulerian method and the Lagrangian
method (Ibrahim & Meguid, 2020). Within the Eulerian framework, flow variables are
conceptualized as a continuous entity within spatially fixed or moving grids. The evolution of these
variables over time is tracked locally within each computational cell. Conversely, the Lagrangian
approach involves the tracing of trajectories and various flow parameters concerning fluid or solid
particles over time, considering the behaviour of each individual particle.
1.2.1 Numerical approaches for solid-fluid interaction
The complexity of granular flows necessitates a nuanced classification, which predominantly
depends on the numerical treatment of both the fluid and solid phases. This classification as shown
in Figure 1.1 gives rise to three main categories in granular flow modelling: (1) Eulerian-Eulerian,
where both phases are treated using the Eulerian framework; (2) Eulerian-Lagrangian, where the
fluid phase is handled using the Eulerian method while the solid phase follows the Lagrangian
approach; and (3) Lagrangian-Lagrangian, wherein both phases are approached from a Lagrangian
perspective (Ibrahim & Meguid, 2020). Each of these categories offers distinct advantages and
limitations, making them suitable for specific scenarios and research objectives. Researchers often
choose these methods based on the intricacies of the problem at hand, aiming to strike a balance
2

between computational efficiency and accuracy in capturing the complexities of solid-fluid
interactions.

Figure 1.1: Solid-fluid interaction simulation methods (Ibrahim & Meguid, 2020)
The Eulerian-Eulerian approach, also known as purely Eulerian approach relies on averaging
of flow variables pertaining to both solid particles and fluids. Among the array of averaging
techniques available, volume averaging emerges as a prominent method, where computations
occurs within finite volumes, delineated as computational cells (Crowe, Schwarzkopf,
Sommerfeld, & Tsuji, 2011). A known model in using Eulerian-Eulerian approach is Two-Fluid
Model (TFM) (L. Yang, Padding, & Kuipers, 2016). This model is a direct product of volume
averaging methodologies, offering a unique perspective on granular flows. In the TFM, the solid
phase within the flowing medium is conceptualized as a fluid-like continuum.
A defining characteristic of the TFM lies in its treatment of fluid and solid phases as one
continuum medium (L. Yang et al., 2016). This intricate concept postulates that both phases exist
as continuous entities within the same averaging volume, albeit partially occupying it. This notion

3

challenges visualization due to the coexistence of the phases within the same space yet occupying
only specific portions of that space.
However, despite its applicability to specific simulation cases, the TFM grapples with
limitations, particularly in its fundamental depiction of particulate dynamics at the microscale. A
critical shortcoming arises from the indirect representation of particle-particle and particle-fluid
interplay using averaged values derived from closures (Crowe et al., 2011). While the averaging
technique could be helpful in gaseous flows involving solid particles of similar size to the
averaging of volume, its application cannot be applied to cases such as dispersed solids within a
fluid. The reason is that TFM lacks the granularity necessary to capture the intricacies of particulate
interactions at a microscale level. In cases where precise microscale descriptions of granular
interactions are crucial, such as in studies of sedimentation, particle aggregation, or multiphase
flows with significant size disparities between particles, more sophisticated and nuanced
modelling approaches are imperative (Ibrahim & Meguid, 2020).
In contrast, fully Lagrangian methods take a markedly different route by adopting tracking of
fluid and solid phase trajectories in granular flows. This Lagrangian perspective offers a
comprehensive understanding of the movement and interaction of individual particles within the
flow. One approach to implementing such tracking involves the utilization of a dynamically
adaptive grid, where the grid adjusts its structure to follow the movement of fluid particles
throughout the simulation process (Ibrahim & Meguid, 2020).
Alternatively, the fluid can be represented by a collection of discrete particles devoid of any
mesh structure. Methods such as Smoothed Particle Hydrodynamics (SPH) and Moving Particle
Semi-Implicit (MPS) methods embrace this meshless or mesh-free paradigm, relying on the
interactions between individual particles to model fluid behaviour accurately (Gingold &
4

Monaghan, 1977). In these approaches, the absence of a predefined grid structure allows for a
more flexible and adaptable representation of fluid dynamics, especially in cases where traditional
grid-based methods may face challenges.
Additionally, there are hybrid methods like the Material Point Method (MPM), which
amalgamate the utilization of a computational mesh for computational calculations (Koshizuka &
Oka, 1996). Despite this grid incorporation, MPM is fundamentally categorized as a particle
method due to its emphasis on tracking particle movements and interactions. This hybrid nature
highlights the versatility within the realm of particulate flow modelling, where methodologies
often blur traditional boundaries, emphasizing the importance of tailoring the approach to the
specific intricacies of the problem at hand (Koshizuka & Oka, 1996). Thus, the landscape of
granular flow simulation continues to evolve, driven by innovative techniques that blend
computational precision with adaptability to capture the rich dynamics of solid-fluid interactions.
1.2.2 Eulerian-Lagrangian methods
Employing Lagrangian tracking in granular flow modelling offers distinct advantages, particularly
in accurately simulating flows of high convection and flows involving free surfaces. Unlike
Eulerian methods, including finite volume and finite difference methods, which are prone to errors
related to discretizing the convective terms, Lagrangian methods circumvent these issues (Ibrahim
& Meguid, 2020). This freedom from discretization errors enhances the precision of modelling,
especially in cases involving rapid fluid movements or the intricate dynamics of free surfaces.
However, it is essential to acknowledge that Lagrangian methods, despite their accuracy
benefits, come with their own set of challenges. Tracking particles on a large scale and employing
semi-implicit solutions, which often necessitate solving Poisson's equation to get the pressure field,
can significantly escalate computational costs (Pozorski & Olejnik, 2023). Additionally, the
5

dynamic adaptation of meshes in Lagrangian methods adds another layer of complexity, further
contributing to higher computational demands. In contrast, the Finite Volume Method (FVM), as
a conventional grid-based method, provides a relatively more streamlined computational approach,
making it advantageous in scenarios where computational efficiency is a primary concern (Ibrahim
& Meguid, 2020).
While Lagrangian methods excel in capturing intricate flow phenomena, their computational
demands prompt researchers to weigh the benefits against the costs, considering factors such as
the scale of the simulation, available computational resources, and the specific intricacies of the
problem at hand.
In the Eulerian-Lagrangian approach, solid particle motion is followed by employing
Newtonian laws of motion, allowing for discrete particle representation without relying on a
computational grid. This method brings the advantage of accurately capturing the discrete
characteristics of solid particles (Ibrahim & Meguid, 2020). Unlike the averaging approach of the
Two-Fluid Model (TFM) that is indirect, Eulerian-Lagrangian approach enables direct
implementation of solid-fluid interaction forces, including drag force and pressure difference
forces. Although the estimation of these forces remains complex, surpassing the challenges of
volume averaging, the Lagrangian approach provides unparalleled insights into microscale
interactions inaccessible through Eulerian-Eulerian methods (Ibrahim & Meguid, 2020).
Various models have been developed for both Lagrangian and Eulerian treatments in granular
flows. For the fluid phase (Eulerian), the Navier-Stokes equations that are locally averaged are
employed commonly (Mao, Wang, Yang, Huang, & Zheng, 2023). Alternatively, the Lattice
Boltzmann method can be employed, modelling fluid flow through streaming and collision
processes on a structured lattice (Sayyari, Ahmadian, & Kim, 2022). For the Lagrangian treatment
6

of solid particles, the Discrete Element Model (DEM) is the model that has been used widely
(Tavarez & Plesha, 2007). These sophisticated approaches collectively refine the state of the
simulation, enabling researchers to explore microscale interactions that were previously beyond
reach using Eulerian-Eulerian methods.
Traditional computational fluid dynamics (CFD) methods, relying on finite element or finite
volume approaches, have been created for multiphase scenarios to replicate intricate fluid-solid
interactions (Chu, Wang, Yu, & Vince, 2009; Peng, Doroodchi, Luo, & Moghtaderi, 2014).
Alongside these established techniques, the lattice Boltzmann method (LBM) has emerged as an
alternative for addressing multiphase fluid dynamics. In contrast to conventional CFD methods
that solve macroscopic conservation equations (mass, momentum, energy) using numerical means,
LBM takes a unique approach by simulating a system of fictitious particles that adhere to the same
conservation laws. LBM offers multiple advantages in terms of intricate boundary scenarios,
handling microscopic interactions, and easier parallelization (Jiang, Matsumura, Ohgi, & Chen,
2021).
The Discrete Element Method (DEM) stands as a widely utilized numerical method for
modelling the flow of discrete materials based on Newton’s principles of motion. The prevalence
of granular flow implementations within industrial processes has spurred the creation of numerous
DEM studies (Cortés & Gil, 2007; Tsuji, 2007; Z. Y. Zhou, Kuang, Chu, & Yu, 2010). Through a
direct approach, DEM has the capability to dynamically simulate the trajectory of each particle in
dense phase flows, accurately tracking both its translational and rotational motions within the
defined domain. This comprehensive simulation accounts for a particle’s interactions,
encompassing contact forces and couples, both among particles and between particles and
boundaries (Almeida, Spogis, & Silva, 2016; El-Emam, Shi, & Zhou, 2019).
7

1.2.3 LBM-DEM coupling method
The development of coupling methods in fluid-particle interactions has been marked by several
studies, starting through 2D simulations of seepage flow (F. Chen, Drumm, & Guiochon, 2011),
piping erosion (Lominé, Scholtès, Sibille, & Poullain, 2013), sedimentation of circular particles
(H. Zhang et al., 2014), and granular column collapse (Soundararajan, 2015). The LBM-DEM
coupling, in particular, stands out for its potential in simulating fluid-particle dynamics. This
preference is attributed not only to LBM's efficiency in parallel computing but also to its versatility
in adopting various coupling techniques, such as immersed boundary (Feng & Michaelides, 2004),
momentum exchange (Ladd, 1994), and immersed moving boundary methods (Noble &
Torczynski, 1998). Moreover, LBM-DEM has seen advancements in accuracy, especially with the
adoption of a two-relaxation-time collision operator, enhanced solid fraction calculation (D. Wang,
Leonardi, & Aminossadati, 2018), and more realistic particle shapes. In exploring the effects of
particle shape, Gardner and Sitar (Gardner & Sitar, 2019) studied the interactions of arbitrarily
shaped polyhedral particles using LBM-DEM. More recently, Rettinger and Rude (Rettinger &
Rüde, 2022) introduced an innovative four-way coupling scheme that significantly enhances
granular flow simulation. This method integrates a combination of the multiple-relaxation-time
LBM with an adaptive strategy to control bulk viscosity near particles, thereby augmenting
simulation stability. Furthermore, the incorporation of polygonal discrete elements within LBM
has enabled the simulation of interactions between convex and concave polygonal particles,
broadening the scope of particle dynamics analysis (M. Wang, Feng, Qu, & Zhao, 2021).
Coupling LBM and DEM has been enabled by accounting for the interaction along the solid
and fluid interface. Noble and Torczynski (Noble & Torczynski, 1998) introduced the concept of
Immersed Moving Boundary (IMB), later combined with LBM and DEM by Cook et al. (Cook,
8

Noble, & Williams, 2004) to simulate particle flow. In this approach, the collision of interface
lattice distribution functions corresponds to the solid fraction, while momentum exchange at the
particle boundary calculates the hydrodynamic force reaction (R. Zhang, Kim, & Dinoy, 2021).
IMB’s compatibility with LBM’s parallel nature and its ease of coding make it a practical choice.
The accurate hydrodynamic force using IMB is easily attainable due to its physical basis. Solid
fraction, crucial for hydrodynamic force calculations, has been the focus of various studies. Owen
et al. (Owen, Leonardi, & Feng, 2011) divided lattice cells into subcells, using their positions
relative to particle radii to determine solid fraction. Galindo-Torres (Galindo-Torres, 2013)
introduced an edge-intersection averaging method to estimate the solid fraction based on
interaction lengths. Another method by Jones and Williams (Jones & Williams, 2017) calculates
the solid fraction directly from volume integrals, applicable to specific interaction cases. A recent
study on nodal solid ratio improvement by Wang et al. (M. Wang, Feng, Owen, & Qu, 2019),
reduced the computing cost by half using Gaussian quadrature for nodal solid ratio calculation. In
a different paper by Wang et al. (M. Wang, Feng, Qu, Tao, & Zhao, 2020), the situation where a
single nodal cell is covered by more than one particle is considered. This situation is particularly
important when an irregularly shaped particle is simulated by a clump of overlapping spheres. The
study proposed a method to resolve the numerical instability that may arise in calculating the
weighting function in such situation. To date, LBM-DEM simulations have been used to simulate
a wide range of engineering scenarios, such as sand production (Y. Han & Cundall, 2013),
sediment settlement (Jiang, Liu, Chen, & Tsuji, 2022), and internal erosion (W. Zhou, Ma, Ma,
Cao, & Cheng, 2020).

9

1.2.4 Simulation of irregularly shaped particles
Most of the abovementioned studies simulated the solid phase as spherical particles. In the field
of granular materials, especially those in geotechnical engineering, the irregular shapes of granular
particles have been recognized as a crucial factor affecting the deformation, strength, and
permeability of granular materials (Zheng, Hu, Tannant, & Zhou, 2021). These particles, diverging
from the traditional spherical form, encompass a wide array of shapes (Zheng, Hu, Tannant, Zhang,
& Xu, 2019), which can be broadly categorized into regular and irregular configurations. Regular
non-spherical particles include ellipsoids, super-ellipsoids, and regular polyhedral, each
possessing distinct geometric characteristics. An important distinction between non-spherical and
spherical particles lies in their particle orientation. Notably, studies by Oschmann (Oschmann,
Hold, & Kruggel-Emden, 2014) investigated elongated particle shapes, such as cylinders, plates,
and cuboids, within fluidized beds, employing numerical simulations to explore the averaged
particle orientation. These investigations show the strong influence of accurate shape
representation on simulation results.
Early investigations into irregular particle shapes employing coupling methods primarily
concentrated on 2D simulations. Han et al. (K. Han, Feng, & Owen, 2007) studied irregular
particles with polygonal shapes using the LBM-DEM coupling method within turbulent flows.
Their findings showed the robustness of the coupled approach up to a certain Reynolds number
threshold. Wachs et al. (Wachs, 2009) conducted relevant research by integrating a DEM solver
with Finite Element methods, comparing the behaviour of disks and triangles in a 2D
sedimentation scenario involving 300 particles within a closed box. This comparative analysis
showed that owing to drag forces, disks settled faster than triangles.

10

As research progressed, a shift towards more realistic simulations of irregular particle shapes
became evident. Shen et al. (Shen, Wang, Huang, & Jin, 2022) conducted a study toward this trend
by employing a resolved CFD-DEM coupling approach to simulate the impact of a dambreak wave
on a rock pile. To represent irregular particle shapes accurately, the study introduced a multisphere clump in the DEM. Moreover, the two-phase interaction involving air and fluid was
resolved through the adoption of the volume of fluid method. In a more recent study, LBM-DEM
was applied to simulate hydraulic plucking using three-dimensional shape of rock particles and
comparing the results with the experiment (Gardner, 2023). The study showed that LBM-DEM
was capable of successfully capturing dominant kinetic modes and transient displacement of the
rock particle from a fractured rock mass. Moreover, Lai et al. (Lai et al., 2023) presentend a signed
distance filed (SDF)-based CFD-DEM for simulating arbitrary-shaped particles in multiphase
flow. The authors verified their proposed model through various benchmark cases. This evolution
in research reflects that to date, there is a research gap specifically focusing on the LBM-DEM
simulation of systems composed of irregularly shaped particles, indicating a novel area in the field.
1.3 Research objective and current study
The main objective of this thesis is to develop a three-dimentional LBM-DEM coupling framework
of a system of irregularly shaped particles using CLUMP to reconstruct the shape of particles, and
then evaluate the accuracy of the method in simulating such particles’ shape. Consequently, the
research tasks for this objective are as follows:
•

Regenerate the geometry of irregularly shaped particles through a clump of overlapping
spheres using CLUMP

•

Investigate the potential and accuracy of LBM-DEM method in simulating irregularly
shaped particles by simulating benchmark cases from the literature
11

•

Study the effect of real shape particles compare to regular shape particles by simulating the
well-known case of granular column collapse.
Therefore, in this study, based on the research objective, an LBM-DEM coupling framework

was developed to simulate the interaction between fluid and solid granular particles with irregular
shapes. Then, according to CLUMP method a MATLAB code was employed to represent
irregularly shaped particles with sub-particles bonded together to form clumps. Two open-source
codes named Palabos and LIGGGHTS were configured for the LBM-DEM coupling of solid
clumps with fluid. Three benchmark cases were considered to validate our model using clumps.
Finally, to study the impact of realistic particle morphology in granular flow simulations, four real
shape Jordan sand particles were chosen to simulate the collapse and transport of a granular column
of irregularly shaped particles immersed in a fluid. Comparison of the results with that from regular
shape particles (spheres) demonstrates the importance of particle shape in granular flow
simulations. This study represents a step towards more realistic modelling of granular flows in
geotechnical engineering.
1.4 Thesis organization
Based on the objective outlined earlier, the organization of the current thesis is as follows. Chapter
2 explains the methodology implemented in this study, including explaining lattice Boltzmann
method, discrete element modelling, and immersed moving boundary. Chapter 3 provides the
benchmark cases used to validate our adopted computational method. Chapter 4 focuses on one
application of current research in simulating a well-known problem in geotechnical engineering,
i.e., the collapse of an immersed granular column. The chapter also investigate the mesh
sensitivity, relaxation time study, and the minimum number of sub-particles that has to be used to

12

reconstruct an irregularly shaped particle. Finally, Chapter 5 provides the conclusion and
limitations of the computational method and some suggestions for future studies.

13

Chapter 2: Methodology
2.1 Introduction
This chapter introduces the computational method adopted throughout this study. The chapter
starts with explaining the Lattice Boltzmann and Discrete Element methods. Then Immersed
Moving Boundary is described as a method for handling solid boundaries. Afterward, the details
of LBM-DEM implementation are provided. Eventually, the CLUMP method is explained for
generating irregularly shaped particles by a cluster of overlapping spheres.
2.2 Lattice Boltzmann method
The lattice Boltzmann method (LBM) is a mesoscopic scale fluid flow simulation method based
on the kinetic theory (S. Chen & Doolen, 1998). In this method, the computational domain is
divided into structured Cartesian lattice nodes, and then a probability distribution function (PDF)
is defined on each lattice node. In other words, 𝑓𝑖 (𝐱, 𝑡) is the PDF of the point 𝐱 at time 𝑡 at
direction 𝑖.
In contrast to conventional CFD methods that solve nonlinear Navier-Stokes partial
differential equations, the LBM takes a different approach. In LBM, the Boltzmann equation
describes how PDFs evolve throughout time. Based on the BGK approximation, named after its
contributors Bhatnagar-Gross-Krook (Bhatnagar, Gross, & Krook, 1954), the well-known
Boltzmann equation turns into Equation (2.1). As a result of this approximation, the nonlinear
Boltzmann equation becomes an explicit linear equation efficient for parallel computing.
1
𝑓𝑖 (𝐱 + 𝒄𝑖 𝛿𝑡 , 𝑡 + 𝛿𝑡 ) − 𝑓𝑖 (𝐱, 𝑡) = − [𝑓𝑖 (𝐱, 𝑡) − 𝑓𝑒𝑞
(𝐱, 𝑡)]
𝑖
𝜏

(2.1)

The governing Equation (2.1) consists of two main processes: streaming and collision. On the
left-hand side (LHS) of the equation, the streaming process occurs. It involves the propagation of
14

the PDFs from one lattice node (𝐱) (Figure 2.1) to its neighboring node (𝐱 + 𝐜𝑖 𝛿𝑡 ) in the direction
𝑖 at the time 𝑡 + 𝛿𝑡 , where 𝐜𝑖 represents the lattice velocity and 𝛿𝑡 is the increment in time t.

Figure 2.1: D3Q19 lattice structure
On the other hand, in Equation (2.1), the right-hand side describes the occurrence of the
collision process. In this process, the PDFs undergo a linear relaxation to the equilibrium
distribution functions (EDFs), 𝑓𝑖𝑒𝑞 (𝑥, 𝑡), using a single dimensionless relaxation time (𝜏). The
EDF employed in this context is the Maxwellian distribution, which uses the Taylor series and can
be expanded regarding the macroscopic fluid velocity (𝐮𝑓 ) (Aidun & Clausen, 2010).
2
𝐜𝑖 . 𝐮𝑓 (𝐜𝑖 . 𝐮𝑓 )
𝑢𝑓 2
𝑒𝑞
𝑓𝑖 = 𝑤𝑖 𝜌𝑓 [1 + 2 +
−
] (𝑖 = 0, … ,18)
𝑐𝑠
2𝑐4𝑠
2𝑐2𝑠

(2.2)

where 𝑤𝑖 is the weight factor related to the corresponding 𝐜𝑖 direction. The speed of sound 𝑐𝑠 in
lattice is 1⁄√3. Note that a D3Q19 lattice structure, as shown in Figure 2.1, is adopted in this
study. The associated lattice velocity and weight factors are shown in Table 2.1.

15

Table 2.1: Values of 𝐜𝑖 , 𝑤𝑖 for a D3Q19 lattice structure (G. C. Yang, Jing, Kwok, & Sobral,
2019).
PDFs, 𝑓𝑖

Lattice velocity, 𝒄𝑖

Weights, 𝑤𝑖

𝑓0

(0,0,0)

1/3

𝑓1 − 𝑓6

(±1,0,0), (0, ±1,0), (0,0, ±1)

1/18

𝑓7 − 𝑓18

(±1, ±1,0), (±1,0, ±1), (0, ±1, ±1)

1/36

After calculating 𝑓𝑖 at each node, the fluid density 𝜌𝑓 and velocity 𝐮𝑓 can be reconstructed
using the following equations according to the conservation laws of mass and momentum (He &
Luo, 1997):
18

(2.3)

𝜌𝑓 = ∑ 𝑓𝑖
𝑖=0
18

(2.4)

𝜌𝑓 𝐮𝑓 = ∑ 𝐜𝑖 𝑓𝑖
𝑖=0

Moreover, the kinematic fluid viscosity, 𝑣𝑓 , is a material property constant and has the
following relationship with the relaxation time, LBM timestep, and lattice spacing (He & Luo,
1997):
1 𝛿𝑥2
𝑣𝑓 = 𝑐𝑠2 (𝜏 − )
2 𝛿𝑡

(2.5)

Pressure as other macroscopic variable is related to fluid density based on the equation of
state (He & Luo, 1997):
𝑝 = 𝑐𝑠2 𝜌𝑓

(2.6)

One of the primary sources of compressibility error in LBM arises from the truncated terms
of Taylor expansion of the equilibrium distribution functions. In order to approximate an
incompressible flow, the Mach number (Ma) must remain below unity (𝑀𝑎 ≪ 1). The necessity
16

for incompressibility in LBM simulations imposes constraints on selecting the LBM relaxation
time. In other words, while based on Equation (2.5), the relaxation time has a bottom threshold of
0.5, the higher relaxation time should be set in a way to ensure that the LBM accurately captures
the behaviour of incompressible flows while minimizing compressibility-related errors (G. C.
Yang et al., 2019).
2.3 Discrete element method
In the discrete element method (DEM), granular materials are represented as a collection of
individual particles, where interactions between neighbouring particles are explicitly computed
using simplified contact laws (Sun, Sakai, & Yamada, 2013). In this study, an open-source DEM
software, LIGGGHTS, was employed to simulate the interactions among solid particles. In each
iteration of the DEM process, Newton’s second law determines the amount of particle translation
and rotation. Each iteration starts with particle-particle contacts assessment, and accordingly, the
corresponding forces are computed. Subsequently, particles’ linear and rotational velocities are
updated. Afterward, new particle positions will be determined based on Newton’s law of motion
(Ding & Xu, 2018).
Figure 2.2 illustrates particle “a” (shown as green) and particle “b” (shown as red) are in
contact with each other, which 𝜹𝑛 denotes the overlap. The overlap can be calculated using the
following equation (G. C. Yang et al., 2019):
(2.7)

𝜹𝑛 = (𝑟𝑎 + 𝑟𝑏 − 𝑟𝑎𝑏 )𝐧

where 𝑟𝑎 and 𝑟𝑏 are radii of particles 𝑎 and 𝑏, respectively. The distance between two particles is
represented by 𝑟𝑎𝑏 , and 𝐧 is the normal unit vector directed toward particle centers.

17

(a)

(b)

Figure 2.2: Interaction of two particles in DEM a) when two particles are in contact, b) the
spring-dashpot model
To calculate the contact force between particles, constitutive contact models can be used, and
this study uses the spring-dashpot model. In this case, the following equation gives the normal
force (Cundall & Strack, 1979):
(2.8)

𝐅𝑛 = 𝑘𝑛 𝛿𝑛 + 𝑐𝑛 ∆𝐮𝑛

where 𝑘𝑛 is stiffness and 𝑐𝑛 is damping coefficient, both in the normal direction. ∆𝐮𝑛 is the
relative normal velocity of two interacting particles. In this study, we chose Hertz contact model
to calculate the stiffness and damping coefficient as follows (Di Renzo & Di Maio, 2004):
4 ∗ ∗
𝐸 √𝑅 𝛿𝑛
3

(2.9)

5
𝑐𝑛 = −2√ 𝛽√𝑚∗ 𝑆𝑛
6

(2.10)

𝑘𝑛 =

where 𝐸 ∗ , 𝑅 ∗ , 𝑚∗ are Young modulus, radius and mass equivalent, 𝛽 is damping ratio, and 𝑆𝑛 can
be calculated using the following equations (Di Renzo & Di Maio, 2004).
1
(1 − 𝜗12 ) (1 − 𝜗22 )
=
+
𝐸∗
𝐸1
𝐸2

(2.11)
18

1
1
1
=
+
∗
𝑅
𝑅1 𝑅2

(2.12)

1
1
1
=
+
∗
𝑚
𝑚1 𝑚2

(2.13)

𝛽=

ln(𝑒)
(2.14)

√ln2 (𝑒) + 𝜋 2

(2.15)

𝑆𝑛 = 2𝐸 ∗ √𝑅 ∗ 𝛿𝑛

In the above equations, 𝐸 is the Young modulus, 𝜗 is the Poisson ration, 𝑒 is the coefficient
of restitution, 𝑚 is the mass, and 𝑅 is the radius of each two interacting particles.
In the same manner, the tangential force is given as follows (Cundall & Strack, 1979):
𝑡𝑐

(2.16)

𝐅𝑡 = −𝑘𝑡 ∫ ∆𝐮𝑡 𝑑𝑡 − 𝑐𝑡 ∆𝐮𝑡
𝑡𝑐,0

where 𝑘𝑡 and 𝑐𝑡 are corresponding stiffness and damping coefficient in the tangential direction,
and ∆𝐮𝑡 denotes the relative tangential velocity of two particles. The integral term in the equation
represents the increment of spring, which accumulates energy resulting from the relative motion
in the tangential direction between contacting particles. This integral accounts for the deformation
of the particle surface (elastic) throughout the time the particles are in contact from 𝑡𝑐,0 to 𝑡𝑐 . The
tangential force exerted between the particles acts in the opposite direction of the displacement in
the tangential direction. Additionally, the tangential force magnitude is restricted by the Coulomb
friction, denoted as 𝜇𝐹𝑛 . The value of 𝜇 corresponds to the smaller of the two friction coefficients
of the contacting particles, and it represents the threshold at which the particles initiate sliding
motion against each other.
Also based on Hertz model adopted in this study, 𝑘𝑡 and 𝑐𝑡 can be calculated as follows (Di
Renzo & Di Maio, 2004):
19

𝑘𝑡 = 8𝐺 ∗ √𝑅 ∗ 𝛿𝑛

(2.17)

5
𝑐𝑡 = −2√ 𝛽√𝑚∗ 𝑆𝑡
6

(2.18)

where the following equations give 𝑆𝑡 and 𝐺 ∗ .
1
2(2 − 𝜗1 )(1 + 𝜗1 ) 2(2 − 𝜗2 )(1 + 𝜗2 )
=
+
∗
𝐺
𝐸1
𝐸2
∗
∗
𝑆𝑡 = 2𝐺 √𝑅 𝛿𝑛

(2.19)
(2.20)

As mentioned, there are different proposed contact models that calculate contact force (𝐅𝑐 =
𝐅𝑛 + 𝐅𝑡 ) as a function of relative velocity and distance. In our study, we adopted the Hertz contact
model (Jing, Yang, Kwok, & Sobral, 2019). To update the velocities of a particle (both linear and
angular), Newton’s second law of motion takes into account the forces (contact force 𝐅𝑐 , gravity
𝐆, fluid drag force 𝐅𝑓 ) and torques (contact torque 𝐓𝑐 , fluid drag torque 𝐓𝑓 ) exerting on the particle.
The linear velocity is updated by considering the net force acting on the particle and dividing it by
its mass, while the angular velocity is updated by considering the net torque acting on the particle
and dividing it by its moment of inertia. By applying Newton’s second law, the particle’s linear
and angular velocities can be appropriately adjusted to account for the combined effects of these
forces and torques.
𝑚𝐚 = 𝐅𝑐 + 𝐆 + 𝐅𝑓

(2.21)

𝐼𝛚̇ = 𝐓𝑐 + 𝐓𝑓

(2.22)

where 𝑚 is mass and 𝐼 is the moment of inertia of particles. 𝐚 and 𝛚̇ are translational and
rotational acceleration of particles, respectively. Eventually, by using Equations (2.21) and (2.22)
via the Verlet method (Verlet, 1967), the particles’ position and orientation can be updated at each
timestep.

20

2.4 Immersed moving boundary method
The Immersed Moving Boundary (IMB) method operates based on the principle of introducing a
new collision operator Ω that depends on the solid ratio 𝜀. While the ideal value of this ratio can
be obtained through a geometric analysis, it often demands significant computational resources
(G. C. Yang et al., 2019). To address this challenge, a cell decomposition method is employed
(Owen et al., 2011), as depicted in Figure 2.3. According to this method, fluid cells that are
partially saturated by solids are subdivided into equal-sized sub-cells, 𝑛𝑠𝑢𝑏 . An algorithm is
applied to the centers of these sub-cells to estimate the solid ratio, 𝜀, by dividing the number of
3
sub-cells located inside the solid boundary over the total sub-cell number 𝑛𝑠𝑢𝑏
with one cell (G.

C. Yang et al., 2019).
For the implementation of IMB, first, each sphere is encased within its own bounding box
(Figure 2.3). The core of this method involves calculating the solid fraction for each cell of the
bounding box in relation to the sphere's center. Specifically, cells within a distance less than the
radius from the center are classified as solid (solid fraction = 1). Conversely, cells beyond the
radius have a solid fraction of 0. In the case, the cell is covered by more that two particles (Figure
2.3), the solid fraction of the cell will be calculated based on the particle that has a higher solid
fraction in the cell. For the partially saturated cells, the method divides them into smaller, equalsized sub-cells, as shown in Figure 2.3 (denoted as 𝑛𝑠𝑢𝑏 ). Then the cell decomposition method is
employed (Owen et al., 2011), as depicted in Figure 2.3. An algorithm is then applied to the centers
of these sub-cells to estimate the solid ratio, 𝜀, by dividing the number of sub-cells located inside
3
the solid boundary over the total sub-cell number 𝑛𝑠𝑢𝑏
with one cell (G. C. Yang et al., 2019). The

solid fraction (𝜀) in these cells is then determined by comparing the number of sub-cells within the
solid boundary to the total number of sub-cells.
21

An important aspect of this process is the assignment of particle IDs. Cells that are fully or
partially solid receive the ID of the occupying particle, while purely fluid cells are assigned an ID
of zero. This approach offers two significant benefits. Firstly, for cells in the overlapping region
of two spheres, the cell is assigned only one particle ID, preventing the duplication of force and
momentum calculations. Secondly, as each cell carries a particle ID, calculating the total
hydrodynamic force on each particle becomes more efficient – the code simply sums the forces
across cells sharing the same particle ID.

Figure 2.3: The process of calculating solid fraction of cells covered by a clump of overlapping
spheres. In this figure, the red box is the bounding box of each sub-particles.

Figure 2.4: The left figure shows how an irregularly shaped particle is modelled by clump. The
right figure depicts how a clump is mapped on an LBM grid. The zoomed-in cell is a cell
partially saturated by solids and is divided into five sub-cells in each direction.
In fluid cells where 𝜀 is equal to 0, the aforementioned fluid hydrodynamic collision occurs,
which is represented by the BGK collision operator Ω 𝑓 , as provided in the Equation (2.1). In solid
cells where the solid ratio is equal to 1, a collision operator given by Noble and Torczynski (Noble
22

& Torczynski, 1998), according to the concept of non-equilibrium bounce-back (Zou & He, 1997),
is employed. This collision operator is denoted as Ω𝑠 , is given by the following expression:
𝑒𝑞
Ω𝑖𝑠 = 𝑓−𝑖 (𝐱, 𝑡) − 𝑓−𝑖
(𝜌𝑓 , 𝐮𝑓 ) + 𝑓𝑖𝑒𝑞 (𝜌𝑓 , 𝐮𝑠 ) − 𝑓𝑖 (𝐱, 𝑡)

(2.23)

where 𝐮𝑠 is the macroscopic velocity of the solid at the location of the lattice node position. The
subscript, −𝑖, means the opposite direction of 𝑖. The purpose of the solid collision operator (Ω𝑠 ) is
to enforce a no-slip boundary condition between two phases of fluid and solid. This is achieved by
setting the PDF 𝑓𝑖 (𝐱 + 𝒄𝑖 𝛿𝑡 , 𝑡 + 𝛿𝑡 ) equal to the equilibrium distribution function (EDF),
𝑓𝑒𝑞
, 𝐮 , plus the bounce-back of the non-equilibrium part in the opposite direction (G. C. Yang
𝑖 (𝜌𝑓 𝑠 )
et al., 2019). The non-equilibrium bounce-back term ensures that the fluid particles reflect off the
solid surface with a velocity corresponding to the macroscopic velocity of the solid phase at that
lattice node.
In case the cell is partially saturated, namely 𝜀 is between 0 and 1; a weighting factor
incorporates both collision operators into account.
(2.24)

Ω = 𝐵Ω𝑠 + (1 − 𝐵)Ω 𝑓

where the weighting factor proposed by Noble and Torczynski (Noble & Torczynski, 1998) can
be calculated based on the relaxation time and the solid ratio as follows:

𝐵(𝜀, 𝜏) =

𝜀(𝜏 − 1⁄2)
(1 − 𝜀) + (𝜏 − 1⁄2)

(2.25)

Afterward, the following equations provide the hydrodynamic force and torque at each node
according to the lattice structure that was adopted earlier.
𝑛

18

𝐅𝑓 = ∑ 𝐵𝑗 ∑ Ω𝑖𝑠 𝐜𝑖
𝑗=1

(2.26)

𝑖=0

23

𝑛

18

𝐓𝑓 = ∑ [𝐵𝑗 (𝐱𝑗 − 𝐱 𝑠 ) × ∑ Ω𝑖𝑠 𝐜𝑖 ]
𝑗=1

(2.27)

𝑖=0

where 𝐱 𝑠 is the center of mass of solid particles and 𝐱𝑗 is the jth lattice cell coordinates. When the
hydrodynamic force and torque are calculated from Equations (2.26) and (2.27), then they are
incorporated again in Equations (2.21) and (2.22) to update solid particles’ velocity and location.
By utilizing these equations, the particle’s motion and position are adjusted based on the calculated
hydrodynamic forces and torques, ensuring the accurate representation of the particle’s dynamics
within the simulation.
2.5 LBM-DEM coupling scheme
The coupled LBM-DEM code (Seil & Pirker, 2017) used in this study utilizes Palabos (Latt et al.,
2021), an open-source C++ library for LBM, as the master program. This master program then
calls an external library, LIGGGHTS (Kloss, Goniva, Hager, Amberger, & Pirker, 2012), for DEM
simulations. The flowchart in Figure 2.5 illustrates the coupling scheme employed in this study.
The computational cycle begins by generating DEM particles and then initializing the fluid field.
Subsequently, interactions of particles are calculated using Equations (2.8) and (2.16).

24

Figure 2.5: Flowchart of LBM-DEM coupling method (G. C. Yang et al., 2019)
In order to ensure stability in the DEM simulation, it is crucial for the DEM time step (t) to
be smaller than a critical time step (𝑡𝑐𝑟 ) that can be obtained by the stiffness and mass of the
particles. Therefore, the calculated critical DEM time step (𝑡𝑐𝑟 ) is generally smaller than the time
step (t) used in LBM simulations, particularly for geotechnical engineering problems where soils
and rocks exhibit high stiffness (G. C. Yang et al., 2019). To synchronize the DEM simulations
with LBM, 𝑁𝑠𝑢𝑏 DEM sub-cycles are performed for each LBM evolution step, resulting in the
following relationship:

∆𝑡𝐷𝐸𝑀 =

𝛿𝑡
𝑁𝑠𝑢𝑏

(2.28)

25

which means that for every LBM timestep 𝛿𝑡 , there are 𝑁𝑠𝑢𝑏 DEM sub-cycles during which
hydrodynamic force and torque are constant.
In an LBM-DEM simulation, once DEM sub-cycles (𝑁𝑠𝑢𝑏 ) are completed, and the updated
positions of the particles are mapped onto the lattice grid. This mapping process identifies the
lattice cells that are covered by solid particles. For partially saturated cells, the solid ratio is
determined using the cell decomposition method, as depicted in Figure 2.3. Based on the type of
the cell (fluid, partially saturated, or solid), the collision process happens. It is worth noting that
the modified collision operator is also present in Equations (2.26) and (2.27), enabling the
calculation of hydrodynamic forces and torques immediately after the collision process, before the
streaming step, thus ensuring high computational efficiency. Subsequently, the hydrodynamic
forces are returned to the DEM to determine the motion of the particles. The resulting PDFs are
then streamed to the neighbouring lattice nodes. Using the redistributed PDFs, the fluid density
(𝜌𝑓 ) and velocity (𝐮𝑓 ) are updated according to Equations (2.3) and (2.4). If necessary, the fluid
pressure field can be obtained using Equation (2.6). After all the steps are taken, one cycle of the
LBM-DEM simulation concludes, and the simulation continues until it reaches the number of
cycles that is preset.
2.6 Irregular shape simulation
The practice of modelling real particles as clumps of overlapping spheres or clusters of nonoverlapping spheres has been widely adopted in many studies, particularly when particle
morphology is a key consideration. Simulating irregularly shaped particles has various
applications in agriculture (Pasha, Hare, Ghadiri, Gunadi, & Piccione, 2016), civil engineering
(Suhr & Six, 2020), chemical engineering (Nan, Ghadiri, & Wang, 2017) as well as geotechnical
engineering (Kodicherla, Gong, Fan, Wilkinson, & Moy, 2020). In this regard, there are a number
26

of methods proposed in the literature, including Wang et al. (L. Wang, Park, & Fu, 2007), Ferellec
and McDowell (Ferellec & McDowell, 2010), Favier et al. (Favier, Abbaspour‐Fard, Kremmer,
& Raji, 1999), and Taghavi (Taghavi, 2011), to name a few. The used open-source code,
LIGGGHTS, does not have functions for generating clumps. Therefore, this study adopted
CLUMP (Angelidakis, Nadimi, Otsubo, & Utili, 2021), a MATLAB code that has modules to
generate and export clumps based on Euclidean distance transform. The code has the capability to
even transform the self-generated clump into surface, that is essential if user would like to compare
the morphological characterisation of the generated volume to the original particle.
In this clump-generating code, the particle morphology undergoes a transformation from a
surface mesh into a voxelated representation. The flowchart of the code is shown in Figure 2.6.
The user can set a minimum dimension to control the number of voxels representing the particle.
Within the voxelated image, all voxels corresponding to the particle are assigned the value zero.
The Euclidean distance transform calculates the minimum distance between each zero voxel to its
nearest non-zero voxel. Afterward, the voxel that has the maximum value of Euclidean distance
transform is set as the first inscribed sphere and the corresponding voxels are set as one. After that,
another Euclidean distance transform calculates the new voxelated image. This process repeats
until one of the two conditions is met. The first condition happens when the number of spheres
passes the maximum user-defined number of spheres. The second condition occurs when the
algorithm wants to generate a sphere that is smaller than the user-defined minimum radius.
Moreover, this approach allows for generating overlapping spheres by selectively setting a certain
percentage of voxels to one, rather than all, forming each new sphere. The extent of overlap is
controlled using a variable overlap, which can take values within the range of [0, 1), providing
flexibility in managing the degree of overlap among the spheres produced. For this study, all
27

shapes are represented by 100 sub-particles, and the rest of the setting is brought in Table 2.2.
Figure 2.7 also provides the illustration of two particle geometries after transformation into a
clump.

Figure 2.6: CLUMP algorithm flowchart

28

Figure 2.7: Illustration of original STL file (left column), corresponding clump (middle column)
and the subsequent shape (right column) of particles after transformation into clump
Table 2.2: Parameters used in CLUMP algorithm for generating clumps
Minimum
dimension

Minimum
radius

Minimum No. of
spheres

Maximum
overlap

400

0

100

95%

The generated clumps consist of overlapping spheres, and it is important to correctly calculate
the center of mass and moment of inertia for every clump before the simulation starts. This
calculation occurs in DEM, and it involves the initial generation of an oriented bounding box for
a given multi-sphere clump (Figure 2.8), followed by the implementation of a Monte Carlo
simulation process (Y. Zhou, Zhou, Li, & Wang, 2017). Within this process, thousands of points
(denoted as N) were randomly distributed within the bounding box. Subsequently, points located
outside the clump were eliminated (red dots), while those within the clump were retained (black
dots). The number of points remaining inside the clump was denoted as n, and the ratio n/N was

29

considered as an approximate volume ratio representing the proportion of the clump within the
bounding box.

Figure 2.8: Monte Carlo algorithm for calculating the center of mass and moment of inertia in
clumps
The next step involved treating the retained points as mass-based voxels. Various kinematic
parameters, including mass center and inertia tensor, were then obtained by integrating these
voxels in accordance with their mathematical definitions (Y. Zhou et al., 2017). Notably, this
approach offers a distinct advantage by excluding the influence of overlapping density.
Consequently, it facilitates the extraction of true kinematic parameters for the clump, providing a
more accurate representation of its mass distribution and inertial properties (Y. Zhou et al., 2017).
This methodology is particularly advantageous in situations where obtaining precise kinematic
parameters without the interference of overlapping density is crucial for the analysis of clump
behaviour.

30

Chapter 3: Validation
3.1 Introduction
This chapter is dedicated to validating LBM-DEM coupling method against known benchmark
cases from the literature. Before this, the chapter starts with finding the right grid size, relaxation
time, and number of sub-particles that ensure accurate result by conducting a sensitivity study on
mesh size, relaxation time, and number of sub-particles needed to reconstruct an irregularly shaped
template. Then, the results for the three benchmark cases, which are 1) settling of a single clump,
2) settling of two clumps, and 3) settling of a disc-shaped particle, are provided. In all cases, the
particle’s shape is reconstructed through CLUMP as a volume of overlapping sub-particle spheres.
3.2 Mesh sensitivity, relaxation time, and sub-particle study
In any computational fluid dynamics simulation, ensuring that the grid size and time step do not
affect the results is crucial. In the LBM, this process occurs through investigating different mesh
resolutions and relaxation time, which is another representation of time step in LBM. The
simulation case chosen for this purpose is the sedimentation of a single sphere in a fluid as shown
in Figure 3.1. The sphere shape is reconstructed with 100 sub-particles. The diameter of the clump
is 15 mm, and the fluid density and viscosity are 1000 𝑘𝑔⁄𝑚3 , and 0.001𝑃𝑎. 𝑠, respectively. The
density of the particle is 1120 𝑘𝑔⁄𝑚3 .

31

Figure 3.1: Domain of sedimentation of a single clump. All boundaries are considered as no-slip
walls.
For the resolution study, Figure 3.2 demonstrates the effect of four different resolutions, 𝑁 =
10,15, 20, 25. It is worth noting that the number 𝑁 in each case represents the number of cells in
one diameter of sphere, which is 15 mm in this case. In order to depict the difference between each
grid size, Figure 3.2 only shows the last two seconds of the sedimentation, as the discrepancies in
this interval is more apparent. According to the figure, while there is an apparent discrepancy
between 𝑁 = 10 and the rest, 𝑁 = 20 is almost identical to 𝑁 = 25 meaning that 𝑁 = 20 would
be the right chose for the grid size resolution. Afterward, by setting N to 20, Figure 3.3 displays
the effect of relaxation time on the simulation outcome. The figure displays that as the value for
the relaxation factor gets closer to 0.5, the velocity trend gets smoother around the time the particle
reaches the bottom. Therefore, as a trade-off between accuracy and computational cost, in this
study we chose 0.56 or less for most of the simulation cases.

32

Figure 3.2: Study the effect of gird resolution on the simulation of sedimentation of a spherical
particle in a tank of fluid.

Figure 3.3: Study the effect of relaxation time on the simulation of sedimentation of a
spherical particle in a tank of fluid.
Moreover, to investigate the impact of the sub-particle number used in generating a clump,
Figure 3.4 illustrates the sedimentation velocity of the same irregularly shaped particle using
33

different number of sub-particles. When only 20 spheres are used to create the geometry of the
irregularly shaped particle, the terminal velocity is off by 0.005 m/s compared with other results.
However, as the number of particles increases, the terminal velocity converges to a certain value.
Therefore, in this study, 100 is chosen as the number of sub-particles needed to generate a clump.

Figure 3.4: Sedimentation velocity of an irregularly shaped particle with different numbers of
sub-particles in the clump. The actual shape of the particle is shown at the side.
3.3 Benchmark Cases
3.3.1 Settling of a single clump
For the first benchmark case, the settling of a heavy single particle in an ambient fluid is simulated.
The single spherical particle shown in Figure 2.7 is regenerated using 100 sub-particles for a
clump. The simulation domain is same as Figure 3.1. The diameter of the clump is 15 mm, and its
density is 1120 𝑘𝑔⁄𝑚3 . The fluid density and viscosity are according to (ten Cate, Nieuwstad,
Derksen, & Van den Akker, 2002).

34

For this simulation, the lattice resolution and relaxation times are 20 and 0.56, respectively.
By considering the diameter of the clump particle as the length scale and the particle’s terminal
velocity as the characteristic velocity, the simulation was performed for four Reynolds numbers
(Re), namely 1.48, 4.1, 11.6, and 31.9. The Reynolds number is calculated as follows:

𝑅𝑒 =

𝑈𝑑𝑝
𝜗

3.1)

where 𝑈 is the terminal velocity of the particle, 𝑑𝑝 is the diameter of the particle, and 𝜗 is the
kinetic viscosity. The results were compared to experimental data, as shown in Figure 3.5.

Figure 3.5: Comparison of current simulation and experimental results (ten Cate et al., 2002) of a
single clump settling in ambient fluid.
In Figure 3.5, for low Reynolds numbers, the velocity of the particle first increases until it
reaches its terminal velocity, where the buoyancy force is equal to gravity. Subsequently, as the
particle gets closer to the bottom, the velocity gradually decreases. For higher Reynolds numbers,
the same process is observed but the particle settles faster. Overall, Figure 3.5 demonstrates that

35

the sedimentation velocity of a spherical particle simulated by a clump is in good agreement with
the experimental results at each of the four Reynolds numbers (ten Cate et al., 2002).
Figure 3.6 presents the settling trajectory of the sphere at different Reynolds numbers.
According to the figure, lower Reynolds number takes much more time for the particle to reach
the bottom; and as the Reynolds number increases, the particle falls faster. Figure 3.6 once again
verifies that the code can accurately simulate both the velocity and the trajectory of the particle
with little discrepancies.

Figure 3.6: Comparing the simulation results with the experiment (ten Cate et al., 2002) for
settling trajectory of the sphere particle at various Reynolds numbers
Figure 3.7 provides an appreciation of the vertical velocity contour of the x-plane that is
located in the middle of the tank. The figure depicts that as the particle gains more speed, the wake
that occurs behind the particle also gets bigger. The wake region has lower pressure than the
surroundings. Therefore, as the wake region grows, the pressure difference between the front and
back of the particle increases, leading to a higher drag force. In this regard, the particle loses its
velocity as it approaches the bottom of the tank.
36

(a) t = 0.72 s

(b) t = 1.48 s

Figure 3.7: Vertical velocity contour of a clump at different times (Re = 11.6)
3.3.2 Settling of two clumps
In the second benchmark validation, the interaction between two clumps falling from a certain
height is investigated. Figure 3.8 shows the dimensions of the simulated domain. All boundaries
are no-slip wall. In this case, the density of particles is 1350 kg⁄m3 with the diameter equal to
12.7 mm. The Young modulus and passion ratio of each particle are set 5 MPa, and 0.4
respectively. For this simulation, the lattice resolution and relaxation times are 20 and 0.52,
respectively. In addition, fluid density and viscosity are 1195 kg⁄m3 and 0.0305
Pa. s, respectively. Figure 3.9 shows the simulation results in comparison with those by Dash and
Lee (Dash & Lee, 2015). The following equations are used to non-dimensionalize the time and
velocity in Figure 3.9.

𝑡𝑐 = √𝑑𝑝 /|𝜌𝑏 ⁄𝜌𝑓 − 1|𝑔

(3.2)

𝑈𝑐 = √|𝜌𝑏 ⁄𝜌𝑓 − 1|𝑔𝑑𝑝

(3.3)
37

Figure 3.8: Domain of the benchmark case settling of two clumps. All boundaries are considered
as no-slip walls.
As shown in Figure 3.9, when two particles start to fall (Figure 3.10a), they initially have the
same speed, but after some time, the trailing particle gets more speed since it is located at the wake
of the leading particle, also called drafting (Figure 3.10b). Then, the particles start to contact each
other (kissing, Figure 3.10c) until they begin to tumble and separate at a later stage (Figure 3.10d).
While there is a slight discrepancy between the two results, Figure 3.9 shows that our simulation
of clump particles can follow all three stages of drafting, kissing, and tumbling, which are
associated with this benchmark. Note that in this simulation, in order to capture the aforementioned
three stages, a small offset is introduced to the initial position of two clumps, in which the position
of leading clump is ( 5.5𝑑𝑝 + ∆𝑥, 5.5𝑑𝑝 + ∆𝑦, 35𝑑𝑝 ) and the position of trailing clump is
(5.5𝑑𝑝 , 5.5𝑑𝑝 , 37𝑑𝑝 ).

38

Figure 3.9: Settling velocity of two clumped particles in comparison with the results from (Dash
& Lee, 2015)

(a) t = 0 s

(b) t = 1.2 s

(c) t = 1.44 s

(d) t = 2.24 s

Figure 3.10: Sedimentation of two clumps at different times

39

3.3.3 Settling of a disc-shaped clump
The third simulated benchmark case is the settling of a disc-shaped clump. The shape of the disc
shown in Figure 2.7 is generated using a clump of 100 sub-particles. The lattice resolution and
relaxation times are 20 and 0.55, respectively. Furthermore, the initial orientation of the disc is set
to 45° as shown in Figure 3.11. According to the literature, the settling process of a disc-shaped
particle depends on the dimension of the disc. Therefore, we used this point to bring another
validation to our code.

Figure 3.11: Domain of the benchmark case disc shape clump. All boundaries are considered as
no-slip walls.
There are two non-dimensional numbers that characterize the settling process of a disc-shaped
particle, the Reynolds number and the non-dimensional moment of inertia, with the following
definition (Willmarth, Hawk, & Harvey, 1964).
𝐼∗ =

𝜋𝜌𝑝 𝑡
64𝜌𝑓 𝑑

(3.4)

In this equation, 𝑡, 𝑑, and 𝜌𝑝 are the thickness, diameter, and density of the disc, respectively.
The study presents that the falling of a disc in a fluid depends on the value of the two numbers (𝐼 ∗
40

and 𝑅𝑒) and could have either stable or unstable falling patterns. Therefore, we specifically set the
properties of the fluid and the disc to simulate two scenarios, including stable falling (case 1) and
periodic oscillation (case 2). Table 3.1 gives the characteristics of the domain for two cases. In
addition, Table 3.2 provides the value of the non-dimensional numbers for each case. It should be
mentioned that in both cases, the thickness of the disc is 0.3 times the diameter.

Table 3.1: Non-dimensional number values for two cases of disc-shaped particle benchmark.
Case 1
Case 2

𝐼∗
0.0191
0.0398

𝑈𝑧,𝑚𝑒𝑎𝑛 (𝑚/𝑠)
0.021
0.1087

𝑅𝑒𝑝
21
267

Table 3.2: Characteristics of two cases related to disc-shaped particle benchmark.
Case 1
Case 2

𝜌𝑝 (𝑘𝑔⁄𝑚3 )
1300
2700

d (mm)
1 mm
2 mm

Box size (mm)
6×6×15
12×12×48

Figure 3.12 shows the settling velocity versus time for the two cases, which correctly
simulated the falling pattern of the disc-shaped particle consistent with those predicted by
(Willmarth et al., 1964). Moreover, Figure 3.13 provides further illustration of a disc falling with
the oscillating pattern. As the figure shows, based on the non-dimensional number of the particle
in Case 1, the particle correctly follows a stable and non-oscillating settling pattern. However, in
Case 2, the particle has more moment of inertia compared to Case 1, plus a higher Reynolds
number, both of which set the particle unstable and with a certain oscillating frequency.

41

Figure 3.12: Settling velocity of disc-shaped particle in two different cases.

(a) t = 0.14 s

(b) t = 0.28 s

(c) t = 0.38 s

Figure 3.13: Z-velocity contour of settling disc-shaped clump at different times (Case 2).

42

Chapter 4: Immersed granular collapse
4.1 Introduction
In this chapter, a study is conducted to demonstrate the potential of LBM-DEM in simulating
a system of irregularly shaped particles. The case study is a pool of fluid with a column of real
shaped particles initially stacked on one side of the domain. As the simulation begins, the column
starts to collapse due to gravity force until it reaches its equilibrium state. This case study is
performed with four different particle shapes in addition to the one with perfectly sphere particles.
The chapter also looks into how morphological descriptors affect a particle’s drag force by
studying the sedimentation of each particle shape. It is worth mentioning that the computational
resource that is used in this study is a workstation with 16 GB RAM, and 11th Generation Intel
CoreTM i7-11700 @ 2.5 GHz × 16 of processor. All the simulations are performed using parallel
computing and 16 cores were used for each run.
4.2 Irregularly shaped particles
Four irregularly shaped particles were chosen to investigate in this study. The particles’ geometry
and their morphological descriptors are taken from a set of Jordan sand collections (Zheng et al.,
2019). The particles are chosen in a way that, based on the Zingg diagram (Figure 4.1) each one
represents certain shapes i.e. disc, sphere, blade, and rod. The morphological descriptors, including
elongation and flatness, are calculated according to Table 4.1. The real shape of the particles is
shown in Figure 4.2. In addition, the morphological descriptor values of the four considered
irregular particles are given in Table 4.2.

43

Figure 4.1: Zingg diagram. (a) Jordan No. 1, (b) Jordan No. 26, (c) Jordan No. 59, (d) Jordan No.
70
Table 4.1: 3D morphological descriptors (Zheng et al., 2021)
Descriptor

Symbol

Definition

Equivalent diameter

𝐷𝑒𝑞

Diameter of a sphere with grain volume 𝑉

Length

𝐿

Breadth

𝐵

Width

𝑊

Largest Feret diameter
Largest Feret diameter at planes normal to the 3D largest
Feret dimeter
Smallest Feret diameter

Elongation

𝐸𝑙

𝐵⁄𝐿

Flatness

𝐹𝑙

𝑊 ⁄𝐵

Sphericity

𝑆

2 ⁄
𝜋𝐷𝑒𝑞
𝐴 (𝐴 is the grain surface area)

44

No. 1

No. 26

No. 59

No. 70

Figure 4.2: Shape of Jordan sands used in this study
Table 4.2: Morphological descriptors of irregularly shaped particles
Jordan Sand No.
1
26
59
70

Flatness
0.94
0.53
0.53
0.80

Elongation
0.71
0.98
0.53
0.47

Sphericity
0.96
0.92
0.76
0.80

A key effect of particles with irregular shapes is observed in their influence on drag force,
where terminal velocity serves as the primary indicator of the drag force. The terminal velocity
signifies the point at which the forces of buoyancy, gravity, and drag neutralize each other,
resulting in the particle moving at a constant velocity thereafter. While the first two
aforementioned forces are only related to the density, the third one is where shape plays its role.
To illustrate this point, Figure 4.3 shows the terminal velocity of the four particles while all having
the same volume and density as 2460 𝑘𝑔⁄𝑚3. The figure points out that when the particle has a
more spherical shape, it experiences less drag force, i.e. higher terminal velocity. However, the
drag force is stronger for those particles with higher flatness, which eventually leads to less
terminal velocity. Figure 4.4 demonstrates that because particles No. 59, and 70 have a more flat
45

shape when they are falling down, they deviate from the central axis of the domain, in addition to
less sedimentation velocity compared to particles No. 1, and 26 which have a more round shape.
Moreover, the figure displays a stronger wake region behind particles No. 1, and 26 as they are
falling with faster velocity compared to the other two particles.

Figure 4.3: Sedimentation velocity of four different Jordan sand particles

46

(a)

(b)

(c)

(d)

Figure 4.4: Sedimentation of four irregularly shaped particles at t=0.12 s. (a) particle No. 1,
(b) particle No. 26, (c) particle No. 59, (d) particle No. 70.
4.3 Simulation of column collapse
To demonstrate the application of the current coupling method in geotechnical engineering, a
simulation of a column of irregular shape particles collapsing in an ambient fluid is presented. The
dimension of the domain is given in Figure 4.5. No slip wall is applied to all domain boundaries.
The modelling properties of the fluid and solid particles are given in Table 4.3. The simulations
are run for four different particle shapes. The morphological descriptions of each of the four
particles used in these simulations are given in Table 4.2 (Zheng et al., 2021). The lattice resolution
and relaxation times are set to 15 and 0.508, respectively. To generate the column of 1600 particles
(clumps), first, particles stacked at one side of the domain under gravity using only DEM
(𝐿𝑖 region). Then by removing the wall that holds the particles to one side, the particles start to
collapse to the other side of the domain.

47

Table 4.3: Characteristics of fluid and solid
Particle

Fluid

Parameter
Equivalent diameter of irregular shape particle
Density, 𝜌𝑝
Young’s modulus, E
Coefficient of restitution, 𝑒
Coefficient of friction, 𝜇
Poisson’s ratio
Density, 𝜌𝑓
Viscosity, 𝜇𝑓

Value
1.1 mm
2468 kg⁄m3
109 Pa
0.65
0.4
0.24
1000 kg⁄m3
0.001 Pa. s

Figure 4.5. Domain dimensions for immersed column collapse. All boundaries are no-slip walls.
The simulation result for particle No. 1 is shown in Figure 4.6. First, the particles near the
front top of the column started to slump down, and larger fluid velocities at t = 0.08 s were near
the front face (Figure 4.6a). As time passed, some of the slumping particles reached the floor and
created waves with high velocities near the toe (t=0.16 s). Some particles near the top front toppled
over and induced relatively large velocities in their vicinity. When t = 0.39 s, the rear and middle
portion of the column approximated their final slope grade, while the front toe had particles with
48

great momentum and kept progressing forward. Eventually at t = 0.52 s, the system reaches the
final state that is the kinetic energy of the system is less than 10−10 kJ.

(a)

(b)

(c)

(d)

Figure 4.6: Collapse of column of irregularly shaped particle No. 1 at different times a) t = 0.08
s, b) t = 0.16 s, c) t = 0.39 s, d) t = 0.52 s
Particle No. 26 has the disc shape geometry according to the Zingg diagram. Unlike particle
No. 1, at the early stage of collapse (Figure 4.7a), the disc shape particle starts to topple near the
front top same as sliding from the very bottom. This early dynamic results in the toppling segment
converging with the sliding particles at the lower region, inducing a pronounced and irregular
vortex formation around the free surface area (Figure 4.7b). This uneven vortex then leads into a
bump at t = 0.33 (Figure 4.7c). Eventually, at the final stage of collapse, some particles roll down
from upward and make the final free surface profile (Figure 4.7d). This sequence illustrates the

49

complex interplay between particle geometry, motion dynamics, and fluid interaction, leading to
distinctive morphological changes in the collapse structure.

(a)

(b)

(c)

(d)

Figure 4.7: Collapse of column of irregularly shaped particle No. 26 at different times a) t = 0.08
s, b) t = 0.16 s, c) t = 0.33 s, d) t = 0.44 s
Particle No. 59, is the blade-shaped particle. The first point observable about this case is that
due to the special shape of the particle, the initial height of the column is higher than other cases
(Figure 4.8a). As the simulation starts, the particles located at the front top topple, with none from
the bottom of the column sliding forward (Figure 4.8b). Since, in this case, the particles are falling
from a higher elevation, by the time the top particles arrive at the bottom, a heavy vortex is all over
the free surface (Figure 4.8c). Afterward, the system's kinetic energy dissipates through
interactions, and the collapse decelerates, culminating in a new equilibrium where the particles

50

settle into a densely packed arrangement, significantly altering the fluid's flow patterns around
them. It is also worth noting that while particle template No. 59 has the least terminal velocity
compared to the other three particles, according to Figure 4.3, due to its special shape, the particle
tends to oscillate when it gets closer to its terminal velocity. Therefore, oscillation could be the
reason this particle template induces more disturbance into the fluid field and takes a longer time
to reach its final equilibrium state.

(a)

(b)

(c)

(d)

Figure 4.8: Collapse of column of irregularly shaped particle No. 59 at different times a) t = 0.08
s, b) t = 0.16 s, c) t = 0.46 s, d) t = 0.62 s
The last simulation of irregularly shaped particle is particle No. 70, which is the rod shape.
Looking at the first stage of this collapse illustrates how fast this particle shape reacted to gravity
force and deformed from its initial static equilibrium state (Figure 4.9a). Comparing Figure 4.9a,
and Figure 4.9b shows how the particles at the top are slumping down and forward, causing
51

disturbance in the fluid velocity. Another point in this process is that this particle template reaches
its final equilibrium state at a faster pace compared to the disc shape particle. The reason could be
the way column collapses, induces the least amount of disturbance in the fluid field. Therefore, the
whole process happens smoothly and in less than half a second. It is also worth noticing that
according to Figure 4.9, particle No. 70 has one of the lowest terminal velocity magnitudes
compared to other particle templates. Thus, it is reasonable to reach the equilibrium states within
less time relative to other particle templates.

(a)

(b)

(c)

(d)

Figure 4.9: Collapse of column of irregularly shaped particle No. 70 at different times a) t = 0.08
s, b) t = 0.16 s, c) t = 0.33 s, d) t = 0.44 s
The last simulation belongs to the column collapse of the sphere particles. This simulation is
provided mainly as a comparison between a regular sphere particle and other irregularly shaped
particles. As shown in Figure 4.10, the sphere shape reaches its final state faster than any other
52

particle templates. The reason is mainly due to the round and symmetric shape of the particle. In
other words, the round shape potentially reduces the interparticle friction, leading to particles
moving more freely in the domain with less interaction from surroundings. This can also be
observed in Figure 4.10c and Figure 4.10d, where those particles that are near the free surface slide
down and make their way to the bottom of the fluid tank with no so much friction from the
surrounding fluid.

(a)

(b)

(c)

(d)

Figure 4.10: Collapse of column of spherical particle at different times a) t = 0.08 s, b) t = 0.16 s,
c) t = 0.25 s, d) t = 0.33 s

53

Figure 4.11: Final deposition of column collapse of spherical particles. 𝐿𝑓 = 3.45053 𝑐𝑚
Figure 4.11 and Figure 4.12 show the final deposition state of the spherical and irregularly
shaped particles. 𝐿𝑓 is the final length of the granular column. By looking at the figures, it is
apparent that all particles, regardless of their shape have eventually almost the same free-surface
profile. However, in terms of 𝐿𝑓 , particle No. 59 has the longest deposition length with 𝐿𝑓 =
3.78061 𝑐𝑚. This is interesting in particular because the length of the deposition of particle No.
59 is even higher than that of a sphere. We believe that the main reason, as mentioned earlier, is
that particle No. 59, due to its shape, builds a stack of higher height; therefore, those particles that
are at the top have higher potential energy that later on converts into kinetic energy. In other words,
the higher level of the column contributes to the higher kinetic energy of the whole system as well
as reaching an equilibrium state at a slower rate. All in all, even though in all cases, the final
deposition profile of particles is roughly the same, the amount of time for each column of particle
as well as the final length of deposition are not the same among the studied case. This observation
once again demonstrates how the particle’s size is important and should be considered in
simulating the geohazard processes in geotechnical engineering.

54

(a)

(b)

(c)

(d)

Figure 4.12: Final deposition of column collapse of a) irregularly shaped particle No. 1 with
𝐿𝑓 = 3.24631 𝑐𝑚, b) irregularly shaped particle No. 26 with 𝐿𝑓 = 3.45091 𝑐𝑚, c) irregularly
shaped particle No. 59 with 𝐿𝑓 = 3.78061 𝑐𝑚, d) irregularly shaped particle No. 70 with 𝐿𝑓 =
3.31487 𝑐𝑚

55

Chapter 5: Conclusion
5.1 Contribution
In this study, an LBM-DEM coupling framework was established to simulate the fluid-solid
interaction considering irregular solid particle shapes. Irregular particles were represented as
clumps using an open MATLAB code and were used in two open-source codes, LIGGGHTS and
Palabos, for coupled LBM-DEM simulations. Three benchmarks, including the settling of a single
clump, the settling of two clumps, and the settling of a disc-shaped clump, were successfully
simulated, and the results are consistent with those in the literature. Notably, the simulated settling
of two clumps was satisfactory in duplicating the interaction of two clumps in the fluid domain,
including drafting, kissing, and tumbling. In addition, the simulation of a disc falling in the fluid
domain can adequately produce stable falling and periodic oscillations that were consistent with
the theoretical predictions. The number of sub-particles equal to 100 for representing an irregular
particle as a clump in the established framework seems to be sufficient in capturing the grain
morphology and their influence on the fluid-solid interaction.
The established methodology was employed to simulate the collapse of a column of
irregularly shaped particles with four different real shape particles. The simulation results can
reasonably reproduce the slumping and trending of granular particles and the changes in flow
velocity in the fluid domain. The simulation results were compared with the column collapse of
spherical particles. It was found that spherical particles possessed larger momentum, travelled
faster, and reached a longer run-out distance, which warrants the consideration of realistic particle
shape in granular flow when interacting with fluids. To the best of the authors’ knowledge, this
study represents the first time coupling LIGGGHTS and Palabos to simulate the interaction of the
fluid with irregularly shaped solid particles. The established methodology is ready to be used for
56

many other geotechnical engineering problems, such as internal erosion in granular packs where
particle shape has an important influence on pore structure and the related fine migration.
5.2 Limitations
Like every study conducted in the science, this study is also associated with some limitations.
While these limitations do not affect the integrity of the current results, addressing them would be
a chance to develop the computation model for a better and more accurate simulation outcome in
the future. The limitations are as follows:
•

In LBM method, this study adopted the single relaxation time BGK collision operator, which
is known to have non-physical viscosity dependence on boundary location mapping. Therefore,
two-relaxation time collision operator could enhance the accuracy of fluid flow simulation part
especially if the model is used to simulate flow in porous media.

•

In Chapter 4, for the sake of reducing computational cost, the grid resolution used for simulating
the case of immersed granular column collapse is N=15, while based on the grid study, N=20
is the resolution that provides the most accurate results.

•

Also in Chapter 4, and again for the sake of reducing the high computational cost, for the disc
and sphere shape particles, 50 sub-particles used to reconstruct the geometry, while based on
Figure 3.4, the 100 sub-particles is considered as the optimum number of sub-particles.

5.3 Future work
Based on the results achieved in this research, the current code is capable of performing the
following simulations:

57

•

Simulating the drag force of real irregularly shaped particles at various Reynolds numbers and
studying a correlation formula that could relate the drag force to geometric descriptors (such as
elongation, flatness, sphericity) and predict the drag force of irregularly shaped particles.

•

Simulating the suffusion phenomenon using irregularly shaped particles as fixed, immovable
coarse particles along with fine spherical particles that are able to move through the pores
according to the applied pressure gradient.

58

References
Aidun, C. K., & Clausen, J. R. (2010). Lattice-Boltzmann Method for Complex Flows. Annual
Review of Fluid Mechanics, 42(1), 439-472. doi:10.1146/annurev-fluid-121108-145519
Almeida, E., Spogis, N., & Silva, M. (2016). Computational study of the pneumatic separation of
sugarcane bagasse. Paper presented at the Proceedings of the 6th International
Conference on Engineering for Waste and Biomass Valorisation, Albi, France.
Angelidakis, V., Nadimi, S., Otsubo, M., & Utili, S. (2021). CLUMP: A Code Library to
generate Universal Multi-sphere Particles. SoftwareX, 15, 100735.
doi:https://doi.org/10.1016/j.softx.2021.100735
Bhatnagar, P. L., Gross, E. P., & Krook, M. (1954). A model for collision processes in gases. I.
Small amplitude processes in charged and neutral one-component systems. Physical
review, 94(3), 511.
Chen, F., Drumm, E. C., & Guiochon, G. (2011). Coupled discrete element and finite volume
solution of two classical soil mechanics problems. Computers and Geotechnics, 38(5),
638-647. doi:https://doi.org/10.1016/j.compgeo.2011.03.009
Chen, S., & Doolen, G. D. (1998). Lattice Boltzmann method for fluid flows. Annual Review of
Fluid Mechanics, 30(1), 329-364. doi:10.1146/annurev.fluid.30.1.329
Cho, G.-C., Dodds, J., & Santamarina, J. C. (2007). Closure to “Particle Shape Effects on
Packing Density, Stiffness, and Strength: Natural and Crushed Sands” by Gye-Chun Cho,
Jake Dodds, and J. Carlos Santamarina. Journal of Geotechnical and Geoenvironmental
Engineering, 133(11), 1474-1474. doi:10.1061/(ASCE)1090-0241(2007)133:11(1474)
Chu, K. W., Wang, B., Yu, A. B., & Vince, A. (2009). CFD-DEM modelling of multiphase flow
in dense medium cyclones. Powder Technology, 193(3), 235-247.
doi:https://doi.org/10.1016/j.powtec.2009.03.015
Cook, B. K., Noble, D. R., & Williams, J. R. (2004). A direct simulation method for particle‐
fluid systems. Engineering Computations, 21(2/3/4), 151-168.
doi:10.1108/02644400410519721
Cortés, C., & Gil, A. (2007). Modeling the gas and particle flow inside cyclone separators.
Progress in Energy and Combustion Science, 33(5), 409-452.
doi:https://doi.org/10.1016/j.pecs.2007.02.001
Crowe, C. T., Schwarzkopf, J. D., Sommerfeld, M., & Tsuji, Y. (2011). Multiphase flows with
droplets and particles: CRC press.
Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies.
Géotechnique, 29(1), 47-65. doi:10.1680/geot.1979.29.1.47
59

Dash, S. M., & Lee, T. S. (2015). Two spheres sedimentation dynamics in a viscous liquid
column. Computers & Fluids, 123, 218-234.
doi:https://doi.org/10.1016/j.compfluid.2015.10.003
Di Renzo, A., & Di Maio, F. P. (2004). Comparison of contact-force models for the simulation
of collisions in DEM-based granular flow codes. Chemical Engineering Science, 59(3),
525-541. doi:https://doi.org/10.1016/j.ces.2003.09.037
Ding, W.-T., & Xu, W.-J. (2018). Study on the multiphase fluid-solid interaction in granular
materials based on an LBM-DEM coupled method. Powder Technology, 335, 301-314.
doi:https://doi.org/10.1016/j.powtec.2018.05.006
El-Emam, M. A., Shi, W., & Zhou, L. (2019). CFD-DEM simulation and optimization of gascyclone performance with realistic macroscopic particulate matter. Advanced Powder
Technology, 30(11), 2686-2702. doi:https://doi.org/10.1016/j.apt.2019.08.015
Favier, J. F., Abbaspour‐Fard, M. H., Kremmer, M., & Raji, A. O. (1999). Shape representation
of axi-symmetrical, non-spherical particles in discrete element simulation using multielement model particles. Engineering Computations, 16(4), 467-480.
Feng, Z.-G., & Michaelides, E. E. (2004). The immersed boundary-lattice Boltzmann method for
solving fluid–particles interaction problems. Journal of Computational Physics, 195(2),
602-628. doi:https://doi.org/10.1016/j.jcp.2003.10.013
Ferellec, J.-F., & McDowell, G. R. (2010). A method to model realistic particle shape and inertia
in DEM. Granular Matter, 12(5), 459-467. doi:10.1007/s10035-010-0205-8
Foster, M., Fell, R., & Spannagle, M. (2000). The statistics of embankment dam failures and
accidents. Canadian Geotechnical Journal, 37(5), 1000-1024. doi:10.1139/t00-030
Galindo-Torres, S. A. (2013). A coupled Discrete Element Lattice Boltzmann Method for the
simulation of fluid–solid interaction with particles of general shapes. Computer Methods
in Applied Mechanics and Engineering, 265, 107-119.
doi:https://doi.org/10.1016/j.cma.2013.06.004
Gardner, M. (2023). Toward a Complete Kinematic Description of Hydraulic Plucking of
Fractured Rock. Journal of Hydraulic Engineering, 149(7), 04023015.
doi:doi:10.1061/JHEND8.HYENG-13193
Gardner, M., & Sitar, N. (2019). Coupled three-dimensional discrete element-lattice Boltzmann
methods for fluid-solid interaction with polyhedral particles. International Journal for
Numerical and Analytical Methods in Geomechanics, 43(14), 2270-2287.
doi:https://doi.org/10.1002/nag.2972
Gingold, R. A., & Monaghan, J. J. (1977). Smoothed particle hydrodynamics: theory and
application to non-spherical stars. Monthly Notices of the Royal Astronomical Society,
181(3), 375-389. doi:10.1093/mnras/181.3.375
60

Han, K., Feng, Y., & Owen, D. (2007). Numerical simulations of irregular particle transport in
turbulent flows using coupled LBM-DEM. Computer Modeling in Engineering Sciences,
18(2), 87.
Han, Y., & Cundall, P. A. (2013). LBM–DEM modeling of fluid–solid interaction in porous
media. International Journal for Numerical and Analytical Methods in Geomechanics,
37(10), 1391-1407. doi:https://doi.org/10.1002/nag.2096
He, X., & Luo, L.-S. (1997). Lattice Boltzmann Model for the Incompressible Navier–Stokes
Equation. Journal of Statistical Physics, 88(3), 927-944.
doi:10.1023/B:JOSS.0000015179.12689.e4
Ibrahim, A., & Meguid, M. A. (2020). Coupled Flow Modelling in Geotechnical and Ground
Engineering: An Overview. International Journal of Geosynthetics and Ground
Engineering, 6(3), 39. doi:10.1007/s40891-020-00223-0
Jansen, F., & Harting, J. (2011). From bijels to Pickering emulsions: A lattice Boltzmann study.
Physical Review E, 83(4), 046707. doi:10.1103/PhysRevE.83.046707
Jiang, F., Liu, H., Chen, X., & Tsuji, T. (2022). A coupled LBM-DEM method for simulating the
multiphase fluid-solid interaction problem. Journal of Computational Physics, 454,
110963. doi:https://doi.org/10.1016/j.jcp.2022.110963
Jiang, F., Matsumura, K., Ohgi, J., & Chen, X. (2021). A GPU-accelerated fluid–structureinteraction solver developed by coupling finite element and lattice Boltzmann methods.
Computer Physics Communications, 259, 107661.
doi:https://doi.org/10.1016/j.cpc.2020.107661
Jing, L., Yang, G. C., Kwok, C. Y., & Sobral, Y. D. (2019). Flow regimes and dynamic
similarity of immersed granular collapse: A CFD-DEM investigation. Powder
Technology, 345, 532-543. doi:https://doi.org/10.1016/j.powtec.2019.01.029
Jones, B. D., & Williams, J. R. (2017). Fast computation of accurate sphere-cube intersection
volume. Engineering Computations, 34(4), 1204-1216. doi:10.1108/EC-02-2016-0052
Kloss, C., Goniva, C., Hager, A., Amberger, S., & Pirker, S. (2012). Models, algorithms and
validation for opensource DEM and CFD–DEM. Progress in Computational Fluid
Dynamics, an International Journal, 12(2-3), 140-152. doi:10.1504/PCFD.2012.047457
Kodicherla, S. P. K., Gong, G., Fan, L., Wilkinson, S., & Moy, C. K. S. (2020). Investigations of
the effects of particle morphology on granular material behaviors using a multi-sphere
approach. Journal of Rock Mechanics and Geotechnical Engineering, 12(6), 1301-1312.
doi:https://doi.org/10.1016/j.jrmge.2020.04.005
Koshizuka, S., & Oka, Y. (1996). Moving-Particle Semi-Implicit Method for Fragmentation of
Incompressible Fluid. Nuclear Science and Engineering, 123(3), 421-434.
doi:10.13182/NSE96-A24205
61

Ladd, A. J. C. (1994). Numerical simulations of particulate suspensions via a discretized
Boltzmann equation. Part 1. Theoretical foundation. Journal of Fluid Mechanics, 271,
285-309. doi:10.1017/S0022112094001771
Lai, Z., Zhao, J., Zhao, S., & Huang, L. (2023). Signed distance field enhanced fully resolved
CFD-DEM for simulation of granular flows involving multiphase fluids and irregularly
shaped particles. Computer Methods in Applied Mechanics and Engineering, 414,
116195. doi:https://doi.org/10.1016/j.cma.2023.116195
Latt, J., Malaspinas, O., Kontaxakis, D., Parmigiani, A., Lagrava, D., Brogi, F., Chopard, B.
(2021). Palabos: Parallel Lattice Boltzmann Solver. Computers & Mathematics with
Applications, 81, 334-350. doi:https://doi.org/10.1016/j.camwa.2020.03.022
Lominé, F., Scholtès, L., Sibille, L., & Poullain, P. (2013). Modeling of fluid–solid interaction in
granular media with coupled lattice Boltzmann/discrete element methods: application to
piping erosion. Numerical and Analytical Methods in Geomechanics, 37(6), 577-596.
doi:https://doi.org/10.1002/nag.1109
Mao, W., Wang, Y., Yang, P., Huang, Y., & Zheng, H. (2023). Dynamics of granular debris
flows against slit dams based on the CFD–DEM method: effect of grain size distribution
and ambient environments. Acta Geotechnica, 18(11), 5811-5838. doi:10.1007/s11440023-01944-y
Nan, W., Ghadiri, M., & Wang, Y. (2017). Analysis of powder rheometry of FT4: Effect of
particle shape. Chemical Engineering Science, 173, 374-383.
doi:https://doi.org/10.1016/j.ces.2017.08.004
Nielsen, A. W., Probst, T., Petersen, T. U., & Sumer, B. M. (2015). Sinking of armour layer
around a vertical cylinder exposed to waves and current. Coastal Engineering, 100, 5866. doi:https://doi.org/10.1016/j.coastaleng.2015.03.010
Noble, D. R., & Torczynski, J. R. (1998). A Lattice-Boltzmann Method for Partially Saturated
Computational Cells. International Journal of Modern Physics C, 09(08), 1189-1201.
doi:10.1142/S0129183198001084
Oschmann, T., Hold, J., & Kruggel-Emden, H. (2014). Numerical investigation of mixing and
orientation of non-spherical particles in a model type fluidized bed. Powder Technology,
258, 304-323. doi:https://doi.org/10.1016/j.powtec.2014.03.046
Owen, D. R. J., Leonardi, C. R., & Feng, Y. T. (2011). An efficient framework for fluid–
structure interaction using the lattice Boltzmann method and immersed moving
boundaries. International Journal of Numerical Methods in Engineering(1-5), 66-95.
doi:https://doi.org/10.1002/nme.2985
Pasha, M., Hare, C., Ghadiri, M., Gunadi, A., & Piccione, P. M. (2016). Effect of particle shape
on flow in discrete element method simulation of a rotary batch seed coater. Powder
Technology, 296, 29-36. doi:https://doi.org/10.1016/j.powtec.2015.10.055
62

Peng, Z., Doroodchi, E., Luo, C., & Moghtaderi, B. (2014). Influence of void fraction calculation
on fidelity of CFD-DEM simulation of gas-solid bubbling fluidized beds. 60(6), 20002018. doi:https://doi.org/10.1002/aic.14421
Pozorski, J., & Olejnik, M. (2023). Smoothed particle hydrodynamics modelling of multiphase
flows: an overview. Acta Mechanica. doi:10.1007/s00707-023-03763-4
Rettinger, C., & Rüde, U. (2022). An efficient four-way coupled lattice Boltzmann – discrete
element method for fully resolved simulations of particle-laden flows. Journal of
Computational Physics, 453, 110942. doi:https://doi.org/10.1016/j.jcp.2022.110942
Sayyari, M. J., Ahmadian, M. H., & Kim, K. C. (2022). Three-dimensional condensation in a
vertical channel filled with metal foam using a pseudo-potential lattice Boltzmann model.
International Journal of Thermal Sciences, 172, 107352.
doi:https://doi.org/10.1016/j.ijthermalsci.2021.107352
Seil, P., & Pirker, S. (2017). LBDEMcoupling: Open-Source Power for Fluid-Particle Systems.
Paper presented at the Proceedings of the 7th International Conference on Discrete
Element Methods, Singapore.
Shen, Z., Wang, G., Huang, D., & Jin, F. (2022). A resolved CFD-DEM coupling model for
modeling two-phase fluids interaction with irregularly shaped particles. Journal of
Computational Physics, 448, 110695. doi:https://doi.org/10.1016/j.jcp.2021.110695
Soundararajan, K. K. (2015). Multi-scale multiphase modelling of granular flows. University of
Cambridge
Stratford, K., Adhikari, R., Pagonabarraga, I., & Desplat, J. C. (2005). Lattice Boltzmann for
Binary Fluids with Suspended Colloids. Journal of Statistical Physics, 121(1), 163-178.
doi:10.1007/s10955-005-8411-1
Suhr, B., & Six, K. (2020). Simple particle shapes for DEM simulations of railway ballast:
influence of shape descriptors on packing behaviour. Granular Matter, 22, 1-17.
Sun, X., Sakai, M., & Yamada, Y. (2013). Three-dimensional simulation of a solid–liquid flow
by the DEM–SPH method. Journal of Computational Physics, 248, 147-176.
doi:https://doi.org/10.1016/j.jcp.2013.04.019
Taghavi, R. (2011). Automatic clump generation based on mid-surface. Paper presented at the
Proceedings, 2nd international FLAC/DEM symposium, Melbourne.
Tavarez, F. A., & Plesha, M. E. (2007). Discrete element method for modelling solid and
particulate materials. International Journal for Numerical Methods in Engineering, 70(4),
379-404. doi:https://doi.org/10.1002/nme.1881
ten Cate, A., Nieuwstad, C. H., Derksen, J. J., & Van den Akker, H. E. A. (2002). Particle
imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere
settling under gravity. Physics of Fluids, 14(11), 4012-4025. doi:10.1063/1.1512918
63

Teng, P., Johansson, F., & Hellström, J. G. I. (2023). Modelling erosion of a single rock block
using a coupled CFD-DEM approach. Journal of Rock Mechanics and Geotechnical
Engineering, 15(9), 2375-2387. doi:https://doi.org/10.1016/j.jrmge.2023.06.001
Tsuji, Y. (2007). Multi-scale modeling of dense phase gas–particle flow. Chemical Engineering
Science, 62(13), 3410-3418. doi:https://doi.org/10.1016/j.ces.2006.12.090
Verlet, L. (1967). Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties
of Lennard-Jones Molecules. Physical review, 159(1), 98-103.
doi:10.1103/PhysRev.159.98
Wachs, A. (2009). A DEM-DLM/FD method for direct numerical simulation of particulate
flows: Sedimentation of polygonal isometric particles in a Newtonian fluid with
collisions. Computers & Fluids, 38(8), 1608-1628.
doi:https://doi.org/10.1016/j.compfluid.2009.01.005
Wang, D., Leonardi, C. R., & Aminossadati, S. M. (2018). Improved coupling of time
integration and hydrodynamic interaction in particle suspensions using the lattice
Boltzmann and discrete element methods. Computers & Mathematics with Applications,
75(7), 2593-2606. doi:https://doi.org/10.1016/j.camwa.2018.01.002
Wang, L., Park, J.-Y., & Fu, Y. (2007). Representation of real particles for DEM simulation
using X-ray tomography. Construction and Building Materials, 21(2), 338-346.
doi:https://doi.org/10.1016/j.conbuildmat.2005.08.013
Wang, M., Feng, Y. T., Owen, D. R. J., & Qu, T. M. (2019). A novel algorithm of immersed
moving boundary scheme for fluid–particle interactions in DEM–LBM. Computer
Methods in Applied Mechanics and Engineering, 346, 109-125.
doi:https://doi.org/10.1016/j.cma.2018.12.001
Wang, M., Feng, Y. T., Qu, T. M., Tao, S., & Zhao, T. T. (2020). Instability and treatments of
the coupled discrete element and lattice Boltzmann method by the immersed moving
boundary scheme. International Journal for Numerical Methods in Engineering, 121(21),
4901-4919. doi:https://doi.org/10.1002/nme.6499
Wang, M., Feng, Y. T., Qu, T. M., & Zhao, T. T. (2021). A coupled polygonal DEM-LBM
technique based on an immersed boundary method and energy-conserving contact
algorithm. Powder Technology, 381, 101-109.
doi:https://doi.org/10.1016/j.powtec.2020.11.081
Wang, T., Zhang, F., & Zheng, W. (2023). Suffusion of Gap-Graded Soil with Realistically
Shaped Coarse Grains: A DEM–DFM Numerical Study. International Journal of
Geomechanics, 23(1), 04022247. doi:10.1061/(ASCE)GM.1943-5622.0002616
Willmarth, W. W., Hawk, N. E., & Harvey, R. L. (1964). Steady and Unsteady Motions and
Wakes of Freely Falling Disks. The Physics of Fluids, 7(2), 197-208.
doi:10.1063/1.1711133
64

Yang, G. C., Jing, L., Kwok, C. Y., & Sobral, Y. D. (2019). A comprehensive parametric study
of LBM-DEM for immersed granular flows. Computers and Geotechnics, 114, 103100.
doi:https://doi.org/10.1016/j.compgeo.2019.103100
Yang, L., Padding, J. T., & Kuipers, J. A. M. (2016). Modification of kinetic theory of granular
flow for frictional spheres, Part I: Two-fluid model derivation and numerical
implementation. Chemical Engineering Science, 152, 767-782.
doi:https://doi.org/10.1016/j.ces.2016.05.031
Zhang, H., Tan, Y., Shu, S., Niu, X., Trias, F. X., Yang, D., Sheng, Y. (2014). Numerical
investigation on the role of discrete element method in combined LBM–IBM–DEM
modeling. Computers & Fluids, 94, 37-48.
doi:https://doi.org/10.1016/j.compfluid.2014.01.032
Zhang, R., Kim, H.-J., & Dinoy, P. R. (2021). Particle Flow Simulation Based on Hybrid IMBDEM-LBM Approach with New Solid Fraction Calculation Scheme. Applied Sciences,
11(8). doi:10.3390/app11083436
Zheng, W., Hu, X., Tannant, D. D., Zhang, K., & Xu, C. (2019). Characterization of two- and
three-dimensional morphological properties of fragmented sand grains. Engineering
Geology, 263, 105358. doi:https://doi.org/10.1016/j.enggeo.2019.105358
Zheng, W., Hu, X., Tannant, D. D., & Zhou, B. (2021). Quantifying the influence of grain
morphology on sand hydraulic conductivity: A detailed pore-scale study. Computers and
Geotechnics, 135, 104147. doi:https://doi.org/10.1016/j.compgeo.2021.104147
Zhou, W., Ma, Q., Ma, G., Cao, X., & Cheng, Y. (2020). Microscopic investigation of internal
erosion in binary mixtures via the coupled LBM-DEM method. Powder Technology, 376,
31-41. doi:https://doi.org/10.1016/j.powtec.2020.07.099
Zhou, Y., Zhou, B., Li, J., & Wang, H. (2017). Study on the Multi-sphere Method Modeling the
3D Particle Morphology in DEM. Paper presented at the Proceedings of the 7th
International Conference on Discrete Element Methods, Singapore.
Zhou, Z. Y., Kuang, S. B., Chu, K. W., & Yu, A. B. (2010). Discrete particle simulation of
particle–fluid flow: model formulations and their applicability. Journal of Fluid
Mechanics, 661, 482-510. doi:10.1017/S002211201000306X
Zou, Q., & He, X. (1997). On pressure and velocity boundary conditions for the lattice
Boltzmann BGK model. Physics of Fluids, 9(6), 1591-1598. doi:10.1063/1.869307

65