Modelling Complex Particle–Fluid Flow with a Discrete Element Method Coupled with Lattice Boltzmann Methods (DEM-LBM)

Abstract

Particle–fluid flows are ubiquitous in nature and industry. Understanding the dynamic behaviour of these complex flows becomes a rapidly developing interdisciplinary research focus. In this work, a numerical modelling approach for complex particle–fluid flows using the discrete element method coupled with the lattice Boltzmann method (DEM-LBM) is presented. The discrete element method and the lattice Boltzmann method, as well as the coupling techniques, are discussed in detail. The DEM-LBM is thoroughly validated for typical benchmark cases: the single-phase Poiseuille flow, the gravitational settling and the drag force on a fixed particle. In order to demonstrate the potential and applicability of DEM-LBM, three case studies are performed, which include the inertial migration of dense particle suspensions, the agglomeration of adhesive particle flows in channel flow and the sedimentation of particles in cavity flow. It is shown that DEM-LBM is a robust numerical approach for analysing complex particle–fluid flows.

1. Introduction

Particulate flows are extensively encountered in nature and industrial processes, attracting tremendous engineering research interests in almost all areas of sciences [1,2,3], including astrophysics, chemical engineering, biology, life science and so on. Due to the complex particle–particle interactions and their interactions with the surrounding media, i.e., gas/liquid environment or electrostatic/magnetic field, the dynamic behaviour of particulate flow is very complicated. Therefore, it is of significant importance to understand the particle dynamics, in order to promote relevant manufacture and production processes.

With the rapid development of the modern modelling technique, abundant numerical studies on the particulate flows have sprung up in recent years. As the name suggests, the particulate flow can behave like a continuous fluid or a rigid solid. Therefore, studies of particulate flows can be classified into three categories based on the types of approach: (1) continuum modelling, which extends the continuum mechanics of single-phase fluid to describe the particle transportation, leading to a representative population balance method [4,5]; (2) developing the kinetic theory of particulate flow [6] based on the averaged equations for multi-particle systems, which generalise the dynamics of particle–particle collision processes; (3) discrete particle modelling, which solves the particle’s motion individually based on certain interaction laws. The third approach is also classified as Lagrangian particle method and has many different variations, including the discrete element method (DEM), dissipative particle dynamics (DPD), molecular dynamics (MD) and Brownian dynamics (BD). These discrete particle methods differentiate themselves with different particle–particle interaction laws and their corresponding length and time scales [3]. DEM is widely used for modelling of granular materials and adhesive particulate flows [7,8,9,10,11,12,13,14,15,16,17], and was first introduced by Cundall and Strack [18] to study the rock mechanics. In DEM, particles are treated as individual entities, and interact with each other through contact and non-contact forces. The translational and rotational motions of the particles are described using the Newton’s second law. The greatest advantage of DEM is that a large amount of dynamic information, which is almost unattainable experimentally, such as the transient forces and torques, can be obtained in great detail.

Apart from the particle–particle interaction, the hydrodynamic interaction associated with the surrounding fluid is also an important issue in the study of particulate flow. Various computational fluid dynamics (CFD) methods at different scales of time and length are developed to model the single-phase fluid flow, ranging from discrete models, e.g., the lattice Boltzmann method (LBM) [19,20,21], to continuum models, e.g., the direct numerical simulation (DNS) [22], the large eddy simulation (LES) [23] and the two-fluid model (TFM) [24]. Hence, a large number of hybrid models are developed for modelling particle–fluid flow [20,25,26,27,28,29,30,31,32,33]. A comprehensive comparison of various hybrid models was discussed in the literature [1,34].

Comparing with other conventional CFD approaches, LBM arises in the last two decades and has become a promising numerical method due to several advantages [19,20]. First, the primary computation procedure, i.e., the so-called collision operator, is performed locally with information only from where the computation takes place, facilitating straightforward implementation in a parallel way with high computational efficiency. Second, the LBM shows significant flexibility in handling complex boundary conditions, as they are treated with a local bounce-back fashion instead of the conventional boundary treatment. Third, the lattice grid adopted in LBM is traditionally orthogonal and smaller than the particle size, providing a high space resolution on the evaluation of the hydrodynamic interactions. As a result, the LBM was coupled with DEM and widely applied in analysing a large range of particle–fluid problems, such as particle suspensions [35,36,37,38], flow through porous media [39,40], multicomponent flows [20,41,42], particle-laden turbulent flows [31,43,44] and magnetohydrodynamics [45].

In this paper, a numerical approach to model complex particle–fluid flow using a coupled DEM-LBM is introduced. The fundamental equations of DEM and LBM, as well as various coupling techniques are presented in detail. The validity and accuracy of the numerical approach are demonstrated with some benchmark tests. Finally, case studies, including the migration, the agglomeration and the sedimentation of particle suspensions, are presented to illustrate the potential and applicability of the numerical approach. It should be noted that the interactions between particles and the surrounding environment, such as electrostatic interactions, liquid bridge effect and so on, are not within the scope of this work, but can be found in the literature [46,47,48,49,50,51,52,53,54].

2. Discrete Element Method (DEM)

In DEM, particle motion is solved individually using the Newton’s equation of motion,

$$\begin{array}{l}m\frac{d{U}_p}{dt}=G+{\displaystyle \sum F},\\ I\frac{d{\mathbf{\Omega }}_p}{dt}={\displaystyle \sum M},\end{array}$$

(1)

where U_p and Ω_p are the translational and rotational velocities of a particle, respectively. m and I are the mass and rotational inertia of the particle. G is the gravitational force. F and M represent the force and torque acting on the particle, which could come from various sources, such as interparticle collision, hydrodynamic interactions with the fluid, electrostatic interactions with charge and so on. As DEM was designed to deal with the numerous collision problems between the distinct particles, the collision forces and torques are the most important aspect. When two particles get in contact with each other, the collision force and torque can be decomposed as

$$\begin{array}{l}{\mathit{F}}_{collision}=F_n\mathbf{n}+F_s{\mathbf{t}}_s,\\ {\mathit{M}}_{collision}=r_p{F}_s(\mathbf{n}\times {\mathbf{t}}_s)+M_r({\mathbf{t}}_r\times \mathbf{n})+M_t\mathbf{n},\end{array}$$

(2)

where F_n and F_s denote the contact forces in the normal and tangential direction, respectively. M_r and M_t are the rolling and twisting resistance torques, respectively. r_p represents the particle radius n, t_s and t_r are the unit vectors in the normal, tangential and rolling direction at the contact point, respectively.

2.1. Normal Contact Force

The most common model for the normal contact force is the Hertz model [55], which can be expressed by a non-linear Hook’s law. The Hertz model was validated and widely applied in the contact of spherical particles with size above millimeters. However, as the particle approaches the micrometer size range, the adhesion effect due to the van der Waals force must be taken into consideration. Extensions of the Hertz model, including the well-known Johnson–Kendall–Roberts (JKR) theory [56] and Derjaguin–Muller–Toporov (DMT) theory [57], were also proposed to describe the contact force in presence of adhesion.

Let us consider two contacting particles at centroid positions x_i and x_j with radii r_p,i and r_p,j, respectively. The normal overlap δ_N between these two particles is derived as δ_N = r_p,i + r_p,j − |x_j− x_i|. Then the normal forces between two particles are given as:

$$\begin{array}{l}{F}_{ne-Hertz}=-k_N{\delta }_N^{1.5},\\ F_{ne-JKR}=-4F_C({\widehat{a}}^3-{\widehat{a}}^{1.5}),\end{array}$$

(3)

where k_N is normal stiffness coefficient in the Hertz model, $\widehat{a}=a/a_0$ is the normalized contact area radius and F_C represents the critical pull-off force obtained from the JKR theory. In the Hertz model, the normal stiffness is analytically derived as $k_N=4/3E\sqrt{R}$, where R is the effective particle radius and E is the effective elastic moduli, defined as:

$$\begin{array}{l}\frac{1}{R}=\frac{1}{r_{p,i}}+\frac{1}{r_{p,j}},\\ \frac{1}{E}=\frac{1-{\sigma }_i^2}{E_i}+\frac{1-{\sigma }_j^2}{E_j},\end{array}$$

(4)

with E_i and E_j being elastic moduli and σ_i and σ_j being the Poisson’s ratios, respectively. Note that the contact area radius a, the effective radius R and the normal overlap δ_N are related to each other via a simple geometrical relation, $a=\sqrt{R{\delta }_N}$. In the JKR theory, the particles will keep in contact due to the presence of van der Waals adhesion even if the external load is zero [56], which gives rise to an equilibrium contact radius a₀,

$$a_0={(\frac{9\pi \gamma R^2}{E})}^{1/3},$$

(5)

where γ is the particle’s surface energy. The critical pull-off force in the JKR theory is then derived as F_C = 3πγR [56]. However, both the Hertz model and the JKR theory are energy conserved, i.e., there is no energy dissipation during the process of particle collision. To consider the energy dissipation during a collision, various dissipation mechanisms were proposed, including the viscoelastic response, plastic deformation as well as the fluid-phase dissipation [58]. For the low-speed impact regime, the viscoelastic dissipation is dominant. The dissipation force F_nd is approximated with a dashpot model, which is proportional to the approaching velocity of the particles, i.e.,

$$F_{nd}=-{\eta }_N{\mathbf{v}}_R·\mathbf{n},$$

(6)

where η_N is the normal dissipation coefficient, and v_R is the relative velocity at the contact point.

2.2. Sliding, Rolling and Twisting Resistance

The tangential contact force results from the relative sliding motion between two particles. Similar to the normal force, the standard sliding model is approximated by a linear spring-dashpot,

$$F_s=-k_T{\xi }_T - {\eta }_T{v}_s·t_s,$$

(7)

where k_T is the spring stiffness coefficient in the tangential direction, v_s is relative sliding velocity, ${\xi }_T={\displaystyle {\int }_{t_0}^t{v}_s(t)·t_s dt}$ is the cumulative displacement in the tangential direction and η_T is the sliding dissipation coefficient. According to the Amonton’s friction law, two contacting particles start to slide relative to each other when the tangential force reaches a critical value F_s,crit and the dynamic friction force remains constant as F_s = −F_s,crit. For non-adhesive particles, the critical value is F_s,crit = μ_s |F_ne|, where μ_s is sliding friction coefficient. Due to the presence of adhesion, the critical value for adhesive particles is larger, F_s,crit = μ_s |F_ne + 2F_C| [3,58].

The rolling and twisting resistances are similarly postulated in the form:

$$\begin{array}{l}{M}_r=-k_R{\xi }_R-{\eta }_R{v}_L\cdot t_r,\\ M_t=-k_Q{\xi }_Q - {\eta }_Q{\mathrm{\Omega }}_T,\end{array}$$

(8)

where ${\xi }_R={\displaystyle {\int }_{t_0}^t{v}_L(t)\cdot t_r dt}$ and ${\xi }_Q={\displaystyle {\int }_{t_0}^t{\mathrm{\Omega }}_T(t)dt}$ are the rolling and twisting displacements, v_L= −R(Ω_p,i − Ω_p,j) × n and Ω_T= (Ω_p,i − Ω_p,j)· n are the relative rolling and twisting velocity, respectively. The rolling direction is t_r = v_L/|v_L|. k_R, k_Q are the rolling and twisting stiffness, respectively, and η_R, η_Q are the rolling and twisting dissipation coefficients, respectively. Correspondingly, the critical values for rolling and twisting resistances are given as:

$$\begin{array}{l}{M}_{r,crit}=k_R{\theta }_{crit}R,\\ M_{t,crit}=3\pi a{F}_{s,crit}/16,\end{array}$$

(9)

where θ_crit is the critical rolling angle in the relative rolling motion between two contacting particles.

3. Lattice Boltzmann Method (LBM)

3.1. Lattice Boltzmann Equation

Fundamentally, the LBM aims to establish a microscopic kinetic model involving essential physics, so as to make the average physical quantities in macroscopic scale follow the expected equations, which is completely different from the conventional numerical methods that directly discretise the macroscopic continuum equations. Hence an important premise lies in that the dynamic behaviour of a fluid in the macroscopic scale is the outcome from the collective behavior of numerous constitutions in the microscopic scale [19].

The LBM originated from lattice gas automata, which utilises a discrete lattice and discrete time. In each lattice, the fluid is described by a packet of microscopic particles, which only move in certain directions with certain discrete velocities. A number of lattice speed models are available, depending on the dimension of the problem. For example, the D2Q9 and D3Q19 are frequently used in the 2D and 3D problems, respectively, as illustrated in Figure 1. The individual particle motion is neglected in the LBM and a single-particle distribution function f_i(x,t) is adopted instead.

The evolution of the particle distribution functions f_i(x,t) follows the lattice Boltzmann equation (LBE) [19] with an external force term,

$$f_i(x+e_i\Delta t,t+\Delta t)=f_i(x,t)+{\mathrm{\Omega }}_i[f_i(x,t)]+F_i\Delta t,$$

(10)

where x is the position of the node, e_i is the unit vector of the lattice speed in the direction i as denoted in Figure 1, Δt is the explicit time step, F_i is the external body force and Ω_i[f_i(x,t)] is the collision operator. Equation (10) can be evaluated in two independent processes: collision and streaming, as shown below:

$$f_i(x,t+\Delta t)=f_i(x,t)+{\mathrm{\Omega }}_i[f_i(x,t)]+F_i\Delta t,$$

(11)

$$f_i(x+e_i\Delta t,t+\Delta t)=f_i(x,t+\Delta t).$$

(12)

Equation (11) represents the collision process, which updates the distribution functions at position x from t to t+Δt based on the relaxation operator Ω_i[f_i(x,t)]. Then the streaming process (Equation (12)) takes place, which simply propagates the updated distribution functions from position x to their nearest neighbor nodes x + e_iΔt according to the lattice speed model. The sequence of collision and streaming becomes irrelevant after a few time steps.

3.1.1. Single-Relaxation-Time Model

The single-relaxation-time (SRT) model refers to the Bhatnagar-Gross-Krook (BGK) approximation of the collision operator [59]. The collision operator is simplified to a linear form [59,60,61,62]:

$${\mathrm{\Omega }}_i=-\frac{\Delta t}{\tau }[f_i(x,t)-f_i^{eq}(x,t)],$$

(13)

where τ is the dimensionless relaxation parameter and f_i^eq(x,t) is the equilibrium distribution function that is defined as:

$$f_i^{eq}=\rho {\omega }_i[1+\frac{e_i\cdot u}{c_s^2}+\frac{{(e_i\cdot u)}^2}{2c_s^4}-\frac{u^2}{2c_s^2}].$$

(14)

In this equation, the weight coefficient ω_i depends on the lattice speed model. For example, ω₀ = 4/9, ω_1,2,3,4 = 1/9 and ω_5,6,7,8 = 1/36 for D2Q9 model, while ω₀ = 1/3, ω_1,…,6 = 1/18 and ω_7,…,18 = 1/36 for D3Q19 model. $c_s=1/\sqrt{3}$ is the lattice sound speed.

The discrete body force term F_i in Equation (10) is given as [63]:

$$F_i=(1-\frac{1}{2\tau }){\omega }_i[\frac{e_i-u}{c_s^2}+\frac{(e_i\cdot u)}{c_s^4}{e}_i]\cdot F,$$

(15)

where F represents the macroscopic body force. The macroscopic fluid properties, including the density ρ, velocity u and pressure p are related to the microscopic particle distribution function, which are determined as:

$$\begin{array}{l}\rho ={\displaystyle \sum _i{f}_i},\\ \rho u={\displaystyle \sum _i{f}_i{e}_i+\frac{\Delta t}{2}F},\\ p=c_s^2\rho .\end{array}$$

(16)

The kinematic viscosity of the fluid ν is related to the relaxation parameter as:

$$\nu =\frac{2\tau -1}{6}.$$

(17)

3.1.2. Multi-Relaxation-Time Model

It was reported that the SRT LBE does not show good stability with small relaxation time τ or with high Reynolds number. To improve the performance of the LBM, the multi-relaxation-time (MRT) LBE was proposed, which can be expressed as [20]:

$$f_i(x+e_i\Delta t,t+\Delta t)=f_i(x,t)-{\displaystyle \sum _{\alpha }{\mathbf{\Omega }}_{i\alpha }[f_{\alpha }(x,t)-f_{\alpha }^{eq}(x,t)]}+[S_i(x,t)-0.5{\displaystyle \sum _{\alpha }{\mathbf{\Omega }}_{i\alpha }{S}_{\alpha }(x,t)}],$$

(18)

where S is the force term. From Equation (15), we have:

$$S_i={\omega }_i[\frac{e_i-u}{c_s^2}+\frac{(e_i\cdot u)}{c_s^4}{e}_i]\cdot F.$$

(19)

Equation (18) can be transformed into the matrix form as:

$$f(x+e_i\Delta t,t+\Delta t)=f(x,t)-M^{-1}\Lambda [m(x,t)-m^{eq}(x,t)]+M^{-1}[{\rm I}-\frac{\Lambda }{2}]\overline{S}(x,t),$$

(20)

where M is the transformation matrix, m and m^eq are the moment spaces of the distribution function f_i and its equilibrium distribution f_i^eq, which can be derived as m = Mf_i and m^eq= Mf_i^eq, respectively. The collision matrix Λ is a diagonal matrix in the moment space.

The matrixes M and Λ vary for different speed models. For D2Q9 model, the matrix M is:

$$M=\left[\begin{array}{ccccccccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill -4\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill 2\hfill & \hfill 2\hfill & \hfill 2\hfill \\ \hfill 4\hfill & \hfill -2\hfill & \hfill -2\hfill & \hfill -2\hfill & \hfill -2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill -2\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -2\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill -1\hfill \end{array}\right].$$

(21)

The equilibrium distribution functions m^eq in the moment space are related to that in the velocity space as follows:

$$m^{eq}=M{f}_i^{eq}=\left[\begin{array}{c}\rho \\ e^{eq}\\ {\epsilon }^{eq}\\ j_x^{eq}\\ q_x^{eq}\\ j_y^{eq}\\ q_y^{eq}\\ p_{xx}^{eq}\\ p_{xy}^{eq}\end{array}\right]=\left[\begin{array}{c}\rho \\ -2\rho +3(u_x^2+u_y^2)\\ \rho -3(u_x^2+u_y^2)\\ \rho u_x\\ -u_x\\ \rho u_y\\ -u_y\\ u_x^2-u_y^2\\ u_x{u}_y\end{array}\right].$$

(22)

The diagonal matrix Λ is Λ = diag(s₁, s₂, s₃, s₄, s₅, s₆, s₇, s₈, s₉). In order to ensure the mass and momentum conservation after collision, it is easily derived that s₁ = s₄ = s₆ = 0. s₈ and s₉ are set as s₈ = s₉ = 1/τ in order to reproduce the same viscosity of SRT model [64]. The other relaxation parameters s₂, s₃, s₅ and s₇ can be decided with more flexibility according to the physical model.

For D3Q19 model, the matrix M is:

$$M=\left[\begin{array}{ccccccccccccccccccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill -30\hfill & \hfill -11\hfill & \hfill -11\hfill & \hfill -11\hfill & \hfill -11\hfill & \hfill -11\hfill & \hfill -11\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill & \hfill 8\hfill \\ \hfill 12\hfill & \hfill -4\hfill & \hfill -4\hfill & \hfill -4\hfill & \hfill -4\hfill & \hfill -4\hfill & \hfill -4\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -4\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 4\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -4\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 4\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -4\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 4\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -2\hfill \\ \hfill 0\hfill & \hfill -4\hfill & \hfill 2\hfill & \hfill 2\hfill & \hfill -4\hfill & \hfill 2\hfill & \hfill 2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -2\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -2\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -2\hfill & \hfill 2\hfill & \hfill 0\hfill & \hfill -2\hfill & \hfill 2\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill \end{array}\right].$$

(23)

The equilibrium distribution functions m^eq in the moment space is:

$$m^{eq}=M{f}_i^{eq}=\left[\begin{array}{c}\rho \\ -11\rho +19\rho (u_x^2+u_y^2+u_z^2)\\ 3\rho -\frac{11}{2}\rho (u_x^2+u_y^2+u_z^2)\\ \rho u_x\\ -\frac{2}{3}\rho u_x\\ \rho u_y\\ -\frac{2}{3}\rho u_y\\ \rho u_z\\ -\frac{2}{3}\rho u_z\\ \rho (2u_x^2-u_y^2-u_z^2)\\ -\frac{1}{2}\rho (2u_x^2-u_y^2-u_z^2)\\ \rho (u_y^2-u_z^2)\\ -\frac{1}{2}\rho (u_y^2-u_z^2)\\ \rho u_x{u}_y\\ \rho u_y{u}_z\\ \rho u_x{u}_z\\ 0\\ 0\\ 0\end{array}\right].$$

(24)

The diagonal matrix Λ is Λ = diag(s₁, s₂,……, s₁₈, s₁₉), where s₁ = s₄ = s₆ = s₈ = 1.0, s₂ = s₅ = s₇ = s₉ = 1.2, s₃ = s₁₁ = s₁₃ = 1.4, s₁₀ = s₁₂ = s₁₄ = s₁₅ = s₁₆ = 1/τ, s₁₇ = s₁₈ = s₁₉ = 1.98.

3.2. Boundary Conditions

In the LBM, the conventional velocity and pressure boundary conditions cannot be applied directly, since these quantities are not primary variables in the LBE. Instead, the boundary conditions are imposed with the distribution functions. In this section, only the ‘no-slip’ wall and the periodic boundary conditions are introduced. It should be noted that particles are treated as moving boundaries, which play important roles in the coupling between LBM and DEM. The interactions between the moving particles and the fluid will be discussed separately later.

The ‘no-slip’ boundary condition is implemented with the so-called bounce-back rule [65,66]. As shown in Figure 2, if f_i is a distribution function entering into the solid stationary wall, the undetermined distribution function that comes out of the wall is simply a reflection of f_i, i.e.,

$$f_{-i}(x,t+\Delta t)=f_i^{+}(x,t),$$

(25)

where f_−i and f_i⁺ denote the distribution functions in the opposite direction of i and after collision, respectively. This simple bounce-back rule ensures a strict ‘no-slip’ condition at the wall position, as the tangential velocity is zero. It should be noted that the accuracy of the simple bounce-back rule is only first order, except that the wall is just in the middle of the link between a solid node and a fluid node, where the accuracy is second order. However, the bounce-back rule is easy to implement and can be extended to obstacles with arbitrary shape, which makes the LBM quite suitable for problems with complex geometries.

Periodic boundary conditions are also straightforward to impose in the LBM. As shown in Figure 2, distribution functions propagating out of the domain stream into the corresponding boundary node at the opposite boundary in the same direction. Take Figure 2 as an example, the inlet and outlet planes are periodic. As a result, the distribution functions f₁, f₅ and f₈ at the outlet boundary will stream into the corresponding node at the inlet, while f₃, f₆ and f₇ at the inlet propagate into the corresponding node at the outlet.

3.3. Unit Conversion in LBM

In the LBM, the physical parameters in the computation are commonly dimensionless, which necessitates a unit conversion. In principle, one only needs to determine the conversion coefficients of three fundamental units, i.e., length, time and mass, and then the conversion coefficients of all the other physical parameters can be sequentially decided.

Normally, the lattice space Δx, the time step Δt, and the fluid density ρ are set to be one in the LBM computation. How to perform the unit conversion will be explained using an example of pipe flow in 3D. Suppose that the pipe size is L × D = 8 mm × 4 mm, where L is the pipe length and D is the pipe diameter, and water flows through this pipe with an average velocity of V = 0.0125 m/s. First, set the number of lattice as, for instance, n_x × n_y × n_z = 100 × 50 × 50, which results in a lattice length of dx = L/n_x = D/n_y = D/n_z = 0.08 mm. Note that the lattice should be square in 2D or cube in 3D, indicating that the lattice length is identical in every dimension. Then the conversion coefficient of length is decided by C_l = dx/Δx = 8 × 10⁻⁵ m. Second, given that the density of water is ρ_f = 1000 kg/m³, the conversion coefficient of density is C_ρ = ρ_f/ρ = 1000 kg/m³. Then the conversion coefficient of mass is determined by C_m = C_ρC_l³ = 5.12 × 10⁻¹⁰ kg. For better numerical stability, the relaxation time τ shall be at least 0.55, implying that the kinematic viscosity ν shall be at least 0.017. Here let us set ν = 0.05. Given that the kinematic viscosity of water is ν_water = 1.0 × 10⁻⁶ m²/s, the conversion coefficient of viscosity is C_ν = ν_water/ν = 2.0 × 10⁻⁵ m²/s. Then the conversion coefficient of time is determined as C_t = C_l²/C_ν = 3.2 × 10⁻⁴ s. So far, the conversion coefficients of the three fundamental units are decided. Other parameters that need to be nondimensionalized can be converted accordingly. As a summary, the unit conversion of this pipe flow problem is listed in

Table 1. Unit conversion of 3D pipe flow.

Physical Parameters	Unit	Real Value	Lattice Value	Conversion Coefficient
Pipe length (L)	m	0.008	100	8.0 × 10−5
Pipe diameter (D)	m	0.004	50	8.0 × 10−5
Fluid density (ρf)	kg/m3	1000	1	1000
Fluid kinematic viscosity (νf)	m2/s	1.0 × 10−6	0.05	2.0 × 10−5
Average fluid velocity (V)	m/s	0.0125	0.05	0.025
Pressure drop (ΔP)	Pa	0.2	3.2 × 10−3	62.5
Time step (Δt)	s	3.2 × 10−4	1	3.2 × 10−4

.

4. LBM-DEM Coupling

For coupling LBM with DEM, it is of significant importance to properly describe the fluid-solid interactions. The first step is to establish a lattice discretisation of the solid particle on the lattice grid. As illustrated in Figure 3, four types of lattice nodes are obtained as: (1) the pure fluid node without connection to any pure solid node; (2) the pure solid node that is not adjacent to any pure fluid node; (3) the fluid boundary node that is linked with both pure fluid node and solid boundary node; and (4) the solid boundary node that is connected with both fluid boundary node and pure solid node. Obviously, it becomes very easy and straightforward to handle any shaped particles with the discrete lattice representation. Hence, there was a variety of techniques with different levels of complexity and accuracy to compute the fluid-solid interaction, including the modified bounce-back (MBB) scheme [65,66], interpolated bounce-back (IBB) scheme [67,68,69,70,71,72] and immersed boundary methods (IBM) [73,74,75,76,77,78,79,80,81].

4.1. Modified Bounce-Back Scheme

The MBB scheme was proposed by Ladd [65,66] in order to improve the simple bounce-back rules to incorporate the influence of moving fluid-solid interfaces. As shown in Figure 3, this method assumes that the exact solid surface lies halfway between fluid boundary node and solid boundary node, where the bounce-back condition is enforced. The particle is treated as a shell full of fluid. Therefore, the bounce-back takes place on both side of the fluid-solid interface.

The first step is to estimate the velocity at the fluid-solid interface based on the solid particle’s velocity, which is given as:

$$v_b=U_p+{\mathbf{\Omega }}_p\times (x+\frac{1}{2}{e}_i\Delta t-X_p),$$

(26)

where X_p is centre of the solid particle. Then the undetermined distribution functions at the fluid boundary node and solid boundary node are computed respectively, according to the modified bounce-back rule,

$$\begin{array}{l}{f}_{-i}(x,t+\Delta t)=f_i^{+}(x,t)-2{\omega }_i\rho v_b\cdot e_i,\\ f_i(x+e_i\Delta t,t+\Delta t)=f_{-i}^{+}(x+e_i\Delta t,t)+2{\omega }_i\rho v_b\cdot e_i.\end{array}$$

(27)

It can be seen that Equation (27) reduces to the simple bounce-back scheme (Equation (25)) when v_b = 0. Once the bounce-back procedure is accomplished, the hydrodynamic force on the solid particle at this particular boundary node is calculated with the momentum exchange method,

$$F_{f,-i}(x+\frac{1}{2}{e}_i\Delta t,t+\frac{1}{2}\Delta t)=2[f_i^{+}(x,t)-f_{-i}^{+}(x+e_i\Delta t,t)-2{\omega }_i\rho v_b\cdot e_i]\cdot e_i.$$

(28)

The total hydrodynamic force and torque are computed by summing up the contributions from every lattice speed direction and boundary node,

$$\begin{array}{l}{F}_f={\displaystyle \sum _{link}{\displaystyle \sum _i{F}_{f,-i}}},\\ M_f={\displaystyle \sum _{link}{\displaystyle \sum _i(x+\frac{1}{2}{e}_i\Delta t-X_p)\times F_{f,-i}}}.\end{array}$$

(29)

Apparently, the discrete representation of the solid particle surface on the lattice grid is in a stepwise fashion. It is reported that the accuracy of MBB is of first order only due to the zig-zag staircases of a curved surface. Large fluctuations of the hydrodynamic interactions are also observed. Using a sufficiently high lattice resolution could provide much more smooth and accurate force evaluations, which, however, will lead to a huge increase of the computational cost. To overcome these drawbacks of MBB, several improved techniques were proposed, which will be introduced below.

4.2. Interpolated Bounce-Back Scheme

To precisely describe the actual boundary of the solid particle, several bounce-back schemes using the spatial interpolation were also proposed [67,68,69,70,71,72]. In all these techniques, the relative location of the exact fluid-solid interface is interpolated with a location parameter q. Hence more accurate and smoother evaluations of the hydrodynamic force and torque are obtained, which is demonstrated to maintain the second-order accuracy. Here a double interpolation scheme proposed by Yu et al. [71,72] is introduced to show how these schemes work. As shown in Figure 4, the undetermined distribution function is the one coming out of the solid boundary f_−i(x_f, t + Δt). In most of time, the solid boundary is not exactly located on a lattice node. Then the relative location of the solid boundary can be described by a weighting parameter q by means of spatial interpolation,

$$q=\left|x_f-x_b\right|/\left|x_f-x_s\right|,$$

(30)

where x_s, x_b and x_f represent the positions of the solid node, the solid boundary and the first fluid node next to the solid boundary, respectively. The weighting parameter q varies between 0 and 1. When q = 0, the solid boundary is exactly located on the nearest fluid node, which is turned into a solid node. When q = 1, the solid boundary lies just on the nearest solid node. Then the density function at the temporary location x_b that propagates exactly to the solid boundary during the streaming process can be interpolated with the existing density functions f_i(x_f, t + Δt) and f_i(x_ff, t + Δt), through a first order form,

$$f_i(x_b,t+\Delta t)=q{f}_i(x_f,t)+(1-q)f_i(x_{ff},t).$$

(31)

Next, a bounce-back operation takes place instantaneously at the solid boundary as:

$$f_{-i}(x_b,t+\Delta t)=f_i(x_b,t+\Delta t)-2{\omega }_i\rho u_b\cdot e_i,$$

(32)

which is the same as the MBB described in Equation (27). In the last step, the unknown density function f_−i(x_f, t + Δt) is determined from f_−i(x_b, t + Δt) and f_−i(x_ff, t + Δt) through a first order interpolation,

$$f_{-i}(x_f,t+\Delta t)=\frac{1}{1+q}{f}_{-i}(x_b,t+\Delta t)+\frac{q}{1+q}{f}_{-i}(x_{ff},t+\Delta t).$$

(33)

It should be noted that the interpolation can also be performed with second order form, where the density function from further fluid node x_fff is needed. The force and torque on the solid particle can be evaluated similarly with a momentum exchange method,

$$\begin{array}{l}{F}_f={\displaystyle \sum _{all x_f}{\displaystyle \sum _i[f_i^{+}(x_f,t)+f_{-i}(x_f,t+\Delta t)]e_i}},\\ M_f={\displaystyle \sum _{all x_f}{\displaystyle \sum _i(x_b-X_p)\times [f_i^{+}(x_f,t)+f_{-i}(x_f,t+\Delta t)]e_i}}.\end{array}$$

(34)

4.3. Immersed Moving Boundary Method

The immersed boundary method (IBM) was proposed by Peskin [73,74] and further implemented in LBM by other researchers [75,76,77,78,79,80]. In this method, extra moving Lagrangian nodes are introduced to represent the solid particle shape, which are assumed to be deformable with a large stiffness. The effects of the immersed particle boundary on the fluid are modelled by restoring forces based on the ‘no-slip’ boundary condition, which tend to keep the particle to its original shape. Then the restoring forces are distributed to their surrounding fluid nodes and considered as the external force terms in the governing equations. The drawbacks of the IBM include the introduction of additional free parameters, the problem-sensitive determination of the spring stiffness and the damping constant, as well as the severely restricted computational time step.

As an alternative, Noble and Torczynski [81] proposed an immersed moving boundary technique based on the local solid fraction in each lattice cell, which is also known as the partially-solid scheme. This method not only overcomes the momentum discontinuity of MBB-based techniques and provides an adequate representation of complex boundaries at relative lower lattice resolutions, but also retains the prominent advantages of the LBM, i.e., the simple linear collision operator and its locality in computation. In this method, the lattice Boltzmann equation is modified to include an additional collision term, which depends on the local volume fraction of solid (see Figure 3 right part),

$$f_i(x+e_i\Delta t,t+\Delta t)=f_i(x,t)-(1-B_n)[\frac{\Delta t}{\tau }(f_i(x,t)-f_i^{eq}(x,t))]+{\displaystyle \sum _s{B}_s{\mathrm{\Omega }}_i^s}+(1-B_n)F_i\Delta t,$$

(35)

where B_n is a total weighting function in each lattice cell, B_s is the weighting function from each solid particle in the same lattice cell, and Ω_i^s is the additional collision term. The total weighting function is calculated by summing up the weighting function of each solid particle that intersect with the same lattice cell such that $B_n={\displaystyle {\sum }_s{B}_s}$. Each particle’s weighting function is determined by its own solid fraction and the relaxation parameter as:

$$B_s({\epsilon }_s,\tau )=\frac{{\epsilon }_s(\tau /\Delta t-0.5)}{(1-{\epsilon }_n)+(\tau /\Delta t-0.5)},$$

(36)

where ε_s is the corresponding particle’s solid fraction and ε_n is total solid fraction in the same lattice cell that is the sum of each particle’s contribution, i.e., ${\epsilon }_n={\displaystyle {\sum }_s{\epsilon }_s}$. It can be seen that when the solid fraction ε_s varies from 0 (a completely fluid cell) to 1 (a completely solid cell), B_s varies from 0 to 1. Equation (35) returns to the original SRT-LBE for pure fluid when B_s = 0, and returns the new collision operator ${\mathrm{\Omega }}_i^s$ plus the distribution from the previous time step when B_s = 1. The new collision operator is given by

$${\mathrm{\Omega }}_i^s=f_{-i}(x,t)-f_{-i}^{eq}(\rho ,u)+f_i^{eq}(\rho ,U_p)-f_i(x,t).$$

(37)

The total hydrodynamic force and torque acting on the solid particle can be evaluated using the momentum exchange method with the additional collision operator over all lattice directions at each node and then over all fluid boundary, solid boundary and internal solid nodes, which are expressed as:

$$\begin{array}{l}{F}_f=-{\displaystyle \sum _n{B}_s({\displaystyle \sum _i{\mathrm{\Omega }}_i^s{e}_i)}},\\ T_f=-{\displaystyle \sum _n[(x_n-X_p)\times B_s({\displaystyle \sum _i{\mathrm{\Omega }}_i^s{e}_i)}}].\end{array}$$

(38)

Note that the implementation of the immersed moving boundary method is straightforward. Only quantities already available in the lattice cell or easily derived are used. No additional data storage or organization is needed, which is a crucial advantage over other moving boundary techniques. However, it is also noticed that the critical step in this method is the calculation of the solid fraction in each lattice cell, which is primarily a geometry problem. Many methods are proposed to estimate the solid coverage ratio of a geometry with a square or a cube, including analytically solution, Monte Carlo sampling technique, cell decomposition method, polygonal approximation, edge-intersection averaging method, linear approximation and so on [82,83]. In the current LBM-DEM numerical framework, the analytical solution is used for 2D problems, while for 3D a linear approximation method is adopted, in order to save computational cost but still retain relatively accurate results.

4.4. Time Steps in the LBM-DEM Coupling

The time step coupling is also a crucial issue in the LBM-DEM coupling [31,82]. As we discussed previously, the time step of LBM is generally set as a dimensionless value of one in the computation. The real physical time is implicitly determined by computational parameters as well as the unit conversion scheme. Based on Li and Marshall’s work [3,58], the DEM time step depends on the interparticle collision time scale,

$$t_c\approx R{(\frac{{\rho }_p^2}{E^2{v}_R})}^{1/5},$$

(39)

which typically varies around 10⁻⁶ ~ 10⁻⁹ s. Therefore, it is believed that the LBM time step is generally larger than the DEM time step. A time step ratio λ is then introduced as Δt_LBM = λΔt_DEM, which can be decided after all the parameters are settled. It should be noted that the time step ratio λ could be either greater or lower than one. If λ < 1, the critical time step can be simply set as Δt_LBM, and the DEM time step can be equal to Δt_LBM. When λ > 1, a number of DEM computation steps are performed within one LBM time step, during which the fluid forces and torques are kept unchanged.

5. Validation

In this section, a few validation studies are performed to demonstrate the accuracy of our LBM-DEM numerical approach, including single-phase Poiseuille flow, gravitational settling of a particle, and the drag force on a stationary particle.

5.1. Poiseuille Flow

Single-phase Poiseuille flows in both 2D and 3D are first analysed using the LBM-DEM. The validation is performed with both SRT and MRT models. As shown in Figure 5, the boundaries in are periodic, and y-direction are bounded by two ‘no-slip’ walls. A constant pressure gradient is imposed in the channel in x-direction to drive the flow. The channel size is H × H for the 2D case and 10 × H × H for the 3D case. A few lattice resolutions are selected in the range H = 11 ~ 101 to examine the convergence of the numerical model. The relaxation parameter is fixed at 0.65. Different channel Reynolds numbers are achieved by changing the magnitude of the body force. Figure 6 and Figure 7 show the normalised velocity profiles for 2D and 3D Poiseuille flow, respectively. The L-2 norm calculation is used to compute the relative error:

$$relative error={\displaystyle \sum _n\sqrt{\frac{{(U-U_{theory})}^2}{U_{theory}^2}}},$$

(40)

where U is the flow velocity obtained from the simulation and U_theory is the analytical solution of the Poiseuille flow [84]. It can be inferred that the numerical results agree very well with the theoretical prediction. The accuracy is second order for the 2D flow, while it is slightly lower than second order for the 3D duct flow.

5.2. Gravitational Settling of a Particle

For gravitational settling of a particle, both 2D and 3D cases are simulated to evaluate the validity of our numerical approach in a dynamic system. This validation is performed with SRT-LBM along with the Hertz model in DEM. For the 2D case, the computational setup is identical to that in Wen et al.’s work [85]. As shown in Figure 8, a cylinder with diameter d_p = 0.1 cm is initially located 0.076 cm away from the left wall of a vertical channel with width 0.4 cm. The fluid density and kinematic viscosity are 1000 kg/m³ and 1 × 10⁻⁶m²/s, respectively. The mass density of the cylinder is 1030 kg/m³. Therefore, the cylinder will settle under the gravitational force and finally reaches a steady state that moves along the centreline at a constant velocity. The computational setup in the lattice unit scheme is as follows. The size of the channel is n_x × n_y = 121 × 1201 and ‘no-slip’ wall boundaries are imposed on all the four faces. The cylinder’s diameter and mass density are d_p = 30 and ρ_p = 1.03, respectively. The dimensionless relaxation parameter is set as 0.6. Note that the IBM is applied in the solid–fluid coupling in this particular problem.

A detailed comparison between the present numerical results and that obtained by the arbitrary Lagrangian–Eulerian technique (ALE) [85] is presented in Figure 9, where very good agreement is reached. Furthermore, the force profiles are quite smooth without any large fluctuations.

For the 3D case, the single particle settling experiment reported by Ten Cate et al. [86] is numerically reproduced. A spherical particle with diameter d_p = 15 mm and density ρ_p = 1,120 kg/m³ is initially released from a height h = 120 mm in a cuboid box with a size of 100 × 100 × 160 mm³. The box is full of silicon oil, and four different types of oil are used in the experiment. The densities and dynamic viscosities of the silicon oil are ρ_f = 970, 965, 962, 960 kg/m³ and ν_f = 373, 212, 113, 58 Pa·s, respectively. The terminal settling velocity of the particle in the four oils are u_ter = 0.038, 0.060, 0.091, 0.128 m/s, respectively, which results in four different particle Reynolds number Re_p = 1.5, 4.1, 11.6, 31.9, respectively. The computational parameters in the dimensionless lattice unit are given below. The domain size is 51 × 51 × 81 and the particle diameter is d_p = 7.5. The density of the fluid is fixed at ρ_f = 1.0, so that the corresponding particle density is ρ_p = 1.155, 1.161, 1.164, and 1.167 for different fluids considered, respectively. The relaxation parameter is set as 0.65. Figure 10 shows the time evolution of particle trajectory and velocity profiles for different Reynolds numbers. It is apparent that the numerical results are in excellent agreement with the experimental results, which further demonstrates the validity and accuracy of the present numerical approach.

5.3. Drag Force on a Stationary Particle

The classical drag force problem is also simulated to verify the accuracy of the force evaluation in our LBM-DEM coupling, which is carried out with the SRT-LBM. As depicted in Figure 11, a circular (spherical) particle is placed in the centre of a rectangular (cuboid) domain and kept stationary. The domain size is L × H = 50d_p × 50d_p in 2D and L × H × H = 20d_p × 10d_p × 10d_p in 3D. The inlet boundary is constant flow with velocity U₀ and the outlet is set as constant pressure boundary. All the other boundaries are set as open boundary (zero gradient). The particle Reynolds number and the drag coefficient are calculated as

$$\begin{array}{l}{\mathrm{Re}}_p=\frac{U_0{d}_p}{\nu },\\ C_d=\frac{8F_d}{{\rho }_f{U}_0^2\pi d_p^2},\end{array}$$

(41)

where F_d is the fluid drag force on the particle. Note that both the IBB and IBM are used in the drag force validation.

Figure 12 shows the drag coefficient C_d as a function of the particle Reynolds number Re_p for both 2D and 3D cases. For the 2D case shown in Figure 12a, the experimental results reported by Tritton are also superimposed for comparison [87]. It is clear that both the results obtained by IBB and IBM agree well with the experiments. Second order accuracy of the force evaluation is achieved for IBM, while the accuracy for IBB is lower than second order but higher than first order. For the 3D case shown in Figure 12a,b widely accepted empirical law of the drag coefficient is introduced to serve as the theoretical prediction [88]

$$C_d=\frac{24}{{\mathrm{Re}}_p}(1+0.15{\mathrm{Re}}_p^{0.687}).$$

(42)

We can see that the simulation results agree well with the above theoretical prediction, which demonstrates the validity and accuracy of our coupled LBM-DEM. Furthermore, it is noticed that accurate drag coefficients are obtained using both IBB and IBM even with a relatively low size resolution d_p = 6, implying that acceptable accuracy can still be retained when the computational cost is reduced.

6. LBM-DEM Applications

In this section, a few numerical examples are presented to demonstrate the capability of LBM-DEM, which include inertial migration of dense particle suspensions, agglomeration of adhesive particles in channel flow, and sedimentation of particle suspensions in a cavity flow.

6.1. Inertial Migration of Dense Particle Suspensions

In an experimental study of pipe flow of a dilute suspension consisting of neutrally buoyant particles, Segré and Silberberg [89,90] first observed that the particle suspensions tend to migrate laterally and focus in an annulus at radial position r ≈ 0.6R, with R being the pipe radius. This phenomenon is named as the tubular pinch effect (or the Segré–Silberberg effect). Due to the lack of theoretical explanation at the time when it was discovered, this interesting phenomenon prompted a great deal of interest to further explore the underlying mechanism. Since then, many studies were carried out on this fascinating phenomenon experimentally [91,92,93,94,95], theoretically [96,97,98,99] and computationally [100,101,102,103,104,105,106]. The lateral focusing of the particle suspensions is found to be closely related to the fluid inertia. It is well recognised that the lateral migration exists in both 2D and 3D flows, and the equilibrium position moves closer to the wall as the channel Reynolds number increases, which is successfully described with the theoretical solution based on the perturbation theory and the asymptotic expansion method. However, most of the previous works are limited to very dilute suspensions with a particle concentration lower than 1%. When the particle concentration is increased, whether the lateral focusing phenomenon still exists remains less explored. Therefore, the coupled LBM-DEM is employed in the current study to explore this problem.

As shown in Figure 13, let us consider a pressure-driven flow of non-Brownian particle suspensions in a 2D channel. Periodic boundary conditions are applied in x-direction and the top and bottom planes in y-direction are set as ‘no-slip’ walls. The particles are neutrally buoyant, i.e., the fluid and the particle have the same mass density. The particles’ initial positions are randomly generated inside the channel at the beginning of the simulation. Driven by the pressure gradient, the particle suspensions will transport and migrate to their equilibrium positions. A steady particle–fluid flow is expected after a sufficient long time of computation. The parameters used in the simulation are as follows. The size of the channel is L × H = 501 × 101. The particle diameter is d_p = 6 and 12. The particle concentration ϕ ranges between 1% and 50%. The channel Reynolds number Re₀ is tuned in the range 4 ~ 100 by varying the pressure gradient. Note that the channel Reynolds number refers to the one under a single phase flow with no particles. The SRT model and the Hertz contact model are used in the LBM and DEM for this particular problem, respectively.

Figure 14 displays the snapshots of the particle suspension with different channel Reynolds numbers, where the particle concentration is fixed at ϕ = 10%. It is observed that only a small part of the particles undergoes the lateral migration at a relatively small Reynolds number Re₀ = 4. As Re₀ increases, a complete migration takes place at both Re₀ = 40 and Re₀ = 100, where all the particles are focused at a certain lateral position. Furthermore, it is found out that the larger Re₀ is, the shorter time it takes to develop into a full migration. For example, it occurs around t = 5.5 × 10⁴ for Re₀ = 40, while it is much earlier around t = 4 × 10⁴ for Re₀ = 100. Figure 15 illustrates the effects of particle concentration on the particle migration. It is clear that with relatively low particle concentration (ϕ = 1% and 10%), the lateral migration is still notable. However, as ϕ increases to 40%, the migration is dramatically suppressed and the particle suspensions appear to jam near the wall.

To quantify the influence of the particle concentration, the degree of inertial migration can be defined as [94],

$$P_f={\displaystyle \sum PDF(\left|y\right|\ge \frac{1}{2}H)},$$

(43)

which estimates the total probability distribution function (PDF) of the particles in the upper and lower quarter of the channel. Since the lateral equilibrium position is around 0.6 according to the Segré and Silberberg effect, the degree of the migration must be P_f = 1 if all the particles are laterally focused. On the other hand, if the particle suspension remains randomly distributed in the channel, P_f is expected to be 0.5. Figure 16 shows the variation of the degree of inertial migration with particle concentration for different d_p and Re₀. It is observed that P_f equals one for all the Re₀ at a low particle concentration (ϕ ≤ 5%), while P_f remains around 0.5 ~ 0.6 at ϕ = 50% for all the Re₀, implying a huge impact of the particle concentration on the lateral migration. A full migration is always accessible when the particle suspension is not dense, but it is completely suppressed with a very large particle concentration for the range of Re₀ in the present study. However, when ϕ is between 5% and 50%, the variation of P_f strongly depends on Re₀. With a fixed ϕ, a larger Re₀ leads to a larger P_f, i.e., a higher degree of lateral migration.

The above discussion reveals that the degree of the migration is determined by the particle concentration and channel Reynolds number. A dimensionless focusing number was proposed to characterise the degree of migration [37],

$$F_c=\frac{{\mathrm{Re}}_0^m}{{\varphi }^n},$$

(44)

where m = 0.36 and n = 2.33 are obtained from the fitting of simulation data. Figure 17 shows the degree of migration P_f as a function of the focusing number F_c, where an empirical fitting is given by an exponential equation $P_f=1-0.5\times {10}^{-0.0043\times {(\log F_c)}^{6.51}}$. It can be seen that the numerical results coalesce onto a single master curve and three distinct regimes can be identified based on the value of F_c: a laterally unfocused regime (F_c < 30), a partially migration regime (30 < F_c < 300) and a fully migration regime (F_c > 300).

6.2. Agglomeration of Adhesive Particles in Channel Flow

The inertial focusing introduced in the previous section has prompted a variety of novel applications in the microfabrication, such as particle filtration, segregation and flow cytometry in the micro channels [107,108,109]. An inevitable problem of particle agglomeration due to the cohesive nature when the size reduces to the micron scale must be taken into consideration. The agglomeration is ubiquitous in nature and industry, and sometimes will cause unfavourable effects in industrial facility. For example, the agglomeration of asphaltenes or gas hydrates results in the pipeline blockage in the oil and gas engineering [110,111]. Therefore, it is worth investigating the fundamental mechanism of micro-particle agglomeration.

As displayed in Figure 18, a stream of dilute particle suspension flowing through a 3D square channel is considered. Periodic boundary conditions are imposed at the inlet and outlet, while other four faces of the channel are set as ‘no-slip’ walls. The flow is driven by a constant pressure gradient. The particles are neutrally buoyant and adhesive, with initial positions randomly distributed in the channel. The computational parameters are set as follow. The dimension of the channel is L × H × H = 201 × 61 × 61. By varying the pressure gradient, the channel Reynolds number ranges between 5.4 and 108. The particle size is d_p = 5 and the particle concentration is fixed at 1%, corresponding to a particle number of 110. The strength of the particle adhesion is represented by the surface energy, which is chosen as 0.075 and 0.0075. The SRT model and the JKR contact theory are used in the LBM and DEM for this problem, respectively.

Figure 19 and Figure 20 display the snapshots of adhesive particle suspensions with different surface energies γ = 0.075 and γ = 0.0075, respectively. It can be found that large sized agglomerates are formed with relatively higher surface energy (γ = 0.075). Some particles even stick on the wall due to strong adhesion. However, when the surface energy is reduced by one order of magnitude (γ = 0.0075), the size of the agglomerates becomes smaller and the particles appear to focus around a certain lateral position with Re = 54 and 108 (see Figure 20c,d), which resembles the Segré-Silberberg effect as discussed in the last section.

Figure 21 shows the lateral PDF of the particle suspensions based on the following definition,

$$f(l)=\frac{N(l,l+\Delta l)}{N},$$

(45)

where N is the total number of particles in the channel and N(l, l+Δl) is the number of particles in the square annulus between l and l + Δl. For comparison, the channel flow with non-adhesive particles is also modelled and the results are included. It can be seen from Figure 21a that a single peak at the lateral position around 0.6 is observed in the PDF for non-adhesive particles, which moves closer to the wall as Re increases and is in consistent with the Segré and Silberberg effect. However, the PDFs are quite distinct for adhesive particles. The single peak distribution is not clearly observed for particles with γ = 0.075 (see Figure 21b). Moreover, there appears to be a peak close to the wall position, where a number of particles are found to stick to wall because of the strong adhesion, as indicated by Figure 19. On the other hand, for much less adhesive particles (γ = 0.0075), the PDFs become similar to those of non-adhesive particles (see Figure 21c). The average lateral position as a function of Re is further shown in Figure 22. Note that the relatively larger value of average lateral position for Re = 5.4 is caused by the statistical effect [38], because the fluid inertia is relatively small to drive the lateral migration and the particles almost remain randomly distributed in the cross-sectional area. Except for the data for Re = 5.4, it is observed that the average lateral position increases monotonically with the increase of Re for both non-adhesive and less adhesive particles (γ = 0.0075), which demonstrates a similar dynamics between non-adhesive and weak adhesive particles. However, for strongly adhesive particles (γ = 0.075), the monotonic increase of the average lateral position disappears because of the agglomeration.

The agglomeration mechanism for such channel flow can be attributed to the competition between the interparticle adhesive contact force and the hydrodynamic force. The adhesive force sticks particles together, which imposes a positive effect on the agglomeration. The hydrodynamic force exerts two opposite effects on the particle suspensions. First, it induces the particle lateral migration, which motivates the interparticle collision and facilitates the agglomeration. Secondly, the fluid inertia also tears particles apart from each other, resulting in the breakup of agglomerates. As a consequence, a new dimensionless adhesion number Ad is proposed based on the balance of the critical pull-off force and the drag force [11,12,38],

$$Ad=\frac{\gamma }{{\mu }_{f}U},$$

(46)

where μ_f is the fluid dynamic viscosity and U is the mean fluid velocity. Using the adhesion number, the degree of the agglomeration can be characterised and represented by the agglomerate ratio that is defined as the ratio of the number of particles involved in agglomerate to the total particle number. Figure 23 shows the relationship between the agglomerate ratio and the adhesion number. A monotonic increase of the agglomerate ratio as a function of Ad is observed. After Ad reaches a critical value, the agglomerate ratio seems to converge to one. Therefore, two different regimes can be identified, i.e., a partial agglomeration regime where the hydrodynamic force dominates over the adhesive force, and a complete agglomeration regime where all the particles are involved in agglomerates.

6.3. Sedimentation of Particle Suspensions in Cavity Flow

Sedimentation is extensively used in industries to separate particles from fluid or other particles with different sizes or velocities, such as separating particles and cells in microfluids [112], dewatering coal slurries [113], and post-treating wastewater [114,115,116]. Furthermore, sediments are also physical pollutants in the river, which could have great influence on the ecosystem and human health [117,118,119,120]. Therefore, it is of great importance to study the sedimentation, which is a very complicated process, due to the complex fluid mechanics involved with the sediment particles. With LBM-DEM, a numerical investigation is then performed to gain some insights to the sedimentation problem.

As shown in Figure 24, we consider a sedimentation system composed of a cuboid channel and a cavity, which is placed in the middle of the channel. The boundaries in x and y directions are periodic, and all the other faces are set as no-slip walls. The particles are initially generated at random positions in the main channel (no particles in the cavity). The particle density is larger than the fluid, so that they can sedimentate into the cavity under the gravity. The size of the main channel is L × W × H = 321 × 81 × 49. The cavity has the same width as the main channel. The length of the cavity is l = 40 and 80, while the height is h = 16 and 32, which gives rise to four configurations of the cavity. The channel Reynolds number, defined as Re = U_xH/ν varies between 28.44 and 88.32. The particle diameter is d_p = 6 and the density is in the range 1.1 ~ 2.0. The total number of particles is fixed at 110, which is equal to an overall concentration less than 1%. The SRT model and the Hertz model are used in the LBM and DEM for this particular problem, respectively.

Figure 25 shows the snapshots of the particle suspensions at different time. The particles are initially located at random positions in the main channel with no particles trapped inside the cavity. Then the particles settle to the bottom of the channel and flow over the cavity. During the transportation of the suspensions, some particles fall into the cavity and stay stationary inside, while some other keep circulating around an orbit in the central vortex inside the cavity. The total number of particles trapped inside the cavity is evaluated after a sufficiently long time, when the transportation of the particles reaches a steady state. From Figure 25 the typical trajectories of particle suspensions in the cavity are analysed, where three distinct dynamic behaviours can be identified: resuspension, deposition and circulation. Resuspension means that the particles fall into the cavity and then come out due to the lift force from the fluid. Deposition happens mostly at the rear edge of the cavity, where the particles deposit to the corner of the cavity and stay motionless. Circulation indicates that the particles are captured by the central vortex inside the cavity and keeps moving along an orbit in a periodic way. Both deposition and circulation are treated as the entrapment of a particle.

Figure 26 shows the trap efficiency η, defined as the ratio of the trapped particle number to the total particle number, as a function of the particle density for different cavity configurations and Reynolds numbers. It is clear that as the particle density increases, the trap efficiency increases continuously from 0 to 1. For a fixed Reynolds number, the largest cavity (h = 32, l = 80) has the highest trap efficiency, while the smallest cavity (h = 16, l = 40) has the lowest trap efficiency. The results of the other two cavity configurations are in between, which indicates that increasing either length or depth leads to a higher trap efficiency. Furthermore, for a fixed cavity configuration, an increase in the Reynolds number results in a decrease in the trap efficiency. With a larger Reynolds number, the particle’s translational velocity in x direction becomes larger, so that it has less time to enter the cavity in the vertical direction. Meanwhile, a larger Reynolds number produces a larger lift force. As a consequence, the particles are less easily trapped in the cavity.

7. Summary

An overview of the coupled discrete element method with the lattice Boltzmann method is presented together with its applications in solid–liquid flows. The fundamentals of DEM and LBM, including the contact mechanics based on Hertz model and JKR theory, the lattice Boltzmann equation with both single-relaxation-time model and multi-relaxation-time model are introduced. Several solid–fluid interaction models, i.e., the coupling between DEM and LBM, are discussed and compared in detail. The validity and accuracy of the numerical method are validated for three classical solid–liquid flow problems, and the results are in good quantitative agreement with those from experiments and other numerical approaches. Three case studies, including the inertial migration of particle suspensions, the agglomeration of adhesive particles, and the sedimentation problem, are also performed to demonstrate the capability of the coupled DEM-LBM to solid–liquid flow problems and potential engineering applications.