Particle Flow Simulation Based on Hybrid IMB-DEM-LBM Approach with New Solid Fraction Calculation Scheme

A new coupling method, immersed moving boundary–discrete element method–lattice Boltzmann method (IMB-DEM-LBM), is proposed to simulate particle flow for application in soil mechanics or coastal engineering. In this study, LBM fluid is simulated on the regular Eulerian grid and Lagrangian particle motion is governed by DEM while IMB couples the two algorithms. The new method is promising and robust as it resolves numerical instability near the particle boundary caused by mesh distortion in the conventional grid method. In IMB, the interface lattice solid fraction determines the distribution function ratio of non-equilibrium bounce back and Bhatnagar-Gross-Krook (BGK) collision. The non-equilibrium bounce back at moving boundary results in the fluid momentum change and contributes to the hydrodynamic force on particle. For numerical stability, this paper introduces the hydrodynamic force calculation concept from IB (immersed boundary method) to IMB, and at the same time, proposes a new solid fraction calculation method for sphere that divides the intersection into simple sector and triangle, as well as calculates the intersection area by vector. With this method, approximate inaccuracy is overcome while complicated integration is avoided.


Introduction
Particle flow exists in various fields, such as submarine granular flows, particle settlement, fluid beds combustion, motion of blood cells, and so on. However, capturing the micro-interaction of particles by experiment is difficult. Simulation is a remarkable technique that can well-represent a phenomenon or system. To acknowledge the hydrodynamic motion and collective behavior at grain scale, particle flow simulation is suggested.
Particle flow simulation with conventional finite element method (FEM) results in mesh distortion as the boundary of particles move. Arbitrary Lagrangian-Eulerian (ALE) and adaptive mesh refinement (AMR) techniques were proposed to reorganize the boundary grid [1]; however, FEM still suffers additional error in force coupling between fluid and particle. Since 1986, the Los Alamos lab has proposed the concept of discrete simulation of fluid dynamic [2]. A number of mesh-free methods such as smoothed particle hydrodynamics (SPH), material point method (MPM), and lattice Boltzmann method (LBM) were rapidly developed in recent years, which bypass the usual necessary grid, thereby avoiding mesh tangling and distortion.
LBM is widely used in hydrodynamic simulation, which replaces fluid macroscopic properties with distribution function, and progresses stream and collision on a background regular grid. Because moving boundary has an influence on background grid flag change, it is easy to simulate the fluid-structure interaction problem such as particle flow with a suitable coupling technique. There are two kinds of fluid-structure interaction boundaries for LBM: (a) stationary boundary, and (b) moving boundary. Ziegler [3] proposed a simple bounce back scheme of stationary boundary which turns the incident distribution function to the opposite direction. Filippova [4] first applied the lattice Boltzmann method to gasparticle flow, which allowed the flow to pass through fixed deposition of particles with irregular surface shape. For the moving boundary, three popular coupling methods exist: (a) modified bounce back method (MBB) [5], (b) immersed boundary method (IB) [6], and (c) immersed moving boundary method (IMB). In this paper, particle motion is governed by discrete element method (DEM), which is coupled with lattice Boltzmann method (LBM) to simulate particle flow.
Ladd [5] proposed MBB-LBM, introducing an additional term to the collision operator, causing momentum change of particle's neighbor fluid lattice, which then causes a reaction force on particle (hydrodynamic force) as a result of momentum change of fluid. When the boundary interface chain is not in the lattice center, the streaming point is not on the lattice. As a result, Bouzidi et al. [7] proposed a method of calculating the streaming point distribution function by spatial interpolation. Yu et al. [8] stated that momentum exchange only happens between solid interface and fluid lattice, and that force calculation only depends on neighbor fluid lattice. Soundararajan [9] combined DEM with MBB-LBM for landside simulation. In addition, Bonaccorso et al. [10] studied the rheology of bijel nano particle in confined flow. The drawback of MBB-LBM is the discontinuance of flag change of lattice, which causes complex interpolation, and fluctuation in population updating and hydrodynamic force on particle, which increases programming difficulty.
Peskin [11] suggested the immersed boundary method to simulate flow membranes interaction by using a spring between different interface chain node, and by interpolating chain node velocity from neighbor fluid lattice. Calculation of the restoring force depended on the chain node's relative displacement by Hooke's law, and the reaction boundary body force on fluid is obtained by interpolation of restoring force of boundary. Feng and Michaelides [12] combined IB and LBM to simulate the fluid-structure interaction (FSI) problem. IB-LBM was then extended to more FSI problems such as symmetric flapping wing movement simulation [13], deformable particles [14], and multicomponent flow spinodal decomposition [15]. However, IB discretized delta function interpolation can be inaccurate [16], and complicated simulations using IB can result in computational expense [6]. In addition, coupling of IB-LBM is not compatible to the parallel nature of LBM [17].
Immersed moving boundary was proposed by Noble [18]. Cook et al. [19] then coupled LBM and DEM to simulate particle flow together with IMB. Interface lattice distribution function collision is associated with solid fraction, and momentum exchange on the particle boundary is introduced to calculate the hydrodynamic force reaction on the particle. The IMB is easy to code, and is compatible with the parallel nature of LBM. In addition, the hydrodynamic force based on accurate physical description can be easily obtained using IMB. Since solid fraction is a critical variable for hydrodynamic force calculation, several papers have focused on accurate solid fraction calculation algorithm. Owen et al. [20] decomposed the lattice cell into n 2 subcells with cell size ∆x sub = ∆x/n, marking subcells' flag as solid if their centers are located within the particle radius domain. Solid fraction is calculated as the sum of all the solid subcells in the corresponding interface cell lattice. Galindo-Torres [21] proposed an edge-intersection averaging method, which calculates the intersection length of particle and interface lattice cell edge, and estimates the solid fraction of interface lattice with ratio of total interaction length and all edges length. However, edge-intersection averaging is less accurate, adding further interaction points to reconstruct convex hull [22], and solid fraction calculation depends on convex hull area. Jones and Williams [23] calculated solid fraction directly from volume integral, which is a linear expression about particle radius. However, the method only works on four and three edges interaction cases.
Based on all of the previous mentioned algorithms, this paper proposes an accurate, understandable, and easy to code solid fraction calculation scheme for sphere particles, Appl. Sci. 2021, 11, 3436 3 of 18 which is highly feasible for all kinds of interaction modes. In the present paper, the IMB-DEM-LBM, which is easy to extend to multi-particle flow, is applied to simulate the settlement and collision of two rigid particles, improving the numerical stability of solid fraction calculation. The paper is organized as follows: Section 2, Section 3, Section 4 present the theory of LBM, DEM, and IMB, respectively, while Section 5 illustrates the new solid fraction algorithm; in Section 6, the new method is verified by comparing the results from IMB-DEM-LBM with the results from analytical solution and experiments.

Incompressible Fluid Mode with D2Q9
In this section, consider a 2D incompressible fluid with standard lattice Boltzmann method (LBM). The fluid domain is divided into regular lattices and the fluid molecular velocity of each lattice is restricted to nine specific velocity vectors (see Figure 1), which is e i in Equation (1). Based on all of the previous mentioned algorithms, this paper proposes an accurate, understandable, and easy to code solid fraction calculation scheme for sphere particles, which is highly feasible for all kinds of interaction modes. In the present paper, the IMB-DEM-LBM, which is easy to extend to multi-particle flow, is applied to simulate the settlement and collision of two rigid particles, improving the numerical stability of solid fraction calculation. The paper is organized as follows: sections 2,3,4 present the theory of LBM, DEM, and IMB, respectively, while section 5 illustrates the new solid fraction algorithm; in section 6, the new method is verified by comparing the results from IMB-DEM-LBM with the results from analytical solution and experiments.

Incompressible Fluid Mode with D2Q9
In this section, consider a 2D incompressible fluid with standard lattice Boltzmann method (LBM). The fluid domain is divided into regular lattices and the fluid molecular velocity of each lattice is restricted to nine specific velocity vectors (see Figure 1), which is in Equation (1). In D2Q9 model standard LBM simulation, the lattice cell size is ∆ = 1, the time step is ∆ = 1, and the density distribution function is , which describes the total fluid molecular mass moving with specific velocity per volume. By summation, is recovered to macroscopic fluid variables, density and momentum density with Equations (2) and (3).
The time derivative of undergoes two processes. The first is streaming, in which with velocity moves to ith directional neighbor lattice cell. The second is collision. Using the second law of thermodynamics, a new set of collide towards the maximum entropy equilibrium distribution function (0) , that is the Maxwell-Boltzmann distribution, which is discretized by second order hermit polynomials corresponding to nine velocity , in which the (0) expression is given by Equation (4). In D2Q9 model standard LBM simulation, the lattice cell size is ∆x = 1, the time step is ∆t = 1, and the density distribution function is f i , which describes the total fluid molecular mass moving with specific velocity e i per volume. By summation, f i is recovered to macroscopic fluid variables, density ρ and momentum density ρu with Equations (2) and (3).
The time derivative of f i undergoes two processes. The first is streaming, in which f i with velocity e i moves to ith directional neighbor lattice cell. The second is collision. Using the second law of thermodynamics, a new set of f i collide towards the maximum entropy equilibrium distribution function f (0) , that is the Maxwell-Boltzmann distribution, which is discretized by second order hermit polynomials corresponding to nine velocity e i , in which the f (0) expression is given by Equation (4).
In collision process, distribution function transfers from chaos state to equilibrium state gradually. Bhatnagar et al. [24] proposed the most popular BGK model with single relaxation time τ, in which f i move towards equilibrium distribution function f To obtain τ, distribution function f i is transformed to the Navier-Stokes equation by Chapman-Enskog analysis, and τ is given by Equation (7), v is lattice dimensionless fluid viscosity, which is scaled from physical fluid viscosity. As necessary, viscosity must be non-negative for numerical stability, thereby τ has to satisfy τ > 0.5.
In this paper, there are two lattice Boltzmann boundary conditions: (a) non-slip resting boundary and (b) velocity boundary. The non-slip boundary rests between the fluid lattice and solid lattice, and satisfies the bounce back rule in Equation (8). f i bounces back to the opposite direction on the half way when streaming to ith directional neighbor lattice cell. In Equation (8), −i denotes the reverse direction of i while the + represents post collision, and x b is boundary lattice cells position.
Boundary with specific velocity such as inlet and outlet (see Figure 2) lacks the streaming source. The half side f i are unknown, namely; f 1 , f 5 , and f 8 . Zou and He [25] reconstructed the unknown distribution function from the channel velocity using Equations (9)- (11).
In collision process, distribution function transfers from chaos state to equilibrium state gradually. Bhatnagar et el. [24] proposed the most popular BGK model with single relaxation time , in which move towards equilibrium distribution function To obtain , distribution function is transformed to the Navier-Stokes equation by Chapman-Enskog analysis, and is given by Equation (7), is lattice dimensionless fluid viscosity, which is scaled from physical fluid viscosity. As necessary, viscosity must be non-negative for numerical stability, thereby has to satisfy > 0.5.
In this paper, there are two lattice Boltzmann boundary conditions: (a) non-slip resting boundary and (b) velocity boundary. The non-slip boundary rests between the fluid lattice and solid lattice, and satisfies the bounce back rule in Equation (8).
bounces back to the opposite direction on the half way when streaming to ith directional neighbor lattice cell. In Equation (8), − denotes the reverse direction of while the + represents post collision, and x is boundary lattice cells position.
Boundary with specific velocity such as inlet and outlet (see Figure 2) lacks the streaming source. The half side are unknown, namely; 1 , 5 , and 8 . Zou and He [25] reconstructed the unknown distribution function from the channel velocity using Equations (9)-(11).  Other boundary condition, such as periodic boundary condition, requires that outlet lattice f i streams to inlet lattice directly. Body force is a common phenomenon, such as gravity and electronic field, and in Poiseuille flow, the pressure difference is presented by body force, adding an additional term associated with body force density F to the BGK model, as shown in Equation (12).

Discrete Element Method
In this paper, particle motion and collision are simulated using the discrete element method. Lagrangian particle's motion follows Newton's second law in Equations (14) and (15). Particle translation is accelerated by hydrodynamic interaction force F f , particles collision force F c , and gravity force G every timestep. Particle rotational acceleration depends on torque M f and M c , derived from F f and F c .
In Equations (14) and (15), U i and Ω i are translational and rotational velocity, respectively; m and I are mass and moment of inertia, respectively; F f is hydrodynamic force on particle calculated by immersed moving boundary, which will be discussed in the next section. x i and x j are two particles' center vector, r i and r j are two particles' radius. v i and v j are two particles' velocity. The collision surface normal and tangential vector of particle i are n = x i − x j / x j − x i and t, respectively. Collision force in the normal and tangential direction is F n and F t . F n and F t are divided into linear elastic force and damping force. For F n , linear elastic force is related to overlap of collision particle δ ij = r i + r j − x j − x i . Damping force is related to the particle's relative velocity δ ij )·n (16) where constant k n is normal stiffness coefficient. constant γ n is normal damping coefficient. Tangential force F t is caused by the rough surface sliding between two particles. Previous time step tangential overlap is forgotten in the next time step. The relative tangential velocity is δ t ij ·dt. The tangential stiffness coefficient is k t = 2 7 k n and the tangential damping coefficient is [26]. The particle i tangential force is given by Equation (17).
The tangential force is limited by the Coulomb's friction law: |F t | < µF n , thereby the final tangential force expression is given by where µ is friction coefficient.

Immersed Moving Boundary for LBM-DEM Coupling
Different from MMB and IB, in IMB, distribution function in the non-conforming boundary lattice is divided into two parts. The first part undergoes non-equilibrium bounce back while the second part collides like normal fluid lattice. The ratio of two parts depends on solid fraction of interface lattice. In Figure 3, cells are divided into three types, and the curve is the particle boundary interface. Rectangle symbol lattice is particle cell, sphere symbol lattice is fluid cell, and pentagon symbol lattice is interface cell. The shaded section represents the fluid part of interface cell.
where is friction coefficient.

Immersed Moving Boundary for LBM-DEM Coupling
Different from MMB and IB, in IMB, distribution function in the non-conforming boundary lattice is divided into two parts. The first part undergoes non-equilibrium bounce back while the second part collides like normal fluid lattice. The ratio of two parts depends on solid fraction of interface lattice. In Figure 3, cells are divided into three types, and the curve is the particle boundary interface. Rectangle symbol lattice is particle cell, sphere symbol lattice is fluid cell, and pentagon symbol lattice is interface cell. The shaded section represents the fluid part of interface cell. Solid fraction of interface lattice cell is , which varies from 0 to 1. This paper proposes a new calculation algorithm for single particle, which will be shown in Section 5. In multi particle collision simulation, solid fraction is calculated by summing all of the overlap area of every particle. For numerical stability, will be replaced with weight function in Equation (19), in which the two parameters have the same physical meaning. The relationship between and is shown in Figure 4, which shows that changes more slowly than when the fluid cell's flag changes, providing a smoother transition of coupling force and a reduction of force fluctuation. Solid fraction of interface lattice cell is ε s , which varies from 0 to 1. This paper proposes a new ε s calculation algorithm for single particle, which will be shown in Section 5. In multi particle collision simulation, solid fraction is calculated by summing all of the overlap area of every particle. For numerical stability, ε s will be replaced with weight function B s in Equation (19), in which the two parameters have the same physical meaning. The relationship between B s and ε s is shown in Figure 4, which shows that B s changes more slowly than ε s when the fluid cell's flag changes, providing a smoother transition of coupling force and a reduction of force fluctuation.  After streaming, interface lattice obtains a new set of from its fluid neighbor lattice as the interface lattice is mixed components including fluid and solid. With the different solid fraction, the corresponding solid part distribution function will bounce back directly, which satisfies the non-equilibrium bounce back rule. If the boundary velocity is u , and the boundary fluid velocity is u, define ( , u ) and ( , u) as the boundary lattice equilibrium distribution function of solid and fluid parts. The non-equilibrium part of distribution function is − (x, ) − − ( , u), which will bounce back to i direction. Con- After streaming, interface lattice obtains a new set of f i from its fluid neighbor lattice as the interface lattice is mixed components including fluid and solid. With the different solid fraction, the corresponding solid part distribution function B s f i will bounce back directly, which satisfies the non-equilibrium bounce back rule. If the boundary velocity is u s , and the boundary fluid velocity is u, define f eq i (ρ, u s ) and f eq i (ρ, u) as the boundary lattice equilibrium distribution function of solid and fluid parts. The non-equilibrium part of distribution function is f −i (x, t) − f eq −i (ρ, u), which will bounce back to i direction. Considering the value change from the origin value f i (x, t), the non-equilibrium bounce back expression is given by Equation (20).
The fluid parts distribution function (1 − B s ) f i , which undergoes BGK collision, is shown in Equation (21).
Combining Equations (20) and (21), an additional term is added to collision operator, and the mix collision expression of interface lattice is given in Equation (22).
conventional collision operator has no effect on momentum change. New term Ω s i introduces the momentum change of fluid in the interface lattice cell. The corresponding reaction force F f on the particle is −B n (∑ i Ω s i e i ). For numerical stability, this paper introduces conception fluid lattice indicator w from MBB method [8]. If the corresponding neighbor cell is fluid, 1 − w = 1, and on the contrary, 1 − w = 0 (see Figure 5). The new hydrodynamic force expression is given by Equation (24).

New Solid Fraction Calculation Algorithm
Solid fraction can be seen as the intersection area of particle geometry and LBM grid. There are three intersection states between particle and lattice, and the relationship between lattice position and particle center introduces more states for solid fraction calculation, which must be considered separately. Thus, it is difficult to calculate solid fraction

New Solid Fraction Calculation Algorithm
Solid fraction can be seen as the intersection area of particle geometry and LBM grid. There are three intersection states between particle and lattice, and the relationship between lattice position and particle center introduces more states for solid fraction calculation, which must be considered separately. Thus, it is difficult to calculate solid fraction cell by cell with different scheme. In the conventional method, intersection area is estimated by approximation of integral because it is difficult to calculate cell by cell.
In the new solid fraction calculation algorithm, the intersection is not calculated by the local cell intersection area, and the calculation progress is taken as a whole, wherein the intersection area depends on the other pre-calculating cells. Before calculation, the possible intersect lattice cell is marked based on the particle center and radius (see Figure 6a), and the marked cell is divided into four parts (A, B, C, D). Corresponding cell index (i, j) is given in each part, and the index origin of each part begins from the farthest cell from the particle center. i and j indicate the row and column (see Figure 6b), and the corresponding intersection area of (i, j) cell is S i,j . , is calculated individually for each of the lattice groups A, B, C, D in an ascending order of i and j (see Figure 7a), in which , depends on all the pre-calculated , , where = 1 ⋯ , = 1 ⋯ . Since , is difficult to calculate directly, a trick is used in which a new variable , is defined, as shown in Equation (25). Equation (25) calculates the sum of all the pre-calculated , , and then all the known , (not include , ) are subtracted from , to get , .
An example is shown in Figure 7b. To calculate , , it is assumed that all the precalculated The whole progress for coding is listed as follows: 1. At every time step, the particle's new position is located, and the intersection lattice group is updated depending on the particle location. Choices of lattices for solid fraction calculation depend on the particle's center lattice location, and the lattices' type and index are updated every time step. 2. Loop A, B, C, D into four parts. S i,j is calculated individually for each of the lattice groups A, B, C, D in an ascending order of i and j (see Figure 7a), in which S i,j depends on all the pre-calculated S m,n , where m = 1 · · · i, n = 1 · · · j. Since S i,j is difficult to calculate directly, a trick is used in which a new variable S pre i,j is defined, as shown in Equation (25). Equation (25) calculates the sum of all the pre-calculated S m,n , and then all the known S m,n (not include S i,j ) are subtracted from S pre i,j to get S i,j .
An example is shown in Figure 7b. To calculate S i,j , it is assumed that all the precalculated S m,n are already known, and only S i−1,j and S i,j−1 are non-zero. S At every time step, the particle's new position is located, and the intersection lattice group is updated depending on the particle location. Choices of lattices for solid fraction calculation depend on the particle's center lattice location, and the lattices' type and index are updated every time step.

2.
Loop A, B, C, D into four parts.

3.
Loop all cells in ascending order of i and j.

4.
Calculate S pre i,j , then subtract all the known S m,n to obtain S i,j .

Simulation
The simulation results of three cases are given in this section. The first case is plane Poiseuille flow, a flow between plane with body force. By comparing the proposed solution and analytical solution, the accuracy of LBM in fluid simulation is verified, which is the foundation in getting an accurate hydrodynamic force on particle.
The second case is drag force on a fixed cylinder. The application of plane flow generates laminar flow with different Reynolds number, in which the corresponding analytical drag coefficient can be calculated. Thereafter, the hydrodynamic force with IMB-LBM is simulated. The simulation and analytical drag coefficient value are compared to illustrate the hydrodynamic force accuracy on static sphere.
The third case is the simulation of settlement of a single particle under the influence of gravity, in which the drag force and the velocity increase as the particle sinks. In this case, the final velocity reaches equilibrium when the drag force is equal to the gravity force.
The fourth case is the simulation of two particles falling freely in a container, releasing two particles in the center line with body force to be able to verify the particle motion: "drafting, kissing, and tumbling" phenomenon.

Poiseuille Flow
Before real particle flow of IMB-DEM-LBM simulation, LBM is applied to Poiseuille flow. The top and bottom boundary are fixed boundary, and the left and right boundary are periodic boundary. Pressure at inlet and outlet accelerate incompressible fluid flow through two plain plates (see Figure 8).

Simulation
The simulation results of three cases are given in this section. The first case is plane Poiseuille flow, a flow between plane with body force. By comparing the proposed solution and analytical solution, the accuracy of LBM in fluid simulation is verified, which is the foundation in getting an accurate hydrodynamic force on particle.
The second case is drag force on a fixed cylinder. The application of plane flow generates laminar flow with different Reynolds number, in which the corresponding analytical drag coefficient can be calculated. Thereafter, the hydrodynamic force with IMB-LBM is simulated. The simulation and analytical drag coefficient value are compared to illustrate the hydrodynamic force accuracy on static sphere.
The third case is the simulation of settlement of a single particle under the influence of gravity, in which the drag force and the velocity increase as the particle sinks. In this case, the final velocity reaches equilibrium when the drag force is equal to the gravity force.
The fourth case is the simulation of two particles falling freely in a container, releasing two particles in the center line with body force to be able to verify the particle motion: "drafting, kissing, and tumbling" phenomenon.

Poiseuille Flow
Before real particle flow of IMB-DEM-LBM simulation, LBM is applied to Poiseuille flow. The top and bottom boundary are fixed boundary, and the left and right boundary are periodic boundary. Pressure at inlet and outlet accelerate incompressible fluid flow through two plain plates (see Figure 8).

Poiseuille Flow
Before real particle flow of IMB-DEM-LBM simulation, LBM is applied to Poiseuille flow. The top and bottom boundary are fixed boundary, and the left and right boundary are periodic boundary. Pressure at inlet and outlet accelerate incompressible fluid flow through two plain plates (see Figure 8). The distance between two plates is H, and fluid viscosity is . Pressure drop is replaced with body force F in Equation (26). The distance between two plates is H, and fluid viscosity is υ. Pressure drop is replaced with body force F b in Equation (26).
The corresponding velocity profile analytical solution resolved from Navier-Stokes equation with boundary condition u(0) = u(H) = 0 is given by Equation (27).
The simulation domain cell number is 60 × 60, and the dimensionless variables are listed as follows: (a) collision time τ = 1, (b) dimensionless viscosity ν dimensionless = 0.167, (c) max dimensionless velocity in x direction is 0.1, (d) velocity in y direction is 0, (e) cell size ∆x = 1, and (f) timestep ∆t = 1. Steady simulation result is obtained at 300 time steps. The velocity profile at x = 30 using the coarse grid solution and analytical solution is shown in Figure 9. Based on the result, the standard error estimate σ est = ∑ (u simulation −u analytical ) 2 N , is σ est = 0.0012. By comparing the simulation result and analytical result, the accuracy of LBM fluid simulation is verified, which is the foundation in obtaining a more accurate hydrodynamic force on particle.

Drag Force on Fix Cylinder
The accuracy of hydrodynamic force on particle is critical for IMB-DEM-LBM particle flow simulation. A fixed cylinder is located in a channel with fixed boundary velocity (see Figure 10). Different inlet velocity generates steady laminar flow with different Reynolds number. The cylinder drag coefficient is expressed analytically by a function of Reynolds number, which is shown in Equation (28) [27,28].

Drag Force on Fix Cylinder
The accuracy of hydrodynamic force on particle is critical for IMB-DEM-LBM particle flow simulation. A fixed cylinder is located in a channel with fixed boundary velocity (see Figure 10). Different inlet velocity generates steady laminar flow with different Reynolds number. The cylinder drag coefficient C d is expressed analytically by a function of Reynolds number, which is shown in Equation (28) [27,28].

Drag Force on Fix Cylinder
The accuracy of hydrodynamic force on particle is critical for IMB-DEM-LBM particle flow simulation. A fixed cylinder is located in a channel with fixed boundary velocity (see Figure 10). Different inlet velocity generates steady laminar flow with different Reynolds number. The cylinder drag coefficient is expressed analytically by a function of Reynolds number, which is shown in Equation (28) [27,28].  The inlet velocity is simulated with Zou-He boundary condition. The drag force on cylinder is F simulation , and the fluid velocity at cylinder is U simulation . The drag coefficient is shown in Equation (29).
The  Figure 11. The simulation force on cylinder is obtained, and then drag coefficient is calculated, as well as the error = (simulation − analysis)/analysis, as shown in Table 1. A comparison between the analytical solution and the proposed solution is shown in Figure 12. As shown, the simulation result agrees with the analytical result. In addition, the simulation has a better performance at high Reynolds number. The simulation error at low Reynolds number is caused by velocity fluctuation near cylinder boundary.
The inlet velocity is simulated with Zou-He boundary condition. The drag force on cylinder is , and the fluid velocity at cylinder is . The drag coefficient is shown in Equation (29).  Figure 11. The simulation force on cylinder is obtained, and then drag coefficient is calculated, as well as the = ( − )/ , as shown in Table  1. A comparison between the analytical solution and the proposed solution is shown in Figure 12. As shown, the simulation result agrees with the analytical result. In addition, the simulation has a better performance at high Reynolds number. The simulation error at low Reynolds number is caused by velocity fluctuation near cylinder boundary.     Error at low Reynolds number is introduced by numerical noise. In low Reynolds number flow, the interaction drag force is not dominant but it helps match the particle and fluid velocity. The particle velocity oscillation is inevitable in the simulation, which is different from the actual physical progress. Using a large drag force helps in decreasing oscillation. This shows the limitation of the coupling method presented in this study; however, this error can be neglected because low Reynolds number condition is rare for particle flow in coastal engineering.

Single Particle Settlement
To demonstrate the feasibility of IMB-DEM-LBM for particle movement, a particle that falls in a container filled with liquid under the effect of gravity is simulated to illustrate velocity accuracy of a moving particle.
The square container's physical size 6 × 10 −4 m, and the domain boundary is no slip bounce back boundary. The gravity acceleration is = 9.8 m s 2 ⁄ , particle radius is = 5 × 10 −5 m, particle density is = 2000 kg m 3 ⁄ , fluid density is = 1000 kg m 3 ⁄ , dynamic fluid viscosity is = 10 −3 Pa • s, conversion factor = 5 × 10 −5 , = 1, the lattice domain size is 60 × 60, and particle radius is 5. The initial particle state is stationary and the fluid velocity is zero. The particle position and velocity profile during settlement is shown in Figure 13. Error at low Reynolds number is introduced by numerical noise. In low Reynolds number flow, the interaction drag force is not dominant but it helps match the particle and fluid velocity. The particle velocity oscillation is inevitable in the simulation, which is different from the actual physical progress. Using a large drag force helps in decreasing oscillation. This shows the limitation of the coupling method presented in this study; however, this error can be neglected because low Reynolds number condition is rare for particle flow in coastal engineering.

Single Particle Settlement
To demonstrate the feasibility of IMB-DEM-LBM for particle movement, a particle that falls in a container filled with liquid under the effect of gravity is simulated to illustrate velocity accuracy of a moving particle.
The square container's physical size 6 × 10 −4 m, and the domain boundary is no slip bounce back boundary. The gravity acceleration is g = 9.8 m/s 2 , particle radius is r = 5 × 10 −5 m, particle density is ρ P = 2000 kg/m 3 , fluid density is ρ F = 1000 kg/m 3 , dynamic fluid viscosity is η = 10 −3 Pa·s, conversion factor C l = 5 × 10 −5 , τ = 1, the lattice domain size is 60 × 60, and particle radius is 5. The initial particle state is stationary and the fluid velocity is zero. The particle position and velocity profile during settlement is shown in Figure 13. The drag force on particle increases as the sinking velocity increases, and the particle velocity reaches equilibrium when the drag force is equal to the gravity force. The analytical equilibrium sinking velocity is given by Equation (30) [29]. Therefore, the analytical sinking equilibrium velocity is = 5.45 × 10 −3 m s ⁄ . The simulation result using the method proposed in this study is shown in Figure 14, which indicates an equilibrium velocity of about 5.75 × 10 −3 m/s. The drag force on particle increases as the sinking velocity increases, and the particle velocity reaches equilibrium when the drag force is equal to the gravity force. The analytical equilibrium sinking velocity is given by Equation (30) [29]. Therefore, the analytical sinking equilibrium velocity is u eq = 5.45 × 10 −3 m/s. The simulation result using the method proposed in this study is shown in Figure 14, which indicates an equilibrium velocity of about 5.75 × 10 −3 m/s.

Particles Settlement and Collision
To demonstrate the feasibility of IMB-DEM-LBM for particle collision in liquid, settlement of two particles in the same path is simulated. When two departed particles realized in the container center line, the particles' initial velocity is zero and both of th start to move due to gravity. The upper particle is in the wake of the lower particle, a is accelerated to catch up with the lower particle [30]. They will then kiss together a settle together as a whole. Little fluctuation in the flow may make them depart from e other. The total progress is called "drafting, kissing, and tumbling" or DKT, in which s ulation of this phenomenon can be done by different techniques such as FEM [31,32], LBM [12]. "Drafting, kissing and tumbling" is a random phenomenon, which occur the settlement and collision of particles. In experiments, the time interval of the kiss progress depends on when the random disturbance happens. Unfortunately, the exp ment cannot catch how and when the disturbance occurs. Hence, the purpose of the s ulation is only to show that these phenomena can be replicated or depicted by the meth presented in this study.
For saving computational power, the two particles' physical initial velocity −0.09 m s ⁄ . The initial velocity allows particles to reach balance velocity with less ti Physical domain is 0.25 × 1 mm, gravity is = 9.8 m s 2 ⁄ , particle radius is 2.5 × 10 −5 m , fluid density is = 2000 kg m 3 ⁄ , and particles mass is 8 × 10 −5 kg. The distance of two particles' center is 0.125 × 10 −5 m. The collision tim = 1, and the lattice domain size is 200 × 800. The particles' position and velocity pro is shown in Figure 15a

Particles Settlement and Collision
To demonstrate the feasibility of IMB-DEM-LBM for particle collision in liquid, the settlement of two particles in the same path is simulated. When two departed particles are realized in the container center line, the particles' initial velocity is zero and both of them start to move due to gravity. The upper particle is in the wake of the lower particle, and is accelerated to catch up with the lower particle [30]. They will then kiss together and settle together as a whole. Little fluctuation in the flow may make them depart from each other. The total progress is called "drafting, kissing, and tumbling" or DKT, in which simulation of this phenomenon can be done by different techniques such as FEM [31,32], IB-LBM [12]. "Drafting, kissing and tumbling" is a random phenomenon, which occurs in the settlement and collision of particles. In experiments, the time interval of the kissing progress depends on when the random disturbance happens. Unfortunately, the experiment cannot catch how and when the disturbance occurs. Hence, the purpose of the simulation is only to show that these phenomena can be replicated or depicted by the method presented in this study.
For saving computational power, the two particles' physical initial velocity is −0.09 m/s. The initial velocity allows particles to reach balance velocity with less time. Physical domain is 0.25 × 1 mm, gravity is g = 9.8 m/s 2 , particle radius is r = 2.5 × 10 −5 m, fluid density is ρ P = 2000 kg/m 3 , and particles mass is mass particle = 8 × 10 −5 kg. The distance of two particles' center is 0.125 × 10 −5 m. The collision time is τ = 1, and the lattice domain size is 200 × 800. The particles' position and velocity profile is shown in Figure 15a-h for the corresponding time steps at t = 0, 40, 80, 130, 190, 230, 260, 300.
In Figure 15a, the two particles begin to settle together, and the particles get closer and closer in Figure 15b,c. In Figure 15d, the two particles kiss each other, and begin to drift together. In Figure 15d,e, the two particles adhere together. In Figure 15f,h, little fluctuation makes the particles depart from each other. The particles' motion in Figure 14 satisfies "drafting, kissing, and tumbling" phenomenon, which occurred in the experiment conducted by Fortes [33].
(e) (f) (g) (h) Figure 15. Two particles settling, drafting, kissing, and tumbling Figure 15. Two particles settling, drafting, kissing, and tumbling. Figure 16 shows the x and y velocity of the two settling particles. In Figure 16b, the whole progress can be divided into three time periods. At timestep 0-75, the upper particle chases the lower one, and at timestep 75-190, the two particles drift together and collide repeatedly. At timestep 190-300, the two particles depart from each other. In drift progress, gravity potential is converted into x momentum, as shown in Figure 16a, and when two particles depart from each other, the x coordinate momentum is exchanged to the fluid around the particles, so all the particles' x velocity decreases, as shown in Figure 16a.
In Figure 15a, the two particles begin to settle together, and the particles get closer and closer in Figure 15b and Figure 15c. In Figure 15d, the two particles kiss each other, and begin to drift together. In Figure 15d and Figure 15e, the two particles adhere together. In Figure 15f and Figure 15h, little fluctuation makes the particles depart from each other. The particles' motion in Figure 14 satisfies "drafting, kissing, and tumbling" phenomenon, which occurred in the experiment conducted by Fortes [33]. Figure 16 shows the x and y velocity of the two settling particles. In Figure 16b, the whole progress can be divided into three time periods. At timestep 0-75, the upper particle chases the lower one, and at timestep 75-190, the two particles drift together and collide repeatedly. At timestep 190-300, the two particles depart from each other. In drift progress, gravity potential is converted into x momentum, as shown in Figure 16a, and when two particles depart from each other, the x coordinate momentum is exchanged to the fluid around the particles, so all the particles' x velocity decreases, as shown in Figure  16a.

Conclusions
This paper presents a promising procedure of IMB-DEM-LBM for particle flow simulation of rigid particles such as soil particles, in which IMB and IB are combined to get the new hydrodynamic force calculation expression. To improve the accuracy of solid fraction calculation, a new method, which can easily be understood and is easy to program, was proposed based on the basic geometric relationship for solid fraction calculation. Three cases were simulated, which includes theoretical comparison and experiment phenomenon verification. The first case was Poiseuille flow, which is a routine test for hydrodynamic simulation. The second case is drag coefficient simulation, which was performed to test the accuracy of hydrodynamic force on particle. By comparing the analytical solution and proposed solution, the new hydrodynamic force expression is verified, which makes the new solid fraction calculation algorithm suitable, feasible, and accurate for IMB. The third case is particles collision, which was performed to extensively test the movement of particles. The simulation results coincide with an experiment conducted in the literature, verifying the new method for more complex boundary conditions, as well as providing steady numerical solution for quick flag change problem. Based on a rigorous theoretical derivation and example verification, the new method presents an outstanding performance on particle-flow simulation. The simulation conducted in this paper is two-dimensional; however, IMB, DEM, and LBM can be easily extended to three dimensions (see Appendix A). The new solid fraction method can be extended to threedimensional problems by replacing the two-dimensional lattice index with three-dimensional lattice index, in which the other conception is the same as the two-dimensional case.

Conclusions
This paper presents a promising procedure of IMB-DEM-LBM for particle flow simulation of rigid particles such as soil particles, in which IMB and IB are combined to get the new hydrodynamic force calculation expression. To improve the accuracy of solid fraction calculation, a new method, which can easily be understood and is easy to program, was proposed based on the basic geometric relationship for solid fraction calculation. Three cases were simulated, which includes theoretical comparison and experiment phenomenon verification. The first case was Poiseuille flow, which is a routine test for hydrodynamic simulation. The second case is drag coefficient simulation, which was performed to test the accuracy of hydrodynamic force on particle. By comparing the analytical solution and proposed solution, the new hydrodynamic force expression is verified, which makes the new solid fraction calculation algorithm suitable, feasible, and accurate for IMB. The third case is particles collision, which was performed to extensively test the movement of particles. The simulation results coincide with an experiment conducted in the literature, verifying the new method for more complex boundary conditions, as well as providing steady numerical solution for quick flag change problem. Based on a rigorous theoretical derivation and example verification, the new method presents an outstanding performance on particleflow simulation. The simulation conducted in this paper is two-dimensional; however, IMB, DEM, and LBM can be easily extended to three dimensions (see Appendix A). The new solid fraction method can be extended to three-dimensional problems by replacing the two-dimensional lattice index with three-dimensional lattice index, in which the other conception is the same as the two-dimensional case.