Pore-Scale Simulations of Particles Migration and Deposition in Porous Media Using LBM-DEM Coupling Method

This paper studies the migration and deposition of suspended particles in porous media. This problem results from the fact that during the operation of a groundwater source heat pump, the recharging process will contribute to the impairment of soil permeability. A coupling lattice Boltzmann method, discrete element method and immersed moving boundary method were used to investigate the migration of particles in porous media. The DKT (Drifting, Kissing, Tumbling) phenomena were employed to validate our program. The coupled effects of concentration, flow rate and pH on the clogging mechanism of the porous media were analyzed. Results show that, due to the repulsive barrier between the particles and porous media, there is a critical velocity. At a low flow rate, the deposition ratio increases with the increase in velocity. Beyond the critical velocity, the deposition ratio decreases when the velocity increases due to higher shear force. Permeability impairment increases with the increase in concentration, especially in the low flow rate condition. Changes in pH mainly affect the repulsive barrier. For a low flow rate, the decrease in repulsive barrier greatly promotes the deposition of particles. Under the condition of favorable deposition, the increase in flow rate reduces the deposition phenomenon. Under the condition of unfavorable deposition, the lower flow rate condition has a lower deposition ratio. The process of particle deposition and the dynamic motion after deposition were observed such as particles gliding over the surface. Accumulated particles in the downstream form bridges and hinder fluid flow. At a high flow rate, strong shear force is more capable of destroying bridges and recovering permeability. Adsorbed particles glide on the surface of the grain and deposit in the downstream. This paper aims to help understanding of the micro-events of particle deposition and the clogging process.


Introduction
Suspended particles migration and deposition in porous media is a common issue in natural phenomena and industry processes. During the operation of the recharge system of a groundwater source heat pump system (GSHP), the problem of internal clogging in the aquifer reduces the efficiency of the recharge system, thus affecting the operation of the GSHP [1,2]. Particles invade the flow channel or pore space under the action of hydrodynamic force and strain force, resulting in reduced porosity and permeability impairment [3][4][5]. Such problems also occur in other fields including treatment of septic tank effluent and wastewater [6,7] and development of oil and gas fields [8,9]. Therefore, understanding the migration and deposition mechanism of particles in porous media and predicting the clogging formation have important academic value and practical significance.
Some studies have been carried out to tackle this problem and put forward several mathematical models. To study particle migration in porous media, some theoretical models has been developed [10,11]. However, estimating permeability value as a function of porosity, grain size, and the structure of pore space can be challenging, which limits the application of theoretical models. At present, the mainstream method to is using column experiments. Alem et al. studied the effect of hydrodynamic forces on particle migration in porous media [12]. The study shows that as the injection volume increases, the larger particles penetrate more deeply into the porous media, and the size of particles observed in the effluent gradually increases. Zhang et al. studied the influence of pore structure and flow rate on the hydrodynamic mechanism, diffusion effect and acceleration effect of particle migration and deposition in saturated porous media [13]. The result shows that with the increase in seepage velocity, the influence of hydrodynamics on the outflow concentration of particles is increasing, while the influence of pore structure is relatively weakened. In addition, several studies have been conducted about the mechanism of suspended particle transportation by experiments [14][15][16]. These column experiments provide important understanding of particle transport and deposition characteristics. Based on experiments, the process of particles deposition in porous media has been detailed before [17][18][19].
To observe the micromechanical process of particle migration, deposition and retention in experiment methods and classical modeling is still a challenge. Although, some microscopic visualization experiment techniques have had rapid development such as particle image velocimetry (PIV) [20] and magnetic resonance velocimetry (MRV) [21], it is still hard to reveal the mechanical properties of interaction between fluid and solid phase. As an alternative, numerical simulations provide a potential way to understand these processes. Computational fluid dynamics (CFD) coupled with the Discrete element method (DEM) method has been an effective tool for studying particles transportation such as fluidized bed system [22,23], particulate fouling [24,25] and some other fields [26][27][28]. By using the coupled method, the local phase velocity, spatial distribution, detailed information flow pattern of liquid phasor, particle residence time and particle trajectory can be studied. Benefiting from simplified governing equations and parallel computing performance, the Lattice Boltzmann method (LBM) coupled with DEM is more efficient than the CFD-DEM method when dealing with moving particle boundaries. Based on this method, some researches have been conducted including particle transport in flows [29,30], solid rheology [31,32] and soil fluidization [33,34]. Actually, there is seldom literature about observing particle migration behavior and clogging formation process in packed granular materials at pore scale. Parvan studied particle retention and clogging in porous media [35]. The result shows that the structure of porous media has great influence on particle migration. Zhou studied particle retention and permeability impairment in porous media [36]. They found that the size of solid particles, flow rate and particle volume fraction have great influence on permeability impairment of porous media. This paper mainly studies the migration and deposition process of particle in porous media and the permeability impairment of porous media. Since fluid-solid interactions in porous media flow occur in pore scale, it seems more appropriate to simulate this process from a microscopic or mesoscopic perspective. Therefore, we attempt to couple the LBM-DEM to further understand the mechanism of particle migration and deposition in pore space. This study provides a microprocess observation and understanding process of GSHP recharge clogging.

Materials and Methods
A modified LBM under a Eulerian coordinate is used to simulate the fluid flow and a DEM soft-sphere model under a Lagrange coordinate is used to deal with the collision problems between solid phase. When dealing with a fluid-solid boundary, the immersed moving boundary method (IMB) is chosen to couple the LBM-DEM model [37]. The model considers the two-way coupled with solid phase and liquid phase and can capture the motion of particle and track the trajectory of individual particle.

Particle Governing Equations Based on DEM
DEM is a computer simulation method in which particles are specified as individual physical and mechanical properties first introduced by Cundall and Strack [38]. The motion Processes 2021, 9,465 3 of 23 of particles in DEM is determined by Newton's second law. Each particle has two types of motion: translation and rotation. The motion of each particle is governed by the law of conservation of linear momentum and angular momentum. The equation of motion is numerically integrated by the finite difference method: where i denotes the particle index; x denotes the position of particle; v denotes the velocity of particle; ω denotes the angular velocity of particle; m denotes the particle mass; I denotes the moment of inertia and ∆t denotes time step. Each particle trajectory is explicitly solved by the following force balance: T where F f and T f denotes the hydrodynamic force and torque induced by fluid flow, respectively; F c and T c denotes the contact forces and the torque resulting from the contact forces, respectively; F ad and F re denotes the adhesive force and repulsive force, separately; F g denotes the gravity force.
In this paper, we calculate the contact force between particles by the soft-sphere DEM model. Soft-sphere DEM simplify the particle interaction into a spring damper process [39]. The interaction force of each pair of contact particles can be expressed as: F n c = − k n δ n + η n V n ij n ij n ij (5) where n and t denote the normal and tangent unit vectors, respectively; F c n and F c t denotes the normal and tangent contact forces, respectively. As Figure 1 shows, when two particles collide, the deformation of the interface is replaced by two overlaps. δ n and δ t denotes the normal and tangent overlap length, respectively; V ijn and V ijt denote the relative velocity at the normal and tangent contact point, respectively; k n and k t denotes normal and tangent stiffness constant, respectively. η n and η t denote normal damping constant and tangent damping constant, respectively and ζ denotes the friction coefficient. The spring coefficients are related to the Young's modulus and Poisson's ratio [40]. Additionally, η is features the restitution coefficient σ and the reduced mass m ij = m i ·m j /(m i + m j ), and in our simulations, we always choose η n = η t : In the static friction stage, the tangential force is derived from the tangential displacement and damping of the particles. When the tangential force crosses the threshold, the friction turns to dynamic friction. At this point the tangential force is a function of the absolute value of the normal force. In addition, the model considering rolling friction, additional torque is added according to the following equations: The spring in the normal and tangential direction accounts for the elastic part of the impact force, while the dashpot accounts for the energy loss during the collision. The duration of the collision event can be calculated by the following [41].
When the contact between particles occurs, the force and torque acting on the particles are calculated by Equations (4)- (8).
Processes 2021, 9, x FOR PEER REVIEW 4 of 23 The spring in the normal and tangential direction accounts for the elastic part of the impact force, while the dashpot accounts for the energy loss during the collision. The duration of the collision event can be calculated by the following [41].
When the contact between particles occurs, the force and torque acting on the particles are calculated by Equations (4)-(8).

Long-Range Force
The long-range force is mainly calculated by DLVO theory in this paper [42]. Van der Waals force and electric double-layer force are considered to be the dominant forces in DLVO theory [43]. Van der Waals interaction is a microscopic scale attraction originating from permanent and induced electrical interactions between two or more atoms or molecules. Van der Waals force between particles can be calculated by Hamaker expression [44]: where Sij denotes the center-to-center distance between the two particles; AH is the Hamaker constant. As can be seen from the formula, Van der Waals energy is generally only related to the size of the particles and the distance between the particles and the surface of the medium. The major repulsive force between particles is electric double-layer force. Analysis expression of force of electric double-layer can be obtained by Debye approximation method: where εr denotes the permittivity of vacuum; ε0 denotes the dielectric constant of the medium; R denotes the gas constant; F denotes the Faraday constant; ψ denotes the zetapotentials of the surface; z denotes the valence of the background electrolyte; hij denotes the distance between particle surfaces and κ is the reciprocal of Debye length, the expression is as follow:

Long-Range Force
The long-range force is mainly calculated by DLVO theory in this paper [42]. Van der Waals force and electric double-layer force are considered to be the dominant forces in DLVO theory [43]. Van der Waals interaction is a microscopic scale attraction originating from permanent and induced electrical interactions between two or more atoms or molecules. Van der Waals force between particles can be calculated by Hamaker expression [44]: where S ij denotes the center-to-center distance between the two particles; A H is the Hamaker constant. As can be seen from the formula, Van der Waals energy is generally only related to the size of the particles and the distance between the particles and the surface of the medium. The major repulsive force between particles is electric double-layer force. Analysis expression of force of electric double-layer can be obtained by Debye approximation method: where ε r denotes the permittivity of vacuum; ε 0 denotes the dielectric constant of the medium; R denotes the gas constant; F denotes the Faraday constant; ψ denotes the zetapotentials of the surface; z denotes the valence of the background electrolyte; h ij denotes the distance between particle surfaces and κ is the reciprocal of Debye length, the expression is as follow: Processes 2021, 9, 465 5 of 23 where e l denotes the electronic basic charge; N A denotes the Avogadro constant and I denotes the ionic strength of solution. DLVO theory is based on the direct collision hypothesis, and does not consider the motion and force caused by collision. For particle size is higher than micron, the force on particle collision cannot be ignored [45]. When S ij approaches zero, there will be numerical singularity problems in Equation (11), which is physically and practically unrealistic. Some research solves the problem by using a cut-off distance (0.1-0.2 nm) between the surfaces of a given particle [46,47]. When simulate adhesive particles contact, the JKR (Johnson, Kendall and Roberts) theory give good predictions of the contact size [48]. In this paper, we use JKR equation to calculated adhesive force F ad when two particles come close to physically contact each other [49]: where w ilj is surface adhesive energy determined by interfacial energy.

Fluid Governing Equations Based on LBM
The Lattice Boltzmann Method (LBM) is a fluid dynamics numerical solution based on the molecular motion theory and statistical mechanics [50]. In LBM, the flow field of fluid is divided into regular grids, and the fluid is regarded as a group of fluid particles at the grid nodes. The flow process of the fluid can be obtained by using the lattice Boltzmann equation to solve the collision function and the migration function of each grid node. The Lattice Boltzmann equation has two main models used to deal with the collision process between fluid particles: Bhatnagar-Gross-Krook model (BGK) based on the theory of single relaxation time and Multiple Relaxation Time model (MRT). Due to the good computational efficiency of BGK, in this paper we use BGK to deal with the collision process. In BGK, the Boltzmann equation reads as: where c i denotes the discrete lattice velocity; f i (x, t) denotes the density distribution function for a fluid node moving with velocity c i at position x and time t; ∆t denotes the time step; τ BGK denotes the relaxation time; f eq i (x, t) denotes the local equilibrium distribution function of f i (x, t), It is usually expressed as a discrete Maxwell distribution function: where w i denotes weight factor; ρ denotes fluid density; u denotes fluid velocity; c s denotes lattice sound velocity. DnQm can be used to divide BGK models with different discrete speeds in the division of grid. Where n denotes the spatial dimension and m denotes the number of discrete velocity directions. In this paper, we use D2Q9 model. In this model c s = c/ √ 3, where c = ∆x/∆t; ∆x denotes the lattice distance; ∆t denotes the time step. In the model, the 9 direction velocity vectors and weighing factors can be expressed as: Processes 2021, 9, 465 6 of 23 According to the above formula, the density distribution function at the lattice point can be obtained. The density ρ, velocity u and pressure P of the fluid can be calculated by the density distribution function: In addition, the kinematic viscosity of the fluid ν is expressed as:

The LBM-DEM Coupled Method
In dealing with the interactions between fluid flow and particle movement, Noble et al. proposed a scheme of immersion boundary method in which an additional collision term is added to the Boltzmann equation to reflect the effect of solids on fluid flow [37]. The method has been proved to have sufficient stability in calculating hydrodynamic forces [51]: where Ω i s denotes the additional collision terms for nodes overlapping with solids; B s is the weight coefficient related to ε s : where ε s is percentage of solids in lattice. As Figure 2 shows, when ε s = 0 the lattice is in the flow field; and ε s = 1, the lattice is located at the internal node of the solid. When ε s is between 0 and 1, the lattice is located at the boundary between solid and fluid. The value B s can be obtained by determining the ε s of each node. Therefore, when the particle invades LBM cell, the process of cell transform from fluid node to solid node is gradual. This improves the stability in calculating fluid flow. The drag force and torque acting on the particle induced by fluid flow can be calculated from the following equations: where n denotes the total number of nodes overlapped by particles; (x r − x j ) denotes the relative position of the rth node to the particle j center. Based on the above algorithm, the coupling model based on LBM and DEM is established. Figure 3 shows the process for the LBM-DEM simulation. After the initial conditions of the flow field and the particles are set, at each step of fluid computational motion, executes the DEM cycle. The DEM time step is commonly limited by the explicit schemes due to stability considerations. To improve the DEM calculation accuracy, an extremely short time step should be considered when particles contact, which is much smaller than time scale of fluid motion. Multi-time scale calculation method is used to solve the difference between solid and fluid time scale. We use the fluid time step to

Model Verification
To verify the accuracy of the coupled LBM-DEM model in this paper, the sedimentation process of two particles falling in the rectangular channel which is known as the DKT phenomenon is simulated. In this model, the channel length is 8 cm, width is 2 cm, and the non-slip boundary condition is applied to all walls. The length of the unit grid is 0.01 cm × 0.01 cm. The particle density is 1.01 g/cm 3 and radius is 0.2 cm. The initial position of first particle is (0.99 cm, 7.2 cm) and the position of second particle 2 is (1 cm, 6.8 cm). The fluid kinematic viscosity is 1 × 10 −6 m 2 /s and the density of the fluid is 1 g/cm 3 .
At the beginning, the fluid is stationary and the particles begin to move by gravity. The Figure 4 shows the state of the two deposited particles at each time, showing the collision and rolling process between the particles. Figure 5 shows the change of longitudinal velocity of particles over time and compared with previous working data [52]. Results have a good match with Feng LBM program but there are some differences with Feng IB-LBM program. The difference of simulation results mainly reflects on the time of collision

Model Verification
To verify the accuracy of the coupled LBM-DEM model in this paper, the sedimentation process of two particles falling in the rectangular channel which is known as the DKT phenomenon is simulated. In this model, the channel length is 8 cm, width is 2 cm, and the non-slip boundary condition is applied to all walls. The length of the unit grid is 0.01 cm × 0.01 cm. The particle density is 1.01 g/cm 3 and radius is 0.2 cm. The initial position of first particle is (0.99 cm, 7.2 cm) and the position of second particle 2 is (1 cm, 6.8 cm). The fluid kinematic viscosity is 1 × 10 −6 m 2 /s and the density of the fluid is 1 g/cm 3 .
At the beginning, the fluid is stationary and the particles begin to move by gravity. The Figure 4 shows the state of the two deposited particles at each time, showing the collision and rolling process between the particles. Figure 5 shows the change of longitudinal velocity of particles over time and compared with previous working data [52]. Results have a good match with Feng LBM program but there are some differences with Feng IB-LBM program. The difference of simulation results mainly reflects on the time of collision

Model Verification
To verify the accuracy of the coupled LBM-DEM model in this paper, the sedimentation process of two particles falling in the rectangular channel which is known as the DKT phenomenon is simulated. In this model, the channel length is 8 cm, width is 2 cm, and the non-slip boundary condition is applied to all walls. The length of the unit grid is 0.01 cm × 0.01 cm. The particle density is 1.01 g/cm 3 and radius is 0.2 cm. The initial position of first particle is (0.99 cm, 7.2 cm) and the position of second particle 2 is (1 cm, 6.8 cm). The fluid kinematic viscosity is 1 × 10 −6 m 2 /s and the density of the fluid is 1 g/cm 3 .
At the beginning, the fluid is stationary and the particles begin to move by gravity. The Figure 4 shows the state of the two deposited particles at each time, showing the collision and rolling process between the particles. Figure 5 shows the change of longitudinal velocity of particles over time and compared with previous working data [52]. Results have a good match with Feng LBM program but there are some differences with Feng IB-LBM program. The difference of simulation results mainly reflects on the time of collision between the two particles. In Feng's calculation, the process lasted about 0.9 s, while in this paper only 0.4 s. This difference depends on the different ways of dealing with boundaries, but the state of the particles before and after contact is close. The accuracy of this coupled LBM-DEM method is acceptable compared with the previous experimental results.
between the two particles. In Feng's calculation, the process lasted about 0.9 s, while in this paper only 0.4 s. This difference depends on the different ways of dealing with boundaries, but the state of the particles before and after contact is close. The accuracy of this coupled LBM-DEM method is acceptable compared with the previous experimental results.

Simulation Method and Conditions
Artificial recharge water in GSHP has the characteristics of low nutrition, low ionic strength and large difference of ionic composition [53]. The particle size of micro-solids in recharging water ranges from millimeters to nanometers. Porous media are constructed by randomly generated circle grains, and then imported into the LBM calculation model by image processing technology. These skeleton grains remain fixed during subsequent numerical calculations. As Figure 6a shows, the length of the model is 750 μm and the width is 400 μm. Simulation area is divided into 750 × 400 grids. The left side is a constant velocity boundary, the right side is a constant pressure boundary, and the top and bottom boundary are the periodic boundary. The porous model porosity is 0.45 and the initial permeability is 5.5 μm 2 . Figure 6b shows the diagram of particle diameter statistical distribution for porous media. In order to accelerate the process of particle aggregation and deposition, a smaller particle size ratio was selected. The particle size ratio dp/d50 = 0.06 (here dp is diameter of particle and d50 is the mean diameter of porous particle grain). In this paper, we use pore volume (PV) through the entrance of porous media as a measure this paper only 0.4 s. This difference depends on the different ways of dealing with bo aries, but the state of the particles before and after contact is close. The accuracy of coupled LBM-DEM method is acceptable compared with the previous experimenta sults.

Simulation Method and Conditions
Artificial recharge water in GSHP has the characteristics of low nutrition, low strength and large difference of ionic composition [53]. The particle size of micro-soli recharging water ranges from millimeters to nanometers. Porous media are constru by randomly generated circle grains, and then imported into the LBM calculation m by image processing technology. These skeleton grains remain fixed during subseq numerical calculations. As Figure 6a shows, the length of the model is 750 μm and width is 400 μm. Simulation area is divided into 750 × 400 grids. The left side is a con velocity boundary, the right side is a constant pressure boundary, and the top and bo boundary are the periodic boundary. The porous model porosity is 0.45 and the in permeability is 5.5 μm 2 . Figure 6b shows the diagram of particle diameter statistical tribution for porous media. In order to accelerate the process of particle aggregation deposition, a smaller particle size ratio was selected. The particle size ratio dp/d50 = (here dp is diameter of particle and d50 is the mean diameter of porous particle grain this paper, we use pore volume (PV) through the entrance of porous media as a mea

Simulation Method and Conditions
Artificial recharge water in GSHP has the characteristics of low nutrition, low ionic strength and large difference of ionic composition [53]. The particle size of micro-solids in recharging water ranges from millimeters to nanometers. Porous media are constructed by randomly generated circle grains, and then imported into the LBM calculation model by image processing technology. These skeleton grains remain fixed during subsequent numerical calculations. As Figure 6a shows, the length of the model is 750 µm and the width is 400 µm. Simulation area is divided into 750 × 400 grids. The left side is a constant velocity boundary, the right side is a constant pressure boundary, and the top and bottom boundary are the periodic boundary. The porous model porosity is 0.45 and the initial permeability is 5.5 µm 2 . Figure 6b shows the diagram of particle diameter statistical distribution for porous media. In order to accelerate the process of particle aggregation and deposition, a smaller particle size ratio was selected. The particle size ratio d p /d 50 = 0.06 (here d p is diameter of particle and d 50 is the mean diameter of porous particle grain). In this paper, we use pore volume (PV) through the entrance of porous media as a measure of time, which is defined as the ratio of the volume of water flowing through porous media to the total pore space volume: where V in denotes the volume of injected liquid; V p denotes the volume of pore space. By using PV instead of time, the influence of pore volume and velocity difference on the Processes 2021, 9, 465 9 of 23 experimental results can be eliminated. At the beginning of the simulation, the fluid flows through the porous media until the flow field is stable, and the time is set to the initial time (PV = 0). Particles continually added from the entrance until PV = 3. The basic parameters of the simulation study are summarized in Table 1.

Data Processing
Permeability is a key factor in porous media, which determines the resistance of fluid flowing through it. The retention of particle changes the porosity of the porous media and preferential flow path, thus affecting the permeability. According to Darcy law, permeability in LBM units can be written as: where k denotes the permeability; μ denotes the viscosity of fluid; L denotes the length of the porous media; ∆P denotes the pressure drop. In order to assess the severity of permeability damage, dimensionless remaining permeability is defined as:  Permeability is a key factor in porous media, which determines the resistance of fluid flowing through it. The retention of particle changes the porosity of the porous media and preferential flow path, thus affecting the permeability. According to Darcy law, permeability in LBM units can be written as: where k denotes the permeability; µ denotes the viscosity of fluid; L denotes the length of the porous media; ∆P denotes the pressure drop. In order to assess the severity of permeability damage, dimensionless remaining permeability is defined as: where k r denotes the percentage of remaining permeability; k 0 denotes the initial permeability of porous media. In breakthrough curves, the normalized concentration C r is a function of PV: where C out and C 0 are the volume concentrations in the effluent and influent, separately. To analyze the difference of particle deposition and migration under different conditions, the deposition ratio is determined as percentage of particles retained (M r ).

Grid-Independence Test
Before we proceed further, the effect of the mesh is considered by a large number of numerical simulations, as shown in the Figure 7, permeability is independent of mesh size as long as the lattice spacing is not smaller than 5 × 10 −7 m. Hydrodynamics of fluid applied to particles are calculated by IMBM. When the lattice spacing is less than 10 −6 , the force on the particle is independent of the lattice spacing. Therefore, we will use this lattice spacing in the following simulation.
where kr denotes the percentage of remaining permeability; k0 denotes the initial permeability of porous media. In breakthrough curves, the normalized concentration Cr is a function of PV: = (28) where Cout and C0 are the volume concentrations in the effluent and influent, separately. To analyze the difference of particle deposition and migration under different conditions, the deposition ratio is determined as percentage of particles retained (Mr).

Grid-Independence Test
Before we proceed further, the effect of the mesh is considered by a large number of numerical simulations, as shown in the Figure 7, permeability is independent of mesh size as long as the lattice spacing is not smaller than 5 × 10 −7 m. Hydrodynamics of fluid applied to particles are calculated by IMBM. When the lattice spacing is less than 10 −6 , the force on the particle is independent of the lattice spacing. Therefore, we will use this lattice spacing in the following simulation.

Coupling Effect of Flow Rate and Concentration
This section mainly discusses the retention characteristics of particles and the permeability damage of porous media under different conditions. Particle migration and deposition in porous media of different inlet flow rate (vin = 4, 7, 12, 24 mm/s) and concentration (C0 = 0.005, 0.01, 0.02) was simulated. The suspension is under condition of low ion and neutral suspension. Figure 8 shows the breakthrough curve of effluent concentration and deposition ratio with different concentration condition. When C0 = 0.005, with the increase in inlet flow rate, the plateau values of Cr were 1.1, 0.5, 0.56, and 0.5, respectively; the deposition ratio was 6%, 64%, 61%, and 59%, respectively. In aqueous solution at normal temperature, both silica and kaolin have a negative charge on the surface. This formation creates an electrical double layer and due to low ionic strength, the electric double-layer

Coupling Effect of Flow Rate and Concentration
This section mainly discusses the retention characteristics of particles and the permeability damage of porous media under different conditions. Particle migration and deposition in porous media of different inlet flow rate (v in = 4, 7, 12, 24 mm/s) and concentration (C 0 = 0.005, 0.01, 0.02) was simulated. The suspension is under condition of low ion and neutral suspension. Figure 8 shows the breakthrough curve of effluent concentration and deposition ratio with different concentration condition. When C 0 = 0.005, with the increase in inlet flow rate, the plateau values of C r were 1.1, 0.5, 0.56, and 0.5, respectively; the deposition ratio was 6%, 64%, 61%, and 59%, respectively. In aqueous solution at normal temperature, both silica and kaolin have a negative charge on the surface. This formation creates an electrical double layer and due to low ionic strength, the electric double-layer repulsive force has a long range of action. In low concentration condition (C 0 = 0.005), the case v in = 4 mm/s has the lowest deposition ratio. Due to the high repulsive barrier, the deposition of particle on porous media is limited at low flow rates. When inlet flow rate v in = 7 mm/s, the particles have enough hydrodynamic force to break through the exclusion barrier and deposit on porous media. Beyond the velocity of 7 mm/s, with increasing inlet velocity, the particles tend to be trapped by porous media grain, but disengaged under stronger hydrodynamic shear and the time reaching the peak value of C r increased. Therefore, there is a critical velocity has highest deposition ratio. This point has been found in previous experiments, especially with high repulsive barrier [55,56]. When C 0 = 0.02, with the increase in inlet flow rate, the plateau values of C r were 0.6, 0.32, 0.28, and 0.59, respectively; the retention rate was 53%, 73%, 76%, and 53%, respectively. As Figure 8d shows, for low velocity (v in = 4.0 mm/s), the deposition ratio increased significantly with increasing of concentration. However, for condition v in = 7.5 mm/s and v in = 12 mm/s, C 0 = 0.01 has the lowest deposition rate. For high velocity (v in = 24.0 mm/s), concentration has little effect on deposition ratio. The interference between particles may promote the process both the deposition of particles and the release of deposited particles. Therefore, the effect of concentration on particle migration and deposition in previous studies is not clear [57,58].  Furthermore, the impairment of permeability was observed. Figure 9 shows the transient remaining permeability. The permeability impairment of porous media is similar at low concentration and the permeability for each case is almost unaffected, even if the deposition ratio is really different. When C 0 = 0.01 and 0.02, the permeability impairment of case v in = 4 mm/s is more serious although the deposition rate is lowest. Figure 10 shows a few snapshots of particle retention in the porous media. It can be found that although the deposition ratio is lowest at low flow rate (v in = 4 mm/s), when the particle bridging occurs, it is difficult to destroy bridge due to low shear force. The high shear force at high velocity is more capable to destroy the stability of the bridge. As a result, with the increasing of flow rate, the permeability impairment decreases even with higher deposition ratio.

Coupling Effect of Flow Rate and pH
In groundwater environment, the suspension is slightly acidic or alkaline due to the action of microorganisms. For different pH conditions, the surface electric charge of particle and porous media grains will change, which will affect the interaction between suspensive particle and porous media, and then affect the deposition of particle. The decrease in pH makes the value of the zeta potential of particle and porous media increase and the absolute value decrease, which is beneficial for particle deposited on porous media. In this paper, we assume that particle have the same zeta potential as porous media. When the pH was 5.5, 7, and 9, the value of zeta was −23, −30, and −40 mV, respectively. Figure 11 shows the potential energy curves between particles and porous media grains under different pH condition. When the pH was 5.5, 7, and 9, the repulsive barrier is 1400 KT, 2800 KT, and 5300 KT, respectively. The secondary potential well is lower than 1 KT and can be neglected under hydrodynamic action. Figure 12a-c shows the breakthrough curve of effluent concentration with different pH condition (pH = 5. 5,7,9). When the suspension is slightly acidic (pH = 5.5), and the repulsive barrier is slow, as the flow rate increased from 4 to 12 mm/s, the plateau values of Cr was 0.35, 0.4, 0.47, and 0.67, respectively; the retention rate was 69%, 66%, 60%, and 55%, respectively. When the suspension is slightly alkaline (pH = 9), as the flow rate increased from 4 to 12 mm/s, the plateau values of Cr was 1.04, 0.86, 0.45, and 0.46, respectively and the deposition rate was 2%, 34%, 61%, and 65%. In the slightly acidic suspension with lower repulsive barrier, the retention rate decreases with the increase in flow rate. The phoneme in alkaline suspension with high barrier is totally opposite. Under acidic conditions, the interaction the double layer electrical force is weak, and the particles can be easily deposited on the surface of porous media. The increasing of flow rate, result in the increases of shear force, which makes the particles more easily detached. However, under alkaline conditions, the lower flow rate makes it difficult for particles deposits on the surface of porous media. The increase in the barrier increases the critical velocity. However, it is worth noting that DLVO theory describes surface interactions qualitatively rather than quantitatively. The repulsive barrier calculated by DLVO theory is often higher than the real condition, especially in the unfavorable condition [56]. As a result, the

Coupling Effect of Flow Rate and pH
In groundwater environment, the suspension is slightly acidic or alkaline due to the action of microorganisms. For different pH conditions, the surface electric charge of particle and porous media grains will change, which will affect the interaction between suspensive particle and porous media, and then affect the deposition of particle. The decrease in pH makes the value of the zeta potential of particle and porous media increase and the absolute value decrease, which is beneficial for particle deposited on porous media. In this paper, we assume that particle have the same zeta potential as porous media. When the pH was 5.5, 7, and 9, the value of zeta was −23, −30, and −40 mV, respectively. Figure 11 shows the potential energy curves between particles and porous media grains under different pH condition. When the pH was 5.5, 7, and 9, the repulsive barrier is 1400 KT, 2800 KT, and 5300 KT, respectively. The secondary potential well is lower than 1 KT and can be neglected under hydrodynamic action. Figure 12a-c shows the breakthrough curve of effluent concentration with different pH condition (pH = 5. 5,7,9). When the suspension is slightly acidic (pH = 5.5), and the repulsive barrier is slow, as the flow rate increased from 4 to 12 mm/s, the plateau values of C r was 0.35, 0.4, 0.47, and 0.67, respectively; the retention rate was 69%, 66%, 60%, and 55%, respectively. When the suspension is slightly alkaline (pH = 9), as the flow rate increased from 4 to 12 mm/s, the plateau values of C r was 1.04, 0.86, 0.45, and 0.46, respectively and the deposition rate was 2%, 34%, 61%, and 65%. In the slightly acidic suspension with lower repulsive barrier, the retention rate decreases with the increase in flow rate. The phoneme in alkaline suspension with high barrier is totally opposite. Under acidic conditions, the interaction the double layer electrical force is weak, and the particles can be easily deposited on the surface of porous media. The increasing of flow rate, result in the increases of shear force, which makes the particles more easily detached. However, under alkaline conditions, the lower flow rate makes it difficult for particles deposits on the surface of porous media. The increase in the barrier increases the critical velocity. However, it is worth noting that DLVO theory describes surface interactions qualitatively rather than quantitatively. The repulsive barrier calculated by DLVO theory is often higher than the real condition, especially in the unfavorable condition [56]. As a result, the critical velocity in this paper is probably higher than previous experiment study.
Processes 2021, 9, x FOR PEER REVIEW 14 of 23 Figure 11. The potential energy curves between particles and porous media grain.    As Figure 12d shows, in the condition of low flow rate, the decrease in the repulsive barrier makes the deposition ratio of the particle increase obviously. The lower barrier can obviously help the deposition of particles. However, when the flow rate is high, the change of pH has little effect on deposition ratio. Figure 13 shows the transient remaining permeability. We note that even there are great differences in deposition ratio under different conditions, the change of permeability impairment is small. The decrease in permeability is mainly due to the retention bridge of particles at the pore throat. We speculate that the production of bridging is mainly affected by concentration.

Particle Micro Deposition Behavior
The advantage of simulation is the ability to trace all particles at pore scale. We wer able to investigate the different microscopic events leading to particle deposition. Whe concentration is low, particle deposition is only related to the interaction between particl and porous grain. At different incident velocities, the force, velocity and trajectory of par ticle is different during collision. As Figure 14 shows, when the particle approach grai the DLVO potential energy diagram sequentially appeared as a secondary potential wel a repulsive barrier and a primary potential well. The simulated environment in this pape is also under the condition of a low ionic concentration environment, and the electrostati shielding effect is weak. The secondary potential well can be neglected under hydrody namic action. In order to reach the primary potential well, the drag force needs to over come the repulsive barrier.
The deposition of particle is mainly related to inertial deposition and interception The flow field is the most direct factors affecting particle deposition. The particle flow along the streamline and having potential contact with the grains of porous media. Figur 15 shows the collision between particle and grain at different flow rate. When vin = 4 mm/ particle cannot overcome the repulsive barrier at low flow rate and hard to deposit on th porous media. For vin = 7.5 mm/s, particle has enough kinetic energy overcome the repu sive barrier and enter the primary potential well. The particles accelerate under the actio of adhesive action after breaking through the repulsive barrier and contact with grain f nally. The particles do not deposit directly in porous media after a collision. After collidin between particle and grain, the particle rebounds in the opposite direction. Due to th energy of the damping dissipates the kinetic energy of the particle, the rebound velocit is smaller than the velocity at the beginning of the collision. However, the particles canno reach the separation point at the end of the first rebound, and are pulled back for a new approach and subsequent rebound. The kinetic energy is gradually dissipated in the pro cess, and the particles are finally deposited on the porous media.

Particle Micro Deposition Behavior
The advantage of simulation is the ability to trace all particles at pore scale. We were able to investigate the different microscopic events leading to particle deposition. When concentration is low, particle deposition is only related to the interaction between particle and porous grain. At different incident velocities, the force, velocity and trajectory of particle is different during collision. As Figure 14 shows, when the particle approach grain the DLVO potential energy diagram sequentially appeared as a secondary potential well, a repulsive barrier and a primary potential well. The simulated environment in this paper is also under the condition of a low ionic concentration environment, and the electrostatic shielding effect is weak. The secondary potential well can be neglected under hydrodynamic action. In order to reach the primary potential well, the drag force needs to overcome the repulsive barrier.
Processes 2021, 9, x FOR PEER REVIEW 16 of 23 throat, the following particle in the center of the channel collide and flocculate with the particle at the edge of the channel. Particles flocculated into large aggregates and occupy the throat, which leads to the permeability impairment. According to above section, under the condition vin = 4 mm/s, as the concentration increased from 0.05 to 0.2, the deposition rate increased from 6% to 53%. Therefore, it can be seen that the deposition rate of particle is low at low concentration and low flow rate, but the deposition of particle at high concentration is very obvious.  The deposition of particle is mainly related to inertial deposition and interception. The flow field is the most direct factors affecting particle deposition. The particle flow along the streamline and having potential contact with the grains of porous media. Figure 15 shows the collision between particle and grain at different flow rate. When v in = 4 mm/s, particle cannot overcome the repulsive barrier at low flow rate and hard to deposit on the porous media. For v in = 7.5 mm/s, particle has enough kinetic energy overcome the repulsive barrier and enter the primary potential well. The particles accelerate under the action of adhesive action after breaking through the repulsive barrier and contact with grain finally. The particles do not deposit directly in porous media after a collision. After colliding between particle and grain, the particle rebounds in the opposite direction. Due to the energy of the damping dissipates the kinetic energy of the particle, the rebound velocity is smaller than the velocity at the beginning of the collision. However, the particles cannot reach the separation point at the end of the first rebound, and are pulled back for a new approach and subsequent rebound. The kinetic energy is gradually dissipated in the process, and the particles are finally deposited on the porous media.
the throat, which leads to the permeability impairment. According to above section, under the condition vin = 4 mm/s, as the concentration increased from 0.05 to 0.2, the deposition rate increased from 6% to 53%. Therefore, it can be seen that the deposition rate of particle is low at low concentration and low flow rate, but the deposition of particle at high concentration is very obvious.   However, in high concentration suspension, the deposition of particle is usually affected by the action of other particles. Figures 16 and 17 show the transportation of some particles in the throat under low velocity (v in = 4 mm/s) and high concentration (C 0 = 0.02) condition. When multiple particles pass through the pore throat, as Figure 16 shows, blue particle near the grain surface transport slower through the pore throat, and the red particle at the center of the throat has higher speeds. The blue particle could not break through the repulsive barrier through hydrodynamic and inertial effects. With the repulsive interference of red particle and the impact of hydrodynamic forces, the blue particle deposit on the surface of grain. In addition, the increase in concentration promotes the flocculation between particles. As Figure 17 shows, due to the velocity gradient in the throat, the following particle in the center of the channel collide and flocculate with the particle at the edge of the channel. Particles flocculated into large aggregates and occupy the throat, which leads to the permeability impairment. According to above section, under the condition v in = 4 mm/s, as the concentration increased from 0.05 to 0.2, the deposition rate increased from 6% to 53%. Therefore, it can be seen that the deposition rate of particle is low at low concentration and low flow rate, but the deposition of particle at high concentration is very obvious.

Dynamics of the Clogging Process and Resuspension
Compared with the previous section, the dynamic behavior after particle deposition is shown in this section. As demonstrated in Figure 18, the contact between particle and porous grain is often upstream, but deposition occurs downstream. During the clogging process, due to hydrodynamic and shear torque, the particles roll with the surface and tend to deposit in the downstream zone of the porous grain. This mainly lies in the assumption that the surface of porous media is smooth and ignoring the roughness of grain surface [59]. Similar particle behavior was also found in some experiments. Particles deposit on the surface or at the secondary potential well and migrate to the downstream as water flows [60]. The movement of particles downstream may be the decisive factor in clogging. Deposition of particles expands the surface area of the collector, making subsequent particles more accessible to contract with grain. As Figure 19 shows, with the accumulation of downstream deposition, the flow throat is gradually encroached until the particles form a bridge, which hinder the fluid flow under the action of the particle bonds. Therefore, the internal structure of porous media plays an important role in the formation of clogging layer. This is similar to the results of previous studies [61]. In their study, the critical concentration of particles was 20%. When the concentration is lower than this value, the clogging layer occurs through continuous deposition, whereas for 20% concentration, sudden bridging formation can occur. However, the particle concentration does

Dynamics of the Clogging Process and Resuspension
Compared with the previous section, the dynamic behavior after particle depositi is shown in this section. As demonstrated in Figure 18, the contact between particle a porous grain is often upstream, but deposition occurs downstream. During the cloggi process, due to hydrodynamic and shear torque, the particles roll with the surface a tend to deposit in the downstream zone of the porous grain. This mainly lies in the sumption that the surface of porous media is smooth and ignoring the roughness of gra surface [59]. Similar particle behavior was also found in some experiments. Particles d posit on the surface or at the secondary potential well and migrate to the downstream water flows [60]. The movement of particles downstream may be the decisive factor clogging. Deposition of particles expands the surface area of the collector, making sub quent particles more accessible to contract with grain. As Figure 19 shows, with the acc mulation of downstream deposition, the flow throat is gradually encroached until the p ticles form a bridge, which hinder the fluid flow under the action of the particle bond

Dynamics of the Clogging Process and Resuspension
Compared with the previous section, the dynamic behavior after particle deposition is shown in this section. As demonstrated in Figure 18, the contact between particle and porous grain is often upstream, but deposition occurs downstream. During the clogging process, due to hydrodynamic and shear torque, the particles roll with the surface and tend to deposit in the downstream zone of the porous grain. This mainly lies in the assumption that the surface of porous media is smooth and ignoring the roughness of grain surface [59]. Similar particle behavior was also found in some experiments. Particles deposit on the surface or at the secondary potential well and migrate to the downstream as water flows [60]. The movement of particles downstream may be the decisive factor in clogging. Deposition of particles expands the surface area of the collector, making subsequent particles more accessible to contract with grain. As Figure 19 shows, with the accumulation of downstream deposition, the flow throat is gradually encroached until the particles form a bridge, which hinder the fluid flow under the action of the particle bonds. Therefore, the internal structure of porous media plays an important role in the formation of clogging layer. This is similar to the results of previous studies [61]. In their study, the critical concentration of particles was 20%. When the concentration is lower than this value, the clogging layer occurs through continuous deposition, whereas for 20% concentration, sudden bridging formation can occur. However, the particle concentration does not reach this value during recharge process. Therefore, the formation of clogging layer should be produced by continuous deposition during recharge process.
Furthermore, in the case of high flow rate, the fluid has strong shear force on deposited particles. Hydrodynamic effects can significantly alter particle deposition, especially detach deposited particles from the primary minimum well [62,63]. As the flow rate increases, more particle will leave the surface and reduce the deposition rate in favorable deposition environment. As Figure 21 shows, after interaction of particle through the throat, the deposited particle is separated from the surface and resuspended. Due to the repulsive force of the deposited particles, suspend particle accelerate or slow in the throat, thus changing the flow field distribution of the throat. The change of flow field increases the lift force of the deposited particles and promotes the resuspend of the deposited particles. Therefore, the interaction between particles may promote the deposition process or facilitate the release of deposited particle. This explains why the deposition rate increases first and then decreases with the increase in particle concentration when vin = 7.12 mm/s. It is important to note that, the effect of increasing hydrodynamic force on the release of sediment particles is clear but not obvious compared to previous research. It mainly due to not considering the complexity of natural environment, such as roughness, which is expected to play a major role both in particle retention and release [64,65].   Furthermore, in the case of high flow rate, the fluid has strong shear force on deposited particles. Hydrodynamic effects can significantly alter particle deposition, especially detach deposited particles from the primary minimum well [62,63]. As the flow rate increases, more particle will leave the surface and reduce the deposition rate in favorable deposition environment. As Figure 21 shows, after interaction of particle through the throat, the deposited particle is separated from the surface and resuspended. Due to the repulsive force of the deposited particles, suspend particle accelerate or slow in the throat, thus changing the flow field distribution of the throat. The change of flow field increases the lift force of the deposited particles and promotes the resuspend of the deposited particles. Therefore, the interaction between particles may promote the deposition process or facilitate the release of deposited particle. This explains why the deposition rate increases first and then decreases with the increase in particle concentration when vin = 7.12 mm/s. It is important to note that, the effect of increasing hydrodynamic force on the release of sediment particles is clear but not obvious compared to previous research. It mainly due to not considering the complexity of natural environment, such as roughness, which is expected to play a major role both in particle retention and release [64,65].   Resuspension is often observed during clogging process, especially in condition of high flow rate. The interaction between particles is much smaller than that between particle and grain which leads to the breakage of particle bonds. As Figure 20 shows, particle bonds break and the throat restored permeability. After resuspension, the agglomerate is deposited downstream or enters deeper into the porous media. Due to large volume, there-suspended agglomerates tend to deposit again. Breakage of particle aggregates usually releases some individual particles, as the yellow particles shown in Figure 20. The stability of the particle bridge is directly related to the magnitude of the hydrodynamic force. This explains that at high flow rates condition, the low permeability decay even with higher deposition rate. Processes 2021, 9, x FOR PEER REVIEW 19 of 23 Figure 20. The breakage of particle bonds and release particle: particle distribution and flow field vector diagram (a) Clogging state; (b) Breakage process of particle bonds; (c) Release particle after breakage of bonds.

Summary and Conclusions
Understanding the migration and deposition mechanism of particles in porous media and predicting the clogging formation has important academic meanings and engineering significance.
Taking IMBM mode as the connection, we coupled the LBM-DEM method to simulate particle migration and retention inside the porous media during recharging groundwater. This method identifies important microscopic events. Based on DLVO theory, the deposition of particles and the permeability attenuation mechanism of porous media are studied. The coupled effects of concentration, flow rate and pH on the clogging mechanism of porous media were analyzed. For single particles, due to the repulsive barrier between particles and porous media, there is a critical velocity. At low velocity, the deposition ratio increases with the increase in velocity. Beyond the critical velocity, the deposition ratio decreases with the velocity increases due to higher shear force. The increase in concentration makes the interaction between particles more frequent. The interaction of particles promotes deposition of particles in low flow rate and also promotes the release of deposited particles in high flow rate. In general, the increase in concentration increases permeability impairment at all velocity conditions, and the increase in concentration in the low flow rate condition is more obvious for permeability impairment. When the concentration is constant, the increase in flow rate improves the permeability impairment, even with a higher deposition ratio. Changes in pH mainly affect the value of the repulsive barrier. Under the acidic condition, the increase in flow rate reduces the deposition rate. Under the alkaline condition, the low flow rate condition has a low deposition ratio due to it being hard for particles to break through the repulsive barrier. For low flow rate, the Figure 20. The breakage of particle bonds and release particle: particle distribution and flow field vector diagram (a) Clogging state; (b) Breakage process of particle bonds; (c) Release particle after breakage of bonds.
Furthermore, in the case of high flow rate, the fluid has strong shear force on deposited particles. Hydrodynamic effects can significantly alter particle deposition, especially detach deposited particles from the primary minimum well [62,63]. As the flow rate increases, more particle will leave the surface and reduce the deposition rate in favorable deposition environment. As Figure 21 shows, after interaction of particle through the throat, the deposited particle is separated from the surface and resuspended. Due to the repulsive force of the deposited particles, suspend particle accelerate or slow in the throat, thus changing the flow field distribution of the throat. The change of flow field increases the lift force of the deposited particles and promotes the resuspend of the deposited particles. Therefore, the interaction between particles may promote the deposition process or facilitate the release of deposited particle. This explains why the deposition rate increases first and then decreases with the increase in particle concentration when v in = 7.12 mm/s. It is important to note that, the effect of increasing hydrodynamic force on the release of sediment particles is clear but not obvious compared to previous research. It mainly due to not considering the complexity of natural environment, such as roughness, which is expected to play a major role both in particle retention and release [64,65].

Summary and Conclusions
Understanding the migration and deposition mechanism of particles in porous media and predicting the clogging formation has important academic meanings and engineering significance.
Taking IMBM mode as the connection, we coupled the LBM-DEM method to simulate particle migration and retention inside the porous media during recharging groundwater. This method identifies important microscopic events. Based on DLVO theory, the deposition of particles and the permeability attenuation mechanism of porous media are studied. The coupled effects of concentration, flow rate and pH on the clogging mechanism of porous media were analyzed. For single particles, due to the repulsive barrier between particles and porous media, there is a critical velocity. At low velocity, the deposition ratio increases with the increase in velocity. Beyond the critical velocity, the depo-

Summary and Conclusions
Understanding the migration and deposition mechanism of particles in porous media and predicting the clogging formation has important academic meanings and engineering significance.
Taking IMBM mode as the connection, we coupled the LBM-DEM method to simulate particle migration and retention inside the porous media during recharging groundwater. This method identifies important microscopic events. Based on DLVO theory, the deposi-tion of particles and the permeability attenuation mechanism of porous media are studied. The coupled effects of concentration, flow rate and pH on the clogging mechanism of porous media were analyzed. For single particles, due to the repulsive barrier between particles and porous media, there is a critical velocity. At low velocity, the deposition ratio increases with the increase in velocity. Beyond the critical velocity, the deposition ratio decreases with the velocity increases due to higher shear force. The increase in concentration makes the interaction between particles more frequent. The interaction of particles promotes deposition of particles in low flow rate and also promotes the release of deposited particles in high flow rate. In general, the increase in concentration increases permeability impairment at all velocity conditions, and the increase in concentration in the low flow rate condition is more obvious for permeability impairment. When the concentration is constant, the increase in flow rate improves the permeability impairment, even with a higher deposition ratio. Changes in pH mainly affect the value of the repulsive barrier. Under the acidic condition, the increase in flow rate reduces the deposition rate. Under the alkaline condition, the low flow rate condition has a low deposition ratio due to it being hard for particles to break through the repulsive barrier. For low flow rate, the decrease in repulsive barrier in the acidic condition greatly promotes the deposition of particles. For high flow rate, particles are more affected by hydrodynamic forces, which make the pH less affected in particle deposition.
In addition, the deposition phenomenon in porous media with micron size particles is more fully understood by simulation. The process of particle deposition and the dynamic motion after deposition were observed such as particles gliding over the surface. Accumulated particles in the downstream form bridges and hinder fluid flow. In high flow rate, strong shear force is more capable of destroying bridges and recovering permeability. Adsorbed particles glide on the surface of the grain and deposit in the downstream. This is almost undetectable in an experimental study as it would require a nanometer spatial resolution and a very high temporal resolution. These micro-phenomena help explain the macro-characteristics of particle migration and deposition.
This method has great potential in studying the migration of particles in porous media. However, there are still some improvements in improving the accuracy and reducing the cost of calculation. Furthermore, groundwater recharging is a very complex microscopic process and the deposition mechanics of particles with different physical properties in porous media vary greatly, which needs further study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to restrictions of privacy.

Nomenclature x
Position of particle/Grid coordinate i, j Particle index/Direction in LBM t, ∆t Time and Time step v p, ω p Velocity and angular velocity of particle m, I Mass of particle and moment of inertia F, T Force and torque on particles F f , T f Hydrodynamic force and torque induced by fluid flow F c , T c Contact forces and Torque resulting from the contact process F ad , F re Adhesive force and repulsive force F g Gravity force n, t