Particle clogging owing to the confinement of the pore structure affects the transport and retention of particles. In this study, fluid–particle flow in a porous medium was modeled using computational fluid dynamics–discrete element method (CFD-DEM) simulations. The motion of particles in the slurry was traced in the porous medium, which enabled particle clogging to be directly precited by the interaction between the particles and pore surfaces by assessing the exact location of each particle. The pressure and flow field of the liquid were investigated, and the variation in flow path owing to particle clogging was predicted. The hydrodynamic study also showed that the Stokes number and particle concentration determined the particle clogging at the pore throats of the porous medium. Keywords : CFD-DEM, particle retention, particle transport, path conversion, porous media 1 Introduction Understanding the transport and retention of particles in porous media is of great significance for various industrial applications, such as particulate filtration in diesel engines in the automotive industry [1,2], water or flue gas treatment using packed granular bed filtration in the processing industry [3,4], and nanoparticle injection for enhanced oil recovery in the petroleum industry [5,6,7]. Particle migration in porous media is constrained by the pore structure. In addition, the structure and topology of the porous medium are altered by particles that become trapped. Clogging owing to the accumulation of particles can reduce the permeability of the porous medium, which can harm productivity [8]. Therefore, the effect of particle retention on porous media is of great concern [9,10,11]. Recent studies investigated the microstructural variation of porous media affected by the filtration process, in addition to the effect of particle retention on the permeability and porosity of porous media studied at the macroscopic level. Chen et al. [12] proposed a “gathered-energized” process of the particle transport in a porous medium from their microfluidic experiments. Three areas were found in the microfluidic structure of the porous media 1) unrestricted flow area, 2) filter cake layers, and 3) particle deep migration area downstream in the porous media. In the filter cake layer, some particles were retained while others entered the deep core of the third layer. Notably, particle clogging at the pore throats exponentially increased the pressure drop in the system and caused pressure fluctuations [12]. This hydrodynamic instability within the porous medium caused the system failure in the process operation. Kokubun et al. [13] investigated the hydrodynamic effects of particle accumulation in porous media. The study showed that clogging causes pressure redistribution and diverts upstream water into neighboring pores. In addition, they noted that clogging occurs in porous media with a heterogeneous pore size distribution, but not in microscopically homogeneous media. Fo et al. [14] studied the fine particle migration and blockage in porous media and their effects on the porosity and fluid velocity. They performed lattice Boltzmann method (LBM)–immersed boundary method (IBM)–discrete 1 Corresponding author. E-mail address: dansuncfd@outlook.com element method (DEM) simulations, wherein the particles were simulated using the Lagrangian method and DEM to trace the motion of each particle, and the fluid was simulated using LBM. Zhou et al. [15] also used LBM-DEM-IBM to simulate particle retention and the associated permeability impairment in porous media owing to size exclusion effects. In their research, particles that were larger than the pore throat blocked the flow path and reduced the permeability of the porous medium. Particle simulation using the Lagrangian method enables the positions of particles to be predicted and traced accurately, and the simulation results provide information on the collection of particles based on the pore structure. This provides evidence for studying the mechanism of particle clogging in porous media. In this study, the migration and retention of particles in a slurry were simulated as they passed through a saturated porous medium using a Eulerian–Lagrangian simulation method, computational fluid dynamics (CFD)-DEM. The fluid in the slurry model was water, and the Knudsen number was sufficiently low for water in the porous medium flow in the continuous region. Therefore, the Navier– Stokes equation was used to model the fluid phase in the slurry. The particle filtration in porous media considered the mechanisms such as size exclusion, bridging, segregation, straining, diffusion, electrostatic attraction [16]. In this study, we first investigated the effects of structural confinement within a porous medium on the migration and retention of particles. In particular, the simulation was designed to simulate the mechanism of particle bridging. A homogeneous porous medium with a uniform pore size was generated to study the effect of particle size on particle clogging. The hydrodynamics of the liquid passing through the porous medium before and after particle clogging were simulated to study how clogging affects the filter cake permeability, as well as the fate and conversion of the flow path through the homogeneous porous medium. 2 CFD-DEM method Particle filtration in a porous medium was simulated using the Eulerian–Lagrangian simulation method, CFD-DEM. In this method, the fluid dynamics of the liquid were simulated by solving the Reynolds- averaged Navier–Stokes equation [17]: డఈ೗ఘ೗ డ௧+∇∙(𝛼௟𝜌௟𝑢௟ሬሬሬ⃑)= 0 (1) డఈ೗ఘ೗௨೗ሬሬሬሬ⃑ డ௧+∇∙(𝛼௟𝜌௟𝑢௟ሬሬሬ⃑𝑢௟ሬሬሬ⃑)= −∇p+∇∙ (𝜏௟ሬሬሬ⃑)−ଵ ఈ೗௏ౙ౛ౢౢ∙ห∑൫𝛼௟𝐹஽,పሬሬሬሬሬሬ⃑൯ ௜ห∙௨೗ሬሬሬሬ⃑ି௨ೞሬሬሬሬ⃑ |௨೗ሬሬሬሬ⃑ି௨ೞሬሬሬሬ⃑| (2) where αl is the volume fraction of the liquid, subscript l represents the liquid phase, ρ is density, u is mean velocity, and 𝜏௟ሬሬሬ⃑ is the stress tensor of liquid flow ( 𝜏௟ሬሬሬ⃑=(𝜆௟−(2 3⁄)𝜇௟)(∇∙𝑢௟ሬሬሬ⃗)𝐼⃗+𝜇௟൫(∇𝑢௟ሬሬሬ⃗)+ (∇𝑢௟ሬሬሬ⃗)்൯). DEM was used to trace the particles under the Lagrangian reference frame. Discrete particle trajectories were predicted by integrating the force balance of the particles. The governing equations of the translational and rotational motions of the particles are as follows: 𝑚ௗ௨೛ሬሬሬሬሬ⃗ ௗ௧=𝐹஽ሬሬሬሬ⃗+𝐹௖ሬሬሬ⃗ (3) 𝐼ௗఠ೛ሬሬሬሬሬሬ⃗ ௗ௧=𝑇௖ሬሬሬ⃗+𝑇௥ሬሬሬ⃗ (4) where 𝐹஽ሬሬሬሬ⃗ is the body force on each particle from the interaction with the fluid phase and 𝐹௖ሬሬሬ⃗ is the contact force due to particle–particle collisions. 2.1 Drag force The drag force on a single particle was calculated as [18,19] 𝐹஽ሬሬሬሬ⃗=ଵ ଼𝐶஽𝜌௟𝜋𝑑௣ଶห𝑢௟ሬሬሬ⃗−𝑢௣ሬሬሬሬ⃗ห൫𝑢௟ሬሬሬ⃗−𝑢௣ሬሬሬሬ⃗൯ (5) where the drag coefficient CD is a function of the particle Reynolds number [20]: 𝐶஽=൬0.62+ସ.଼ ඥୖୣ౦൰ଶ 𝛼௚ଶିఞ (6) 𝜒= 3.7−0.65 exp ቂ−൫1.5−log ଵ଴Re௣൯ଶ/2ቃ (7) and Re p is the Reynolds number of the particle, calculated by Re௣=ఓୀఘ೒ఌ೒ห𝐮𝐠ି𝐮𝐬ห ఓ೒ (8) 2.2 Contact force In DEM, contact forces are simulated using a spring–dashpot model [21]. The spring generates repulsive forces pushing particles apart, and the dashpot represents viscous damping and allows the simulation of collision types other than perfectly elastic. A parallel linear spring–dashpot model represents the normal force, and a parallel linear spring–dashpot model in series with a slider represents the tangential force with respect to the normal vector of the contact plane. In both cases, the spring accounts for the elastic part of the response, and the dashpot accounts for the energy dissipation during the collision. A linear spring contact model was used here to define the contact force [22,23]. The normal component of the contact force was expressed as 𝐹௡= −𝑘௡𝑑௡−𝑐௡൫𝑢௡,௜−𝑢௡,௝൯ (9) Tangential shear forces occur upon collision and increase the torque. The Coulomb friction law was used to determine the tangential contact force: 𝐹௧= minቄቚ𝑘௧∫൫𝑢௧,௜−𝑢௧,௝൯𝑑𝑡௧ ௧೎,଴+𝑐௧൫𝑢௧,௜−𝑢௧,௝൯ቚ,𝜇௖𝐹௡ቅ (10) where 𝑑௡ and 𝑑௧ are the overlaps in the normal and tangential directions, respectively, at the contact point; μc is the Coulomb friction coefficient; k is the spring stiffness; and c is the damping coefficient. The torque due to friction at the contact point was calculated by 𝑇௖ሬሬሬ⃗= r௖×𝐹௧ (10) In addition, a directional contact torque model was applied as a rolling friction model [24]: 𝑇௥ሬሬሬ⃗=𝑅ఓ𝑘௡∆𝑑௡ఠೝ೐೗ |ఠೝ೐೗|𝑟௣ (11) where Rµ is the rolling friction model parameter and the relative angular velocity at contact is defined as 𝜔௥௘௟ሬሬሬሬሬሬሬሬ⃗=௥೔ఠഢሬሬሬሬሬ⃗ା௥ೕఠണሬሬሬሬሬ⃗ ௥೔ା௥ೕ (12) The simulation was carried out using the open-source software CFDEM, and the detailed physical model and numerical method can be found in the publication by Goniva et al. [23]. 3 Simulations and analysis In the simulation, the porous medium was modeled as homogeneously packed spherical particles, as shown in Figure 1(a). The diameter of the packed particles was 1 mm, and the distance between two neighboring particles was 1.05 mm, as measured from the centers of the spheres. Thus, the narrowest gap in the pore throat was 50 μm. The aim of this study was to investigate the particle bridging mechanism of particle retention in porous media; therefore, the injected particles should be smaller than the minimum pore-throat gap (<50 μm diameter). This is less than 1/20 of the diameter of the packed particles of the porous medium. Consequently, the three-dimensional porous medium was simplified into a two-dimensional geometry in the model, because the flowing particles are trapped at the cross-sections between the centers of neighboring spherical packed particles, where the throat gap is the narrowest. The simulation was performed in three dimensions with a periodic boundary condition in the perpendicular direction, and the thickness of the two-dimensional structure shown in Figure 1(a) was assumed to be 50 μm, which is equivalent to the minimum pore-throat gap. In the simulation, the density of the liquid was 1.0 × 103 kg/m3, and the viscosity was 1.23 × 10−3 kg/m‧s. The density of the particles was 2.4 × 103 kg/m3, and the particle size distribution is shown in Figure 1(b). Figure 1. Diagram of porous medium and particle size distribution in the simulation Figure 2. Fluid dynamics of water without particles injected in the porous medium Figure 3. Velocity distribution of fluid at the pore throat 3.1 Fluid dynamics in the porous medium Figure 2 depicts the dynamics of fluid flow through the porous medium. The fluid flow was pressure- driven with a constant pressure drop of 1.0 × 104 Pa. In the absence of injected particles, the pressure distribution was uniform, with a gradual reduction of pressure from the inlet on the left to the outlet on the right, as shown in Figure 2(a). Figure 2(b) shows the velocity of the liquid in the porous medium. The liquid velocity was highest at the pore throats, reaching a maximum of 1.2 m/s. Figure 3 depicts the velocity of the liquid at the cross-sections of the narrowest gap of the pore throat, as labeled in Figure 2(b). 3.2 Particle transport in the porous medium When polydisperse particles are injected with water into the porous medium, the particles travel along the flow path of the liquid. As shown in Figure 4, when the mean diameter of the particles was 15 μm (i.e., when the particle size was relatively small), no particles were trapped in the porous medium, and all particles passed through the porous medium as the time progressed. Figure 4(a), (b), and (c) respectively show the particles as they entered the porous medium, as they approached the middle of the porous medium, and when the particle flow was fully developed. Figure 4. Particle transport in the porous medium when the mean particle diameter was 15 μm By dividing the porous medium into different sections in the lateral direction, as shown in Figure 1(a), the properties of the particles in each section were statistically analyzed. Figure 5(a) and (b) show the particle size distributions in each section. There was little difference in the mass fraction of the smallest particles (5-µm diameter) in each section; however, the mass fraction of larger particles (25- µm diameter) was lower at the entrance of the porous medium (section A) than near the outlet (sections D and E). This is because larger particles have a lower velocity, which means they have a greater resident time than smaller particles. Figure 5 shows the lateral velocity of the particles normalized by the maximum velocity of the liquid plotted in Figure 3. The lateral velocity was averaged for each section. There was little difference in the velocities of the smaller particles (<15 μm) in each section, whereas larger particles (35 μm) had larger differences in velocity. The peak velocities in each section corresponded to different particle sizes; for example, in section A, the 15 μm particles had the highest velocity, whereas in sections B and E, the 20 µm particles had the highest velocity. The Stokes number (St) was calculated for each size of particles in the porous medium. As shown in Figure 3, when St equaled 1, the particle diameter was between 15 and 20 µm. When St was less than 1, the particles were mainly manipulated by the liquid flow, and there were no significant differences for the different sections for particles smaller than 15 μm. When St was greater than 1, the inertia of the particles dominated the particle movement, and some particles’ positions were affected by the flow path of the liquid and the surfaces of the porous medium. For example, the larger particles that collided with the pore walls bounced backwards, which reduced the downstream lateral velocity. Figure 5. Particle distribution and velocity along the lateral direction of the porous medium Figure 6. Retention of particles with mean diameters of a) 25 μm and b) 35 μm in the porous medium 3.3 Particle retention in the porous medium Particle clogging occurs when the diameter of the particles injected into the porous media is increased. Figure 6 shows the predicted particle distributions in the porous medium when the mean particle diameter was increased to 25 and 35 µm. The particles were trapped at the entrances of the pores via bridging. As shown in Figure 6(b), more pore throats were blocked by the particles when the mean particle diameter was larger. Figure 7 shows the particle distribution over time when the mean particle diameter was 35 μm. At t1, particle clogging occurred near the pore entrance at a random location. As the time progressed, the pore throat behind it became blocked. Meanwhile, the separation efficiency of the particles in the porous medium increased, as shown in Figure 8. When particle clogging occurs, the structure of the porous media is altered. Figure 9 shows the variation of pressure distribution with time as more particles block the pore spaces, compared with the uniformly distributed pressure drop in Figure 2(a). As shown in Figure 9, the filtered particles in the porous medium affect the pressure distribution. Because the fluid moves by pressure-driven flow when there is a consistent pressure drop in the system, when the first particles become trapped near the entrance of the pore throat (i.e., at time t1), the local pressure drop increases. Figure 7. Retention of particles with a mean diameter of 35 μm at a) t1, b) t1 + 0.02 s, and c) t1 + 0.04 s Figure 8. Separation efficiency over time The pressure also sharply dropped at each clogged region, typically in section A, where more particles were trapped at the bottom of the porous medium. As shown in Figure 9 Figure 10 , the local pressure variation in even 0.02 s intervals was significant; therefore, the pressure varies frequently as particles become trapped during the operation. 3.4 Flow path conversion due to particle retention The velocity of particles in porous media varies during the process because particle clogging causes porosity variation and pressure redistribution and fluctuation. Figure 10 shows the predicted flow stream variation over time. The distribution of the mean velocity of the fluid in the porous medium varies from Figure 10(a) to Figure 10(c). This is because the permeability is correlated to the porosity in the porous medium, which is affected by particle clogging. The velocity distribution in Figure 10(a) is at the same as the distribution in Figure 7(a) and Figure 9(a), when some particles blocked the pore throats near the pore entrance. The particle velocity was reduced downstream the two clogged throats, and the mass flow rate of the liquid in the other two throats increased. At the next time step, when the pore space in section B was blocked, as shown in Figure 9(b), the flow velocity downstream this clogged space was reduced. This main flow path in the last step at t1 was blocked, as shown in Figure 10(a); thus, the main fluid steam was diverted to the path below. Figure 9. Pressure variation due to particle clogging at a) t1, b) t1 + 0.02 s, and c) t1 + 0.04 s Figure 10 . Flow path conversion with particle clogging at a) t1, b) t1 + 0.02 s, and c) t1 + 0.04 s Figure 11(a) shows the fluid velocity distribution in the cross-section at the narrowest point of the throat, as labeled in Figure 10(a) and Figure 11(b). Owing to particle clogging before the cross-section, the highest predicted velocity is in the middle path, and the lowest is in the bottom path. Therefore, St is highest in the middle path and lowest in the bottom path. At time t2, the bottom path was blocked before the middle path, even though the St value was higher in the middle path. This is because the pore spaces were blocked upstream of the high-St middle path; thus, fewer particles could approach the middle path compared with the bottom path. At time t3, even though the local number density of particles was lower, the particles blocked the middle path owing to the high St. These microscopic studies reveal that the flow conversion caused by particle clogging determines the efficiency of particle filtration downstream. Figure 11 . Velocity variation downstream of the particle 4 Conclusion In this study, the transport and retention of particles in a slurry passing through a porous medium were numerically studied using CFD-DEM. The flow of particles with diameters smaller than the throat size of a homogenous porous medium was modeled to study the hydrodynamics of fluid–particle flow in the porous medium, particularly when particle clogging occurred by particle bridging. In the simulation, when St was lower than 1, the particles passed through the porous medium with little difference between particles of different sizes. However, when St was greater than 1, the larger particles caused particle clogging due to bridging. The pressure and fluid velocity distributions were predicted to become heterogeneous when particle clogging occurred. This provides evidence for the ability of particle clogging to redistribute the pressure distribution within the system and vary the flow path at the microscale, thus reducing the permeability of the porous medium at the macroscale. In addition to the value of St in the pore space, the simulation shows that the particle concentration also affects particle clogging.