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

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.


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

Discrete Element Method (DEM)
In DEM, particle motion is solved individually using the Newton's equation of motion, 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, ChemEngineering 2020, 4, 55 3 of 33 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 F collision = F n n + F s t s , M collision = r p F s (n × t s ) + M r (t r × n) + M t n, (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.

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: where k N is normal stiffness coefficient in the Hertz model,â = 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 √ R, where R is the effective particle radius and E is the effective elastic moduli, defined as: 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 = √ Rδ 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 0 , 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., where η N is the normal dissipation coefficient, and v R is the relative velocity at the contact point.

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, where k T is the spring stiffness coefficient in the tangential direction, v s is relative sliding velocity, ξ T = t t 0 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: where ξ R = t t 0 v L (t) · t r dt and ξ Q = t t 0 Ω 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: M r,crit = k R θ crit R, M t,crit = 3πaF s,crit /16, (9) where θ crit is the critical rolling angle in the relative rolling motion between two contacting particles.

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 fi(x,t) follows the lattice Boltzmann equation (LBE) [19] with an external force term, where x is the position of the node, ei is the unit vector of the lattice speed in the direction i as denoted in Figure 1, Δt is the explicit time step, Fi is the external body force and Ωi[fi(x,t)] is the collision operator. Equation (10) can be evaluated in two independent processes: collision and streaming, as shown below: .
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[fi(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 + eiΔt according to the lattice speed model. The sequence of collision and streaming becomes irrelevant after a few time steps.

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]: (13) where τ is the dimensionless relaxation parameter and fi eq (x,t) is the equilibrium distribution function that is defined as: In this equation, the weight coefficient ωi depends on the lattice speed model. For example, ω0 = 4/9, ω1,2,3,4 = 1/9 and ω5, 6  The discrete body force term Fi in Equation (10) is given as [63]: The evolution of the particle distribution functions f i (x,t) follows the lattice Boltzmann equation (LBE) [19] with an external force term, 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: 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.
ChemEngineering 2020, 4, 55 6 of 33 The discrete body force term F i in Equation (10) is given as [63]: 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: The kinematic viscosity of the fluid ν is related to the relaxation parameter as:

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]: where S is the force term. From Equation (15), we have: Equation (18) can be transformed into the matrix form as: 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: ChemEngineering 2020, 4, 55 The equilibrium distribution functions m eq in the moment space are related to that in the velocity space as follows: ρ e eq ε eq j eq x q eq x j eq y q eq y p eq xx p eq xy The diagonal matrix Λ is Λ = diag(s 1 , s 2 , s 3 , s 4 , s 5 , s 6 , s 7 , s 8 , s 9 ). In order to ensure the mass and momentum conservation after collision, it is easily derived that s 1 = s 4 = s 6 = 0. s 8 and s 9 are set as s 8 = s 9 = 1/τ in order to reproduce the same viscosity of SRT model [64]. The other relaxation parameters s 2 , s 3 , s 5 and s 7 can be decided with more flexibility according to the physical model.

of 33
The equilibrium distribution functions m eq in the moment space is: The diagonal matrix Λ is Λ = diag(s 1 , s 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., 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 f1, f5 and f8 at the outlet boundary will stream into the corresponding node at the inlet, while f3, f6 and f7 at the inlet propagate into the corresponding node at the outlet.

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, nx × ny × nz = 100 × 50 × 50, which results in a lattice length of dx = L/nx = D/ny = D/nz = 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 Cl = dx/Δx = 8 × 10 −5 m. Second, given that the density of water is ρf = 1000 kg/m 3 , the conversion coefficient of density is Cρ = ρf/ρ = 1000 kg/m 3 . Then the conversion coefficient of mass is determined by Cm = CρCl 3 = 5.12 × 10 −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 −6 m 2 /s, the conversion coefficient of viscosity is Cν = νwater/ν = 2.0 × 10 −5 m 2 /s. Then the conversion coefficient of time is determined as Ct = Cl 2 /Cν = 3.2 × 10 −4 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. 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 1 , f 5 and f 8 at the outlet boundary will stream into the corresponding node at the inlet, while f 3 , f 6 and f 7 at the inlet propagate into the corresponding node at the outlet.

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 −5 m. Second, given that the density of water is ρ f = 1000 kg/m 3 , the conversion coefficient of density is C ρ = ρ f/ ρ = 1000 kg/m 3 . Then the conversion coefficient of mass is determined by C m = C ρ C l 3 = 5.12 × 10 −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 −6 m 2 /s, the conversion coefficient of viscosity is C ν = ν water /ν = 2.0 × 10 −5 m 2 /s. Then the conversion coefficient of time is determined as C t = C l 2 /C ν = 3.2 × 10 −4 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.

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].

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].

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

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: 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, 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, The total hydrodynamic force and torque are computed by summing up the contributions from every lattice speed direction and boundary node, 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.

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, 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, Next, a bounce-back operation takes place instantaneously at the solid boundary as: which is the same as the MBB described in Equation (27). In the last step, the unknown density function 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, ChemEngineering 2020, 4, 55 12 of 33 Next, a bounce-back operation takes place instantaneously at the solid boundary as: x u e (32) which is the same as the MBB described in Equation (27). In the last step, the unknown density It should be noted that the interpolation can also be performed with second order form, where the density function from further fluid node xfff is needed. The force and torque on the solid particle can be evaluated similarly with a momentum exchange method,

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.

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), 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 = s B s . Each particle's weighting function is determined by its own solid fraction and the relaxation parameter as: 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., ε n = s ε 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 Ω s i plus the distribution from the previous time step when B s = 1. The new collision operator is given by 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: 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.

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, which typically varies around 10 −6~1 0 −9 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.

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.

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. Figures 6 and 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: 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.

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 3 and 1 × 10 −6 m 2 /s, respectively. The mass density of the cylinder is 1030 kg/m 3 . 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.

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 dp = 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 3 and 1 × 10 −6 m 2 /s, respectively. The mass density of the cylinder is 1030 kg/m 3 . 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 nx × ny = 121 × 1201 and 'no-slip' wall boundaries are imposed on all the four faces. The cylinder's diameter and mass density are dp = 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 dp = 15 mm and density ρp = 1,120 kg/m 3 is initially released from a height h = 120 mm in a cuboid box with a size of 100 × 100 × 160 mm 3 . 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 3 and νf = 373, 212, 113, 58 Pa·s, respectively. The terminal settling velocity of the particle in the four oils are uter = 0.038, 0.060, 0.091, 0.128 m/s, respectively, which results in four different particle Reynolds number Rep = 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 dp = 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. 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 3 is initially released from a height h = 120 mm in a cuboid box with a size of 100 × 100 × 160 mm 3 . 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 3 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.

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 where Fd is the fluid drag force on the particle. Note that both the IBB and IBM are used in the drag force validation. Figure 11. Schematic of flow past a stationary particle. Figure 12 shows the drag coefficient Cd as a function of the particle Reynolds number Rep 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

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 0 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 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.

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 Re , where Fd is the fluid drag force on the particle. Note that both the IBB and IBM are used in the drag force validation. Figure 11. Schematic of flow past a stationary particle. Figure 12 shows the drag coefficient Cd as a function of the particle Reynolds number Rep 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 Figure 11. Schematic of flow past a stationary particle. 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 = 24 Re p (1 + 0.15Re 0.687 p ). (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.
ChemEngineering 2020, 4, 55 19 of 33 case shown in Figure 12a,b widely accepted empirical law of the drag coefficient is introduced to serve as the theoretical prediction [88] = +  (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 dp = 6, implying that acceptable accuracy can still be retained when the computational cost is reduced. Figure 12. Drag coefficient as a function of particle Reynolds number for (a) 2D and (b) 3D cases.

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.

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

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.

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 0 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.
ChemEngineering 2020, 4, 55 20 of 33 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 dp = 6 and 12. The particle concentration ϕ ranges between 1% and 50%. The channel Reynolds number Re0 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 13. Schematic of pressure-driven flow of dense particle suspensions. 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 Re0 = 4. As Re0 increases, a complete migration takes place at both Re0 = 40 and Re0 = 100, where all the particles are focused at a certain lateral position. Furthermore, it is found out that the larger Re0 is, the shorter time it takes to develop into a full migration. For example, it occurs around t = 5.5 × 10 4 for Re0 = 40, while it is much earlier around t = 4 × 10 4 for Re0 = 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.   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 0 = 4. As Re 0 increases, a complete migration takes place at both Re 0 = 40 and Re 0 = 100, where all the particles are focused at a certain lateral position. Furthermore, it is found out that the larger Re 0 is, the shorter time it takes to develop into a full migration. For example, it occurs around t = 5.5 × 10 4 for Re 0 = 40, while it is much earlier around t = 4 × 10 4 for Re 0 = 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.
ChemEngineering 2020, 4, 55 20 of 33 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 dp = 6 and 12. The particle concentration ϕ ranges between 1% and 50%. The channel Reynolds number Re0 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 13. Schematic of pressure-driven flow of dense particle suspensions. 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 Re0 = 4. As Re0 increases, a complete migration takes place at both Re0 = 40 and Re0 = 100, where all the particles are focused at a certain lateral position. Furthermore, it is found out that the larger Re0 is, the shorter time it takes to develop into a full migration. For example, it occurs around t = 5.5 × 10 4 for Re0 = 40, while it is much earlier around t = 4 × 10 4 for Re0 = 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], 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 0 . It is observed that P f equals one for all the Re 0 at a low particle concentration (φ ≤ 5%), while P f remains around 0.5~0.6 at φ = 50% for all the Re 0 , 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 0 in the present study. However, when φ is between 5% and 50%, the variation of P f strongly depends on Re 0 . With a fixed φ, a larger Re 0 leads to a larger P f , i.e., a higher degree of lateral migration. 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. Figure 14. Snapshots of the migration process of particle suspensions with different channel Reynolds numbers, (a) Re0 = 4, (b) Re0 = 40, and (c) Re0 = 100. The particle diameter is dp = 6 and the concentration is ϕ = 10%. Figure 15. Snapshots of the particle positions with different particle concentrations, (a) φ = 1%, (b) φ = 10%, (c) φ = 40%. The particle diameter is d p = 6 and all the snapshots are captured at the steady state.
To quantify the influence of the particle concentration, the degree of inertial migration can be defined as [94], ), 2 f P PDF y 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 Pf = 1 if all the particles are laterally focused. On the other hand, if the particle suspension remains randomly distributed in the channel, Pf is expected to be 0.5. Figure 16 shows the variation of the degree of inertial migration with particle concentration for different dp and Re0. It is observed that Pf equals one for all the Re0 at a low particle concentration (ϕ ≤ 5%), while Pf remains around 0.5 ~ 0.6 at ϕ = 50% for all the Re0, 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 Re0 in the present study. However, when ϕ is between 5% and 50%, the variation of Pf strongly depends on Re0. With a fixed ϕ, a larger Re0 leads to a larger Pf, i.e., a higher degree of lateral migration. Figure 16. The degree of inertial migration Pf as a function of particle concentration for different particle sizes dp and channel Reynolds numbers Re0. The data was originally reported in [37].
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], where m = 0.36 and n = 2.33 are obtained from the fitting of simulation data. Figure 17 shows coalesce onto a single master curve and three distinct regimes can be identified based on the value of Fc: a laterally unfocused regime (Fc < 30), a partially migration regime (30 < Fc < 300) and a fully migration regime (Fc > 300). Figure 16. The degree of inertial migration P f as a function of particle concentration for different particle sizes d p and channel Reynolds numbers Re 0 . The data was originally reported in [37].
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], 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 × 10 −0.0043×(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).
ChemEngineering 2020, 4, 55 22 of 33 Figure 17. The degree of inertial migration as a function of the focusing number. The data was originally reported in [37].

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 dp = 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.

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.
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 dp = 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.   Figures 19 and 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. Figures 19 and 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, (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  Figure 21 shows the lateral PDF of the particle suspensions based on the following definition, 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.
ChemEngineering 2020, 4, 55 25 of 33 to drive the lateral migration and the particles almost remain randomly distributed in the crosssectional 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 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 to drive the lateral migration and the particles almost remain randomly distributed in the crosssectional 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 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 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], 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.
ChemEngineering 2020, 4, 55 26 of 33 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.

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 = UxH/ν varies between 28.44 and 88.32. The particle diameter is dp = 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.

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 x H/ν 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 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 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.
ChemEngineering 2020, 4, 55 28 of 33 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.

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.

Acknowledgments:
The authors acknowledge Duo Zhang and Nicolin Govender in the University of Surrey for their helpful suggestions and fruitful discussion in developing the numerical approach.

Conflicts of Interest:
The authors declare no conflict of interest.

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.