CFD–DEM Simulation of Sand-Retention Mechanisms in Slurry Flow

: The primary motivation of this paper is to investigate the sand-retention mechanisms that occur at the opening of sand ﬁlters. Various retention mechanisms under various conditions are explored that have a particulate ﬂow with a low concentration of sand particles (called slurry ﬂow) such as particle shape, size, and concentration. The computational ﬂuid dynamic (CFD)– discrete element method (DEM) model is applied to predict the retention mechanisms under steady ﬂow conditions of the well-bore. By using coupled CFD–DEM (CFD to model the ﬂuid ﬂow, and DEM to model the particle ﬂow), the physics involved in the retention mechanisms is studied. The coarse grid unresolved and the smoothed unresolved (reﬁned grid unresolved) coupling approaches implemented in STAR-CCM+ (SIEMENS PLM) are used to transfer data between the ﬂuid and solid phases and calculate the forces. The ﬁlter slots under investigation have different geometries: straight, keystone, wire-wrapped screen (WWS) and seamed slot and the particles are considered with different shapes and different aspect ratios and size distributions. The ﬂow regime is laminar in all simulations conducted. The CFD–DEM model is validated from the perspectives of particle–ﬂuid, particle– particle, and particle–wall interactions. Veriﬁcation of the CFD–DEM model is conducted by mesh sensitivity analysis to investigate the coupling resolution between the CFD and DEM. By simulation of numerous slurry ﬂow scenarios, three retention mechanisms including surface deposition, size exclusion, and sequential arching of particles are observed. However, the concentration of particles is too diluted to result in multiparticle arch formation. In the simulations, various conditions are tested to give us an insight into the parameters and conditions that could affect the occurrence of the retention mechanisms. As an example, the importance of the gravity force and interaction forces on retention mechanisms are conﬁrmed at the microscale in comparison with others forces involved in retention mechanisms such as the drag force, lift force, cohesive force, buoyancy force, and virtual mass force. C.F.L.


Introduction
Sand destroys oil production equipment, such as pipes and pumps [1]. Sand control technology reduces energy costs for the oil and gas industry [2]. By using sand filtration [3], with the same amount of energy, more oil will be produced [2]. Particle-retention mechanisms that occur at filter openings help with sand control [4] and supports sand filtration [5]. In this research, various retention mechanisms under various conditions are explored under particulate flows with low concentrations of sand particles. These conditions are applied through cases without and with porous media on top of the filter opening, various concentrations of particles, the effect of fluid flow, shape of the particles, particle size distribution, the role of the physical forces, as well as the filter opening shape and size. The following gaps served as motivation for this study. The majority of studies in this field only focused on investigation of the retention mechanisms in the porous media [6] and at pore-scale [7], not at the filter openings [8,9]. This research tried to fill that gap by engaging with the retention mechanisms at the filter opening and at microscale [8]. Similar studies have not investigated retention mechanisms that have fluids with high viscosity [10]. This study aims to offer an insight into the occurrence of the retention mechanisms in slurry flow with heavy oil as the carrier phase [11]. This study also makes a contribution to exploring new parameters that affect the particle-retention mechanisms, such as the role of physical forces [12,13]. Similar studies have not investigated retention mechanisms by numerical simulation [14]; instead, they have applied experimental studies [15]. This paper aims to offer an insight into these particle-retention mechanisms through numerical simulations.
The research outcomes of this study will be a computational fluid dynamics (CFD)discrete element method (DEM) model, capable of predicting occurrence of retention mechanisms in slurry flow, as well as a breakdown of the physics involved in the retention mechanisms involved in slurry flow.

Filtration
Filtration is an operation that separates solids from fluids (liquids or gases) by having a medium through which only the fluid can pass [5]. This medium is called a filter. In a filter, oversized solids in the fluid are retained while fluid and small particles pass through. Filtration happens in nature and engineering applications [3]. For instance, in the human body, kidney filtration removes waste from the blood; during water treatment, unwanted materials are removed by the filter medium; and in the oil industry, sand filters are used to prevent and control sand production, and to increase oil production while using the same amount of energy.
The filter open flow area-that is, the area of filter open to flow-and mesh size or opening size of the filter are the factors considered when selecting a filter. Solid particle buildup and plugging [15,16] at the entrance and throughout the length of the filter opening are the common failure [17] mechanisms in filters. Plugging is the reduction in the open flow area of the filter, mainly due to trapped solid particles. Plugging is an issue because it decreases the filter's efficiency. Particle buildup [18], another undesirable phenomenon, also happens due to the transport of very small particles in the system.

Retention Mechanisms
Retention mechanisms are important because they play key roles in preventing plugging and stop the solid particles from passing through the filter. Investigation of the conditions and factors that might affect retention mechanisms, such as force distribution (interaction forces, drag, buoyancy, pressure gradient, gravity, etc), particle shape and size, particle concentration, fluid properties (such as non-Newtonian viscosity [19]), and the geometry of the opening can help to improve the performance of the filters. Figure 1 presents the main particle-retention mechanisms at a granular scale while having particles larger than 100 µm [5]: 1.
Size exclusion mechanism results due to a large particle size in comparison to the constriction. It is also known as straining.

2.
Surface deposition mechanism occurs due to forces existing between the particles and the wall. The roughness [20] of the wall could cause surface deposition mechanism as well. 3.
The sequential bridging mechanism happens following surface deposition. 4.
Multiparticle bridging mechanisms could happen due to the high concentration of particles (mechanical bridging) or because of flow convergence at constriction, called hydrodynamic bridging. The problem of interest in this paper is the simulation of the particulate flow at the filter opening and investigation of the particle-retention mechanisms that have slurry flow with a special interest in the occurrence of multiparticle arching. The CFD-DEM model is applied to predict the retention mechanisms under steady flow conditions of the well-bore. By using coupled CFD-DEM (CFD to model the fluid flow, and DEM to model the particle flow), the physics involved in the retention mechanisms is studied.
The first phase of this work involves setting up a numerical model and conducting a rigorous verification and validation of the model to assure that the model is predictive from the perspective of particle-particle, particle-wall, and particle-fluid interactions. Following validation of the CFD-DEM model, best model setup practices for the particulate flow study are obtained. In various scenarios, the possibility of the retention mechanisms' occurrence and multiparticle arch formation are investigated by modeling the slurry flow with various test conditions.

Particulate Fluid Flow Modeling
This section is an overview of the particulate flow modeling approach applied in this paper. It emphasizes on the details of the mathematical and computational models behind the computational fluid dynamics-discrete element method (CFD-DEM technique) [21]. The DEM technique is a particular case of the Lagrangian method to model particular flow behavior.
In DEM, the motion of each particle is analyzed incorporating the fluid dynamic forces, the contact forces, and the moments due to the neighboring particles. Solid phase in the particulate flow can also be simulated on a discrete scale or a continuum scale. The discrete element method can model the discrete nature of the granular substance.
The CFD-DEM approach is the most attractive one for simulating particulate fluid flow at the particle scale. Some advantages of CFD-DEM over other approaches [22] are less computational cost in comparison with Lattice Boltzmann (LB)-DEM and DNS-DEM and the possibility of capturing the particle-particle interactions compared to TFM. The details of the advantages and disadvantages of these different models have been discussed in Zhu et al. [23]. DEM was originally proposed by Cundal et al. [24] in 1979. DEM simulations are important as they are capable of tracking the trajectory of the particles dynamically, considering various forces acting on individual particles. These features are not easily accessible from experimental techniques. The process of tracking particles is often slow and costly to achieve using experiments. Fundamental laws of Newtonian mechanics govern the interaction between discrete particles in DEM [24]. Researchers have used DEM to study particle packing, hopper flow, mixing, and granulation, which are all related to particulate flow over the past two decades [25]. Zhu et al. [26] published a comprehensive review of the theoretical developments on DEM and its applications. The governing equations in DEM for describing translational and rotational motions of individual solid particles are Newton's laws of motion. It relies on the motion of each particle individually as well as particle interactions.
Coupled CFD-DEM, the approach that is applied in this paper, was employed firstly by Tsuji et al. [27] and has gained popularity rapidly due to advances in computational power. The advantage of CFD-DEM is that it is applicable for various sizes of particles and geometries [28]. Considering the relative size between the particles and the numerical grid, there are two CFD-DEM methods: unresolved and resolved CFD-DEM approaches. The unresolved method is applicable for small particles relative to the grid used for the flow. In this approach, DEM is applied to calculate particle motion without resolving the detailed flow around each particle. Resolved methods are suitable for larger particles relative to the numerical grid and also for dealing with complicated geometries where small grids are necessary to increase the accuracy.
Sand filtration operation is a multiphase flow process of a granular-fluid system, where CFD-DEM is a promising approach [23]. In this paper, to simulate particle migration and retention mechanisms in sand screen devices, a fully coupled CFD-DEM model is applied. The concentration of particles in the particulate flow is an important factor in the selection of the methodology.

Overview of the DEM Theory
Discrete element method (DEM) model is a Lagrangian modeling methodology usually applied to model dense particle flow [24]. The governing equations in DEM for describing translational and rotational motions of individual solid particles are Newton's laws of motion. Both forms of Newton's second law of motion, linear momentum equation and angular momentum equation, are applied as the governing equations. DEM relies on the motion of each particle individually as well as particle interactions, and considers a group of computational points where each point is associated with a physical particle. The behavior of the particle is governed by the particle body force, surface force, and contact force [29]. Particles are modeled as the rigid body having two types of motions: translation and rotation. The relationship between the two particles is investigated by introducing two types of interaction forces: inelastic interaction force and elastic interaction force [29]. In the DEM approach, the contact force acting on a solid particle is calculated using a springdashpot system. Dashpots and springs are installed between the two particles to represent the inelastic collision between the two particles and the elastic interaction, respectively. Energy is dissipated to heat and sound during an inelastic collisions between particles. The distinct characteristic of the DEM is that the interparticle contact forces are included in the equations of motion. These forces cannot be ignored for flows including many interacting particles. A DEM calculation is computationally expensive, since it requires integration on small time steps to resolve surface contacts between particles. There is no mesh required for DEM as it is required for CFD and each particle is a computational point [30]. The DEM particles can have various shapes and volumes [31].
STAR-CCM+ from Siemens PLM is the commercial software mainly used in this research. Classical mechanics is used in STAR-CCM+ to model DEM and is based on soft-particle formulation, where particles can develop an overlap, but do not deform [31]. The overlap during the collision is decomposed into normal and tangential directions and introduced into contact models [32]. STAR-CCM+ provides three contact models to calculate the contact force in DEM: Hertz Mindlin, Linear Spring, and Walton Braun [31]. The Hertz Mindlin model has been applied in the majority of this paper.

Mathematical Model of Particulate Flow: Governing Conservation Equations
In this section, fluid flow conservation equations, solid flow conservation equations, and coupling approaches are discussed.

Fluid Flow Conservation Equations
To simulate fluid flow, the conservation equations for mass (continuity equation) and momentum (Navier-Stokes equations (NS)) are solved. In the case of particulate flow, the CFD cells (fluid mesh) are not fully occupied by fluid. NS considers that the volume occupied by the fluid in each cell is dependent on the volume of the solid particles in that cell. The continuity equation is written as follows: where f is the fluid void fraction and u is the fluid velocity.
Comparing Equation (1) to the incompressible continuity equation ∇ · u = 0, the presence of a time-dependent and space-dependent void fraction results in the nondivergence free velocity field to ensure the mass conservation. This means that the divergence of the velocity is not zero. The momentum equation is written as follows: where p is pressure, τ f is fluid shear stress tensor, ρ f is fluid density, and g is gravity. F p f , as a force, is the momentum transfer term, also called the momentum coupling term. n p is the number of particles and f p f ,i is the individual force acting on the particle i in the presence of the surrounding fluid including the drag f d,i , pressure gradient f ∇p,i , shear stress f ∇.τ,i , virtual mass f vm,i , Basset f B,i , Saffman lift force f Sa f f ,i , and Magnus lift force f Mag,i of particle i. ∆V in Equation (3) is the volume of the computational cell or elements. ⊗ is the tensor product of two vectors. This form of the momentum equation is generic and applicable to any concentration of particles and any Reynolds number in the nonturbulent flow regime and incompressible flow [32].

Solid Flow Conservation Equations
DEM equations track the position and velocity of each particle and account for the interactions between the particles [33]. The particles are allowed to have overlaps. The overlap must be less than 5% of the particle diameter and is decomposed into normal and tangential directions in contact models. Contact models simulate the contact forces in solid-solid collisions and accounts for elastic and inelastic (dissipative) collisions.
The conservation equation of linear momentum for a particle of mass m p is written as Equation (6) in the Lagrangian framework. The change in momentum is balanced by surface and body forces that act on the particle.
x i is the instantaneous linear displacement of the particle i, v i is the instantaneous velocity of the particle i, F s is the forces acting on the particle surface, and F b represents the body forces. F d is the drag force that is a particle-fluid interaction force. F p is the pressure gradient force also a particle-fluid interaction force. F vm is the virtual mass force that is an unsteady force. F l is the lift force (Magnus and Saffman), F b is the Basset force, and F τ is the shear stress force. F d , F p , F vm , F l , F b , and F τ are all surface forces. The momentum transfer from the continuous phase to the particle is represented by these surface forces F s in the coupling process. F g is the gravity force and F mg is the magnetic force-both are body forces. DEM modeling introduces extra body force F c representing interparticle interaction due to particle contacts with other particles and domain boundaries. F c is called the contact force. The contact force that represents interparticle and particle-boundary interaction is calculated for each particle considering the contacts that particle i has. This force is exerted at the contact point between the two particles or between the particle and the boundary.
f c,ij is the contact force between particle i and particle j. f c,iw is the contact force between particle i and wall w. The contact forces between the two particles or particle and wall are divided into elastic and dissipative forces that have normal and tangential components.
F nc is the noncontact force on particle i due to other particles such as electrostatic, van der Waals forces, or cohesive forces (such as liquid bridges).
Particles' angular momentum must also be conserved. The conservation equation of angular momentum for a particle is written as Equation (10) in the Lagrangian framework.
I i is the moment of inertia of particle i, and ω i is the angular velocity of particle i. T c,ij and T r,ij are contact torque and rolling friction torque applied on the particle i due to the collision with particle j. T c,iw and T r,iw are the contact torque and rolling friction torque applied on the particle i due to the collision with the wall section w.
To calculate the void fraction, there are various methods proposed in the literature such as the particle counting method, particle cloud, conservative particle cloud, moment preserving approach, two-grid formulation, and divided approach [32].

Coupling Approaches and Grid Resolution
The coupling refers to momentum, heat, and mass exchanges between the continuous phase and the dispersed phase. The important contribution to the particle-fluid momentum exchange is established by the drag force, dependent on the granular volume fraction. In contrast to the drag force acting on a single sphere, the granular volume fraction must be considered for the drag force calculation in the CFD-DEM coupling.

Zero-Way, One-Way, and Two-Way Coupling
There are three general ways of coupling approaches between the solids and the fluid flow: • In zero-way coupling or uncoupled, the dispersed phase and continuous phase evolve independently. • In one-way coupling, the dispersed phase feels the influence of the continuous phase and the effect of the dispersed phase on the fluid is not taken into account. • In two-way coupling, the effects of the dispersed phase on the continuous phase and the effects of the continuous phase on the dispersed phase such as displacement, interphase momentum, mass, and heat transfer are considered.
In all three coupling approaches, the particle-particle interactions and particle-domain interactions have been considered.

Coupling and Grid Resolution
From the perspective of the momentum exchange between fluid and particles, CFD-DEM coupling can be of two types: unresolved and resolved. This distinction is based on the method used for calculation of the forces on the particles. Depending on the method used, the resolution of the mesh must be adapted to the size of the particles.
• Unresolved: In this approach, DEM is applied to calculate particle motion without resolving the detailed flow around each particle. Submodels are used to calculate the forces on each particle and the same models are used as momentum sources in the CFD simulation. In the unresolved approach, the fluid flow is solved without resolving the flow around the particle. Generally, it means that the mesh is coarser than the particle size. Then, the coupling is established by submodels consisting of expressions for the drag force and other forces [32]. • Resolved: Resolved methods simulate the detailed flow around the particles and calculate the forces between the flow and particles without submodels. Resolved methods require a grid much finer than the particle size, and they can also deal with complicated geometries where small grids are necessary to increase the accuracy.
In STAR-CCM+, DEM simulations coupled with CFD are always unresolved. However, the software allows for continuous mesh refinement, including meshes that are finer than the particle size. The particle equations of motion are based on contact theory and Lagrangian tracking. Particle position is tracked by the Particle Centroid field function, and its shape is also a property of the particle. From these, contact detection algorithms are used to determine when a particle is contacting either another particle or a wall. This is represented by an overlap of the particle with the opposing particle. The force exerted between the particles is a function of the size of the overlap, as well as material properties such as Young's Modulus. Larger overlap means a larger repulsive force between the particles. This calculation is done on a particle by particle basis, with wall contacts using the surface mesh for its geometry. When DEM is coupled with CFD, one or more additional forces are applied to the particle from the CFD results. Usually, this includes a pressure gradient force and a drag force. These forces are calculated based on the data in the cell that the particle centroid is located in. These flow-based forces are applied to the particle as momentum source terms and, together with the contact forces, are used to determine the motion of the particle. If the two-way coupling model is activated, the equal and opposite momentum source term is applied in the CFD model to the cell in which the particle centroid resides. This unresolved approach will never show a flow redirecting around the surface of a particle because the particle shape is not resolved.
In STAR-CCM+, when the particle size is smaller than the grid size, the approach will be called the "coarse grid unresolved method" or "coarse unresolved" for short in this paper. The case with grid refinement-that is, when the particle size larger than the grid size-will be called "smoothed source on refined grid unresolved method" or "smoothed unresolved" for short in this paper.
The resulting impact on the flow is an effect similar to a porous zone in the cell where the centroid is. While this impact is confined to a single cell in coarse grids (coarse grid unresolved approach), in the case of fine grids (refined grid unresolved approach), this source term needs to be adjusted for two-way coupling. If a particle is larger than a cell, the forces applied to the particle can be much larger than the fluid in the cell will realistically experience. This can cause unphysical velocities and instability. In STAR-CCM+, this issue of the refined grid unresolved approach is solved by using the Volume Source Smoothing technique, which is activated in the Lagrangian solver settings. This technique creates a secondary grid composed of cells from the original grid, which is used exclusively for the two-way coupling. The smoothing technique ensures that the momentum source is spread across several cells occupied by the particle. During smoothing, the "volume fraction" scheme is used to calculate the void fraction and forces.

Verification and Validation
The CFD-DEM model should be verified and validated. The verified and validated model will be applied to investigate the particulate flow in a single opening of the sand screen applied in the oil industry for the purpose of sand control during oil production. This verification and validation process includes solving benchmarking problems to assure that the developed model is capable of capturing the physics of the particle-fluid, particle-wall, and particle-particle interactions.

Dimensionless Numbers for Particulate Flow Study
Prior to solving problems, it is worth mentioning the two main dimensionless numbers applied to study particulate flow. The first is the Reynolds number that is defined separately for dispersed phase (particles) as well as continuous phase, the second one is Stokes' number.
Equations (11) and (12) present the particle and fluid Reynolds' number, respectively, and Equation (13) presents the Stokes number. For large Stokes number values, the trajectory of the particles is independent of the fluid. As the value of the Stokes number reduces, the influence of the fluid on the particles' trajectory increases. For a Stokes number less than 1, the particles will follow the fluid streamlines closely. In the above equations, ρ f is the fluid density and ρ p is the dispersed phase (particle) density. V f and V p are velocity of the continuous phase and disperse phase, respectively. µ c is the continuous phase viscosity, D p is particle diameter, and L is the continuous phase characteristic length.

Flow through the WWS: CFD Model Verification and Validation
As an initial phase of the verification and validation process, computational fluid dynamic (CFD) was applied to investigate the details of the interaction between the nearwell-bore flow and the flow of fluids through and along the wire-wrapped screens, see Figure 2. In Figure 2, radial probe through the whole study domain is shown to collect pressure and velocity data in the simulations. The probe goes from consolidated porous media to unconsolidated porous media, then wire-wrapped screen, annulus, orifice and pipe (pie-shape), respectively. The study started with the single-phase flow simulation of water in WWS. The CFD model in the pipe and wire-wrapped screen governed by the Navier-Stokes equations was coupled to the porous media model outside of the well, governed by the modified Darcy's law: the Ergun equation. A laminar-steady-state CFD model was applied to investigate the flow behavior in the wire-wrapped screen in the steam-assisted gravity drainage (SAGD) [34] production well and study the effect of the particle size distribution in porous media on the flow in wire-wrapped screen. The particles in porous media were assumed to not be transportable. Consequently, flow in WWS and pipe was single-phase flow of oil. This work was a feasibility study on CFD application in wire-wrapped screens, as well as a parametric study that gave us an assessment of the practicality of CFD and the work-flow for the next complex models.
A comparative study of the coupled porous media with WWS model was conducted on "fine" and "coarse" oil sands using conditions that simulate SAGD production well with CFD. Analysis of the pressure and velocity fields along the WWS and pipe confirmed that having a single-phase flow of bitumen flow in WWS showed no dependence on the particle size distribution in porous media. This result, shown in Figure 3, confirms that by having various PSDs and no moving particles, PSD will not affect the flow in WWS, annulus, orifice, and pipe. The calculations conducted could be considered as validation for the CFD model. In Figure 3 from left to right, the pressure data was collected through the probe shown in Figure 2 and plotted. The lowest pressure belongs to pipe and highest pressure belongs to the consolidated porous media.
With no moving particles in porous media, it was expected that no pressure and velocity difference would be seen in the WWS, annulus, and pipe, as obtained by the simulations. Further, a mesh independence study was applied that resulted in a verified model. Details can be found in [35].

Particle-Fluid Interaction Validation
In this part of the study, a CFD-DEM model was set up and tested using three classic flow problems for the purpose of validation. A single particle (spherical shape) settling in a Newtonian fluid was simulated and it's terminal velocity resulted by CFD-DEM modeling was compared with the theoretical terminal velocity. A suspended (neutrally buoyant) particle was simulated by assigning the same density to the particle and the fluid.
Additionally, the CFD-DEM model was validated by the problem of the settlement of two particles in tandem with a Newtonian fluid.
In the first case study, a single particle was injected into the tank at an initial velocity of 0.014 m/s. It flowed down under gravity, drag, and buoyancy forces until it reached terminal velocity, where the forces were at equilibrium. Reynolds number of the particle was 72. The simulated terminal velocity by the CFD-DEM model was in good agreement with the theoretical particle terminal velocity obtained by balancing the gravity, buoyancy force, and viscous force. The terminal velocity obtained analytically was 0.02601 m/s. In addition, the terminal velocity of a nonspherical particle was obtained by the CFD-DEM model and compared with the calculated velocity by theory [36]. Details can be found in [35]. Figure 4 presents the single-particle terminal velocity simulated by the CFD-DEM model and the one calculated by theory. The fluctuations in the velocity are due to the wake street behind the particle. Having Re = 72, it is expected to have a vortex shedding behind the falling particle in the fluid [37,38].
The Strouhal number represents the ratio of inertial forces due to the local acceleration of the flow to the inertial forces due to the convective acceleration [39]. In Case 1, by the periodic motion of the wake behind the particle, the Strouhal number S tr is given by the following expression [39]: where f is the frequency of the flow oscillations ( 1 s ), L is the characteristic length (m), and U is particle velocity (m/s). In this case, L = 2d p = 0.005 [37], U = 0.02601 m/s, and f = 1 0.2 s . The Strouhal number value calculated based on the assigned values is 0.9612 for the falling particle in the quiescent fluid. In the second test case, a single particle was suspended. Having equal density for the continuous and dispersed phases, it was expected that the terminal velocity of the particle would be zero and the particle would become suspended. The particle was injected into the tank at 0.014 m/s and its velocity reached 7.7 × 10 −7 m/s. The difference of this value with zero could be justified with round-off error. If the simulation ran for a longer time, it would even get closer to zero. By having equal densities for the fluid and solid phases, the CFD-DEM model showed that the particle was eventually becoming suspended, which was equivalent to the zero terminal velocity.
The third problem was the case of settling two particles in tandem in the tank filled with quiescent water. The leading particle is injected and, after one second, the trailing particle is injected. It is expected that the trailing particle accelerated and moved faster toward the leading particle. The reason for this acceleration is the smaller drag force that the trailing particle feels due to moving fluid behind the leading particle [38], which results in a lower relative velocity. The Reynolds number of the trailing particle was expected to be less than the leading particle as well, due to the lower relative velocity. The CFD-DEM simulation confirms the following points: 1.
The velocity of the trailing particle was increasing at a higher rate (accelerating) and it approached the leading particle, as expected.

2.
Reynolds of the trailing particle was less than the leading particle due to the lower relative velocity between the continuous phase and the trailing particle.

3.
The drag force on the trailing particle was less than the drag force on the leading particle due to the lower relative velocity between the fluid and the trailing particle.
The outcome of working on these three scenarios was the validation of the developed CFD-DEM model from the perspective of particle-fluid interaction.

Particle-Wall Interaction Validation
In this part, the goal was to track the trajectory of a falling particle bouncing on a wall simulated by the CFD-DEM model and to compare the results with the experiment.
Considering the quiescent fluid of silicon oil in a tank made of glass, a steel sphere fell in oil toward the bottom glass wall. The mass density of the steel sphere was ρ p = 7800 kg/m 3 . Young's modulus of elasticity E was 214 × 10 9 Pa and Poisson's ratio v was 0.3. The mass density of silicon oil was ρ f = 970 kg/m 3 and the dynamic viscosity was µ = 0.1 Pa · s (at T = 20 • C). The size of the steel sphere was d p = 0.0053 m, Re number was 30, and St number was 55. According to the available information, three problems were investigated as follows: 1.
The falling particle trajectory was simulated with the CFD-DEM model and compared with the experiment [40].

2.
Domain sensitivity analysis was performed.

3.
Mesh sensitivity analysis was applied to compare the effect of the coupling resolution (refined grid unresolved vs. coarse grid unresolved).
The outcome of working on these three scenarios was validation of the developed CFD-DEM model from the perspective of particle-wall interaction. Details can be found in [35].
Additionally, polyhedral particle-wall interaction was investigated through investigation of the trajectory of a polyhedral particle bouncing on the wall.

Particle-Particle Interaction Validation
The problem of drafting-kissing-tumbling of two particles [41] was replicated with the CFD-DEM model to investigate the capability of the model to capture the physics involved in particle-particle interaction. The benchmark problem considered here concerned the simulation of the motions and interaction of two falling identical balls in a vertical tank with a square cross-section. The computational domain size was 0.1 cm × 0.1 cm × 0.4 cm. The diameter d of the two balls was 1 6 mm. The initial translational and angular velocities of the balls were zero. The density of the incompressible fluid was ρ f = 1000 kg/m 3 and the density of the balls was ρ s = 1140 kg/m 3 . The fluid was quiescent. The boundary condition for the tank walls and balls was no-slip. Fluid had a viscosity of 0.01 poise (0.001 Pa · s). The movement of the two balls was simulated and the distance between them was measured and compared with the direct numerical simulation (DNS) results [41]. For verification and validation of the model, different time steps and grid spacing were tested to obtain the distance between the two particles at drafting, kissing, and tumbling stages, then the simulation results were compared with the DNS results.
In this validation case, the effect of the coupling resolution, also known as mesh sensitivity analysis; first-and second-order time discretization schemes; first-and second-order convection discretization schemes; various volume fractions of each cell in the coupling stage; and the effects of various CFD and DEM time steps were tested.
According to the simulation results, the predicated distance between the two particles in three stages of drafting-kissing-tumbling was in good agreement with the DNS results. It confirmed the capability of the CFD-DEM model of capturing the physics involved in the fluid-particle interaction and particle-particle interaction, see Figure 5. In addition, mesh sensitivity analysis was completed by comparing numerous refined grid unresolved and coarse grid unresolved cases. It was confirmed that while the coarse grid unresolved cases follow the trend of the experimental results, the refined grid unresolved cases were in better agreement with the experimental results at the scale of this research, which was microscale.
The verification and validation of this model have provided an understanding of simulating the particle-particle interaction, particle-wall interaction, and particle-fluid interaction, considering spherical and polyhedral particle shapes. Details can be found in [35].

Recommendations on Model Setup for Particulate Flow Study
The first phase of this work involved setting up a numerical model and conducting a rigorous verification and validation of the model to assure that the model was predictive from the perspective of particle-particle, particle-wall, and particle-fluid interactions. Literature was investigated to apply the right particulate flow modeling approach at the microscale, considering the computational resources. It was found that by the CFD-DEM model, the collision of hundreds of particles and the interparticle forces could be simulated. The CFD-DEM models were set up and validated from the perspectives of particle-particle, particle-wall, and particle-fluid interactions with the benchmark problems as well as the available experimental work from the collaborating labs and literature. Following the verification (mesh independence analysis) and validation practices, a list of best practices was achieved to set up the model and study particulate flow and particle-retention mechanisms at filter openings as follows:

1.
A mesh should be designed in a way that gives adequate resolution in regions where spatial gradients are high, such as slot entrance. Considering the computational cost and accuracy, the best practices of mesh resolution for a steady-state particulate flow case was achieved by smoothed unresolved coupling, where the size of the particle was at least twice the mesh size. Moreover, mesh resolution for a transient particulate flow case was achieved by smoothed unresolved coupling, where the size of the particle was at least five times larger than the mesh size. The quality of the robust simulation model was guaranteed by grid independence and error estimation.

2.
Any recirculation across a flow boundary should be avoided.

3.
The mesh should be aligned with the flow to improve the accuracy and increase the rate of convergence.

4.
A grid sensitivity study with two or more meshes must be applied.

5.
The unsteady simulation had to start from a converged steady solution and then convert to a transient one.

6.
Second-order time discretization was preferred, and a spatial discretization that was based on a second-order upwind differencing scheme. 7.
The segregated solver was preferred to the coupled solver for transient cases considering the accuracy and the computational cost. 9.
Freezing the solvers helped initial simulation stability, computational efficiency, and reducing total run time specifically early in the simulation. Convergence issues could arise because of the lack of a good physical initial condition for transient simulations. 10. The preferred boundary condition combination was to use no-slip wall and symmetry for the domain boundaries. 11. The preferred mesh is hexahedral rather than polyhedral considering accuracy and computational cost. 12. A velocity inlet boundary condition type for the upstream, and a pressure outlet boundary condition for the downstream boundary were preferred. 13. It was revealed that a Young's modulus lower than the actual one could reduce simulation time as DEM time-step was dependent on Young's modulus. A reduction of up to three orders of magnitudes could be applied with sensitivity analysis. 14. The Hertz-Mindlin contact model could predict well how particles interact during the collision. 15. The boundary condition was set at continuum level by default and had to be set up for each individual boundary of particles. 16. Maximum porosity could be specified as the limit after which injection of new particles is stopped. This feature could help build the porous regions with particles. 17. Two-way coupling is recommended to study retention mechanisms, where each phase is affected by the other phase-for example, if drag force is included in the simulation, a momentum source is added to the continuous phase. The reason for using two-way coupling is all the interactions involved in the collision are included in the simulation. 18. Volume fraction of the disperse phase in the continuous phase is dependent on cell size, where the maximum volume fraction is 1 in case it has larger particles than cells. The solution could be unstable if the disperse phase produces high sources of momentum. Sources could be smeared/smoothed in more cells through DEM solver to avoid instability. The maximum volume fraction specified in the simulations in the case of smoothed unresolved coupling was 0.98. Using 1 resulted in the simulation's instability. 19. By applying the soft-sphere model, the CFD-DEM model was predictive. In the softsphere model, particles could overlap but not deform. The soft-sphere formulation is based on the contact forces that are established during contact between particles [31]. In this model, contact force is proportional to the overlap and multiple contacts are allowed at the same time.

Sand-Retention Mechanisms at the Openings of the Sand Filters
Simulation of the sand-retention mechanisms (surface deposition, size exclusion, sequential arching, and multiparticle arching) and the investigation of the factors that affect these mechanisms will be helpful to enhance the criteria for sand screen selection and sand screen design in the oil industry [4,42,43]. Sand-retention mechanisms can be explored through simulation of the particulate flow around the sand screen device. Sand screens act as an obstacle to prevent the solid particles from flowing into the oil production well and rely on the four mentioned mechanisms.
Four types of sand-filter slots (openings) are examined in this research: straight, keystone, wire-wrapped screen (WWS), and seamed slot ( Figure 6). The geometry of these slots are different and those differences could affect their sand-retention performance.

Investigation of the Sand-Retention Mechanisms with Slurry Flow
A slurry is a diluted mixture of a fluid and solids denser than the fluid suspended in the fluid phase [44]. The size of solid particles may vary from 1 µm up to hundreds of millimeters. In this case, the solid particles are sand particles at the scale of 200 µm. The fluid phase is the carrier phase transporting solids. In slurry flow, solid phase or dispersed phase could have various concentrations of particles resulting in various slurry concentrations. The slurry concentrations normally range from 6% to 40% for sand slurries [44].
The slurry flow can occur during the initial phase of a SAGD process, called warmup [18], followed by formation of a loosened packed-bed [45][46][47], called the unconsolidated packed-bed, in the gap between the sand filter and the oil-sand reservoir.
The objective was to investigate the sand-retention mechanisms that occurred with the slurry flow with a special interest in multiparticle arching occurrence. Various parameters and conditions of the slurry flow were tested. The geometries and dimensions of the simulation domain, shown in Figure 6, were chosen to replicate the shapes of a slurry flow experiment [48] for qualitative comparison. Table 1 presents the computational and mathematical setup of the various slurry cases studied. The slurry flow cases considered here consist of Athabasca oil and sand particles on top of a single filter opening. Considering the low velocity of the fluid and particles as well as the low concentration of the particles, various time-steps were tested for the simulations and it was found that by setting a DEM time-step that is 10 times larger than the CFD time-step, the overall computational time was much faster without the loss of accuracy.
Particle-fluid interaction force stands for summation of the drag, buoyancy, virtual mass, and lift force in this study. The cohesive force was tested and was found to be negligible at the scale of tested particles (microscale).
Particle normal and tangential stiffness ( N m ), filter normal and tangential stiffness ( N m ), rolling resistance coefficients, normal and tangential damping coefficients, as well as the coefficients of friction were applied directly from the STAR-CCM+ database default values for material properties for steel (filter walls) and glass (particles). None of these values were specified explicitly in the software.

Slurry Cases without and with Porous Media
A structured porous region on top of the slot could have affected particles' interactions and consequently the sand-retention mechanisms. In order to observe the occurrence of various retention mechanisms, one case consists of an open flow area in front of the slot and the other case consists of a structured porous medium with fixed circular obstacles. This case condition and simulation is presented in Table 1. Figure 7 shows two sand-retention mechanisms that were observed: surface deposition and size exclusion in this case. Here is a list of observations by simulating cases with and without the porous region:

1.
A porous region on top of the slot reduced the number of particles passing through the slot in the specified time compared with the case without the porous region.

2.
Surface deposition and size exclusion were observed as retention mechanisms.

3.
There was no formation of a multiparticle arch in either case. It was obvious from the simulations that the concentrations of the particles are too diluted to result in multiparticle arch formation at the slot.

Slurry Cases with Various Concentration of Particles
Slurry concentration is the ratio of the volume of solids and the total volume of the mixture, usually given as a percentage. In these simulations, a range of sand particles' concentrations (10%, 20%, 25%, 30%, 35%, 36%, 38%, 40%) were tested to investigate the effect of particles' concentration on sand-retention mechanisms, as summarized in Table 2. The following findings were observed from the simulated scenarios: 1.
The number of particles passing through the slot within the time specified (11 DEM timesteps = 1.1 s) initially increased at the beginning of the flow and then reduced when particles started to deposit, as shown in Table 2.

2.
Surface deposition was clearly observed, see Figure 7. A size-exclusion mechanism occurred at the pores in porous media, as the size of the particles was larger than the pore size at some locations.

3.
Despite various concentrations of particles within the range for slurry flow, the multiparticle arch and sequential arching did not occur in any case with spherical particles.

Slurry Cases with and without Fluid Flow
To assess the effect of the fluid flow on the retention mechanisms, several simulation scenarios with flow of different fluids were conducted with and without fluid (i.e., without any fluid interaction, as in vacuum). Three fluids were tested: air (µ f = 1.849 × 10 −5 Pa · s and ρ f = 1.184 kg/m 3 ), Athabasca oil (µ f = 0.0136 Pa · s and ρ f = 915.2 kg/m 3 ), and water (µ f = 8.8871 × 10 −4 Pa · s and ρ f = 997.561 kg/m 3 ).
In the first case, the fluid solver was forced to stop-also called the frozen or no-fluid case-to remove the fluid stress, normal, and shear on particles and to investigate the occurrence of the sand-retention mechanisms in the extreme case of no fluid, as in a vacuum. The concentration of the slurry flow was 35% for these scenarios. The problem setup is presented in Table 1 and the corresponding nondimensional numbers are listed in Table 3. The study was conducted on WWS with and without porous media on top of the slot. Here is a list of findings that were observed out of the scenarios with and without fluid flow:

1.
For the case with the most viscous fluid (Athabasca oil), less particles passed through the slot in the specified time, Table 4.

2.
The presence of viscous fluid (water and oil) slowed down the flow of the particles, compared with the cases of no fluid and air, resulting in less particles passing through the slot per unit of time.

3.
There was no multiparticle arch formed in any case. The reason for the difference between the two rows in Table 4 is the existence of the porous media on top of the slot. The porous media could act as a barrier and sometimes as a trap for particles; accordingly, less particles can pass through the slot in the case of porous media.

Slurry Cases with Spherical and Polyhedral Particles
In all previous cases, particles were spherical. In this section, polyhedral particles, with aspect ratios larger than 1 and with sharp corners, were tested with slurry flow. Three sizes were tested for the particles within the possible range, corresponding to 3 to 5 times smaller than the slot opening (200, 275, and 375 microns). The slot shape is WWS and the cases were tested with and without porous media on top of the slot. The problem setup is the same as that presented in Table 1, except that the particles are polyhedral with aspect ratio equal to 1.6. For the nonspherical scenarios, the same force models as presented in Table 1 were applied.
The following findings were observed in this study: 1.
Surface deposition was observed on the filter surface and at the pores in porous media. Size exclusion and sequential arching or sequential bridging were observed in pores as retention mechanisms, shown in Figure 8.

2.
In the case of the uniform-size polyhedral particles, a multiparticle arch was also not formed at the slot.

Slurry Cases with Spherical and Polyhedral Particles
In all previous cases, particles were spherical. In this section, polyhedral particles, with aspect ratios larger than 1 and with sharp corners, were tested with slurry flow. Three sizes were tested for the particles within the possible range, corresponding to 3 to 5 times smaller than the slot opening (200, 275, and 375 microns). The slot shape is WWS and the cases were tested with and without porous media on top of the slot. The problem setup is the same as that presented in Table 1, except that the particles are polyhedral with aspect ratio equal to 1.6. For the nonspherical scenarios, the same force models as presented in Table 1 were applied.
The following findings were observed in this study: 1.
Surface deposition was observed on the filter surface and at the pores in porous media. Size exclusion and sequential arching or sequential bridging were observed in pores as retention mechanisms, shown in Figure 8.

2.
In the case of the uniform-size polyhedral particles, a multiparticle arch was also not formed at the slot.

Particle Build-Up Comparison with Experiment from Literature
The results from the uniform-size polyhedral slurry flow simulation on a WWS opening revealed particle build-up around the slot entrance, Figure 9. They were compared with the particle shadowgraph velocimetry (PSV) experimental results conducted by Kinsale et al. [48,50] on a straight slot. Similar surface deposition and particle build-up around the slots were observed in both simulation and experiment. It could confirm the validity of the simulation results. The scale of the porous media is similar in the experiment and simulation. The experimental particle scale was 40 µm and the scale of the particle used in CFD-DEM simulation was 100 µm. Surface deposition was also observed in both the simulation and experiment.
Despite the significant difference between the particle sizes in simulation and experiment, the retention mechanisms are similar, and the CFD-DEM model can capture the same retention mechanisms, i.e., particle build-up and surface deposition. It can be concluded Figure 8. Three sand-retention mechanisms were observed: surface deposition, size exclusion, and sequential arching in this case with slurry flow of polyhedral particles and porous media on top of the slot (CFD-DEM simulation).

Particle Build-Up Comparison with Experiment from Literature
The results from the uniform-size polyhedral slurry flow simulation on a WWS opening revealed particle build-up around the slot entrance, Figure 9. They were compared with the particle shadowgraph velocimetry (PSV) experimental results conducted by Kinsale et al. [48,50] on a straight slot. Similar surface deposition and particle build-up around the slots were observed in both simulation and experiment. It could confirm the validity of the simulation results. The scale of the porous media is similar in the experiment and simulation. The experimental particle scale was 40 µm and the scale of the particle used in CFD-DEM simulation was 100 µm. Surface deposition was also observed in both the simulation and experiment.
Despite the significant difference between the particle sizes in simulation and experiment, the retention mechanisms are similar, and the CFD-DEM model can capture the same retention mechanisms, i.e., particle build-up and surface deposition. It can be concluded that the difference between the geometry of the two slots had no significant effect on the retention mechanisms and particle build-up before the slots.

Slurry Cases with Particle Size Distribution
The particle size distribution (PSD) of a slurry defines the relative amount, typically by mass, of particles present in the slurry according to size. In the simulations, PSD was tested to investigate the effect of nonuniform particle distribution on sand-retention mechanisms. Particle size distribution was considered as D 10 = 50 µm, D 50 = 200 µm, D 90 = 400 µm. The rest of the problem setup is similar to that presented in Table 1.
The following is a list of the observed findings: 1. Similar to the previous scenario, three sand-retention mechanisms were observed: surface deposition, size exclusion, and sequential arching, as shown in Figure 10.

2.
The average number of particles passing through the slot at the specified equal time was less for the case with PSD than for the uniform case.

3.
With polyhedral particles and PSD, multiparticle arch did not form on the slot opening for cases with and without porous media. that the difference between the geometry of the two slots had no significant effect on the retention mechanisms and particle build-up before the slots.

Slurry Cases with Particle Size Distribution
The particle size distribution (PSD) of a slurry defines the relative amount, typically by mass, of particles present in the slurry according to size. In the simulations, PSD was tested to investigate the effect of nonuniform particle distribution on sand-retention mechanisms. Particle size distribution was considered as D 10 = 50 µm, D 50 = 200 µm, D 90 = 400 µm. The rest of the problem setup is similar to that presented in Table 1.
The following is a list of the observed findings: 1. Similar to the previous scenario, three sand-retention mechanisms were observed: surface deposition, size exclusion, and sequential arching, as shown in Figure 10.

2.
The average number of particles passing through the slot at the specified equal time was less for the case with PSD than for the uniform case.

3.
With polyhedral particles and PSD, multiparticle arch did not form on the slot opening for cases with and without porous media.

Slurry Cases: Testing Physical Forces Exclusion
There are various forces playing roles in the simulation of particulate flow passing through the filter slot, which are drag force (fluid stress), buoyancy force, gravity force, virtual mass force, cohesive force, lift force, and solid interaction forces (particle-particle, particle-domain). The summation of the three main forces acted by the fluid on the particles: drag force, buoyancy, and virtual mass were considered as the particle-fluid interaction force. The problem setup is similar to that presented in Table 1. Forces were excluded one by one and the cumulative number of particles passing through the slot in the specified time interval (10 DEM time-steps = 1 s) was tracked, see Table 5. Table 5. Cumulative number of particles passing through the slot in the specified time interval (1 s) while a force is removed from the simulation.

Omitted Force
WWS Without Porous Media WWS with Porous Media drag force 7 3 buoyancy force 14 5 gravity force 0 0 virtual mass force 12 5 cohesive force 13 5 lift force 12 6 particle-fluid interaction force 6 4 particle-particle interaction force 6 3 particle-wall interaction force 5 3 The main forces that caused particles to pass through the slot were gravity and interaction forces (particle-wall interaction, particle-fluid interaction, particle-particle interaction) since, by omitting them, there are less particles passing through the slot. Note that the particle-fluid interaction combines the main fluid forces (drag, buoyancy, and virtual mass force). In the slurry case, the particle concentration is low and there was no formation of multiparticle arch in any case.

Slurry Cases: Testing Various Slot Shapes
In this series of simulations, four different slot geometries were tested under conditions similar to those presented in Table 1. The objective was the investigation of the possibility of multiparticle arch formation with different slot geometries, Figure 6. These cases were tested with a porous region on top of the slot opening and with/without fluid. The simulation results confirmed that there was no multiparticle arch formed, neither mechanical nor hydrodynamic. Various geometries (keystone, WWS, straight, and seamed slots) could not help with multiparticle arch formation while having slurry flow. Table 6 presents the average number of particles passing through each slot in the specified time interval, which is 10 DEM time-steps = 0.1 s.

Slurry Cases: Testing Various Slot Widths and Lengths
In these simulations, various slot widths and lengths were tested to investigate the possibility of the multiparticle arch formation with different slot widths and lengths. Values of 6d p , 5.5d p , 5d p , 4d p , 3d p , 2.5d p , and 2d p were tested for the length and width of the slots, where d p is the particle diameter. The flow regime was laminar in all cases studied. In any of the cases, no mechanical or hydrodynamic multiparticle bridges were formed. Cases were also explored, whereby particles were deliberately injected into the arch-shape locations initially to support the arch formation. This strategy did not help either. The results confirmed that the slurry flow does not result in multiparticle arch formation.
Three sand-retention mechanisms were observed: surface deposition, size exclusion, and sequential arching.

Conclusions
The first phase of this work involved setting up a numerical model and conducting a rigorous verification and validation of the model to assure that the model was predictive from the perspective of particle-particle, particle-wall, and particle-fluid interactions.
Following validation of the CFD-DEM model, best practices of model setup for particulate flows were obtained and described in detail. In various scenarios, the possibility of the retention mechanisms' occurrence and multiparticle arch formation were investigated by modeling the slurry flow with various test conditions. Three retention mechanisms including size exclusion, surface deposition, and sequential arching were observed. However, no multiparticle arch was formed at the slot. The conclusion is that the multiparticle arching phenomenon does not occur during the slurry flow in various types of sand screen openings under the various conditions tested. These results occurred even for a narrow slot.
It was clear from the majority of simulated slurry scenarios that the concentration of the particles was too diluted to result in multiparticle arch formation. The investigations on slurry cases confirm the finding that the stress from the packed-bed of solid particles on top of the slot opening might be essential to form the multiparticle arch as the efficient retention mechanism.
The simulation results could also suggest the importance of the gravity, particlefluid interaction force, and particle-particle interaction force for modeling the transport of particles around and into the slot in the slurry flow. The current work provides a platform to start the next phase of the research-modeling and investigating of retention mechanisms, particularly multiparticle sand arch formation with a packed-bed of sand particles on top of the slot and in the gap between the sand filter and oil-sand reservoir.