Bonding Strength Effects in Hydro-Mechanical Coupling Transport in Granular Porous Media by Pore-Scale Modeling

The hydro-mechanical coupling transport process of sand production is numerically investigated with special attention paid to the bonding effect between sand grains. By coupling the lattice Boltzmann method (LBM) and the discrete element method (DEM), we are able to capture particles movements and fluid flows simultaneously. In order to account for the bonding effects on sand production, a contact bond model is introduced into the LBM-DEM framework. Our simulations first examine the experimental observation of “initial sand production is evoked by localized failure” and then show that the bonding or cement plays an important role in sand production. Lower bonding strength will lead to more sand production than higher bonding strength. It is also found that the influence of flow rate on sand production depends on the bonding strength in cemented granular media, and for low bonding strength sample, the higher the flow rate is, the more severe the erosion found in localized failure zone becomes.


Introduction
Sand production is commonly observed during the extrusion of hydrocarbons from reservoirs, which makes a great of trouble in oil or gas production [1].It was reported that 70% of the world's gas or oil reservoir has the potential risk to produce sand during their production life [1,2].As a result, the petroleum industry spends millions of dollars each year to find ways to solve this problem [3].
In the past few decades, some achievements on sand production have been made by experimental studies, among which two mechanisms about sand production are convincing [4,5].Firstly, the initial sand production is evoked by localized failure near the cavity owing to concentration of external stress, which is called as mechanical instability.Secondly, this localized failure zone is eroded by flow, which is named hydro-mechanical instability.Furthermore, these two mechanisms are often coupled with each other, which leads to the complexity of sand production.However, owing to the limitation in observing and measuring the micro-information about fluid and sand grains interaction by experiment, the micro underlying physics of these mechanisms is still not clear.
Therefore, numerical studies are of great importance.For a reliable numerical model for sand production, three aspects should be considered: an efficient method to capture the elastic and plastic behavior of granular media, a good description of fluid flow, and an accurate way to consider the fluid-solid interaction.In the previous numerical studies of sand production, continuum model and discrete element method are two popular approaches [3].For continuum method, granular media is regarded as a deformable pore-elastic or pore-elastoplastic solid described by the constitutive equation [6].However, continuum models have difficulty completely capturing the performance of granular media due to its inherent discrete nature.Meanwhile, little micro information can be obtained in continuum method to explore the mechanism at particle scale [7].As an alternative method, the discrete element method (DEM) [8] is an efficient way to overcome those shortcomings.Firstly, DEM can capture the discrete nature of rock by treating it as an assembly of discrete particles.Secondly, DEM is based on the particle scale and the influence of micro parameters on the macroscopic behaviors can be considered, which is important for studying mechanism of sand production.Owing to these advantages, DEM was widely used to study sand production [9][10][11][12][13].
In DEM simulation of sand production, the description of fluid dynamics plays an important role.In earlier works, flow in granular media was described by Darcy's law, which was coupled with DEM by some empirical formulas [9,10,13].However, in Darcy-DEM models, flow in granular media cannot be calculated accurately, so CFD-DEM model was developed, where Navier-Stokes equations based on the local mean values over a computation cell were solved [11,12].Zhou et al. [11] employed this approach to simulate the liquid-induced erosion in a weakly bonded granular media and captured the main features of the sand erosion.Compared with Darcy-DEM model, some progress has been made in CFD-DEM model, but empirical equations are again usually introduced to calculate hydrodynamic force on solid.In addition, due to the fact that hydrodynamic force acting on sand plays a critical role in sand production, the fluid-particle interaction must be evaluated at the particle scale or smaller [14], which means that CFD-DEM model is also not accurate enough to explore the underlying physics of sand production, especially about the fluid and sand grains interaction.
In the past three decades, the lattice Boltzmann method (LBM) has gained much popularity due to its efficiency in dealing with complex-boundary problems [15][16][17][18].It has also been coupled with DEM to simulate particle-fluid system at pore scale [19][20][21].In LBM-DEM model, the flow is calculated at pore scale, so that assumptions about the dependence between porosity and permeability are no longer needed.Therefore, the LBM-DEM model is outstanding among other numerical methods, and has been applied to study sand production process [7,14,22,23].
Boutt et al. [7] and Ghassemi et al. [22] simulated sand production in two dimensions (2D) with a coupled LBM-DEM model, and discussed the influence of the confining pressure and flow rate on sand production.In their models, the granular media is un-cemented, so there is no bonding effect between particles, and the only resistance for particle motion is friction.However, in cemented reservoir, the bonding effect is a primary property determining the macroscopic behavior of granular media, such as sandstone, and the failure of cemented granular media is generally owing to the bond breakage [24].It was also reported that behavior of un-cemented granular media is evidently distinct from cemented granular media owing to the bonding effects between grains [25], so the conclusion obtained from un-cemented granular media may be not suitable for cemented granular media.Moreover, only limited literatures were found to simulate sand production of cemented granular media at pore scale, although some sand production experiments have been carried out for a better understanding of sand production process in cemented reservoir [2,26].To bridge this gap, we study sand production in cemented reservoir, with special focus on the influence of bonding effect.In this paper, the LBM and DEM in a two-dimensional domain are coupled to simulate fluid-solid system, and in order to consider the effect of bonding strength, a contact bond model is incorporated in LBM-DEM framework.
The rest of this paper is organized as follows.Firstly, the numerical methods are introduced in Section 2, which includes the lattice Boltzmann method, discrete element method with contact bond model, and coupling of these two methods.In Section 3, the LBM-DEM code is validated by simulating the sedimentation of circular particles in a channel.Section 4 is the numerical results and discussion, which contains two parts: biaxial compression simulation and sand production simulation where the influence of bonding strength on sand production is discussed.Section 5 concludes the paper.

Numerical Methods
This section is a brief introduction of numerical methods used in our simulation: the lattice Boltzmann method for fluid phase and the discrete element method with contact bond model.The approach to couple these two methods is described at the end of this section.
The basic variable in LBM is density distribution function f i .A widespread LBM implementation in the literature is the lattice Bhatnagar-Gross-Krook (LBGK) model, where the collision operator is simplified as a linearized version by assuming that the collision operator relaxes the local particle distribution functions at a single rate [15].The evolution equation of "LBGK" model can be written as [15]: where f i is the distribution function in the ith discrete velocity direction e i , f eq i is the corresponding equilibrium distribution function, δ t is the time step, and τ is the relaxation time related to the fluid kinematic viscosity: υ where c " δ x {δ t is the lattice speed with δ x representing the lattice size, i.e., the space step.For a two-dimensional nine-speed (D2Q9) model, they are e o " p0, 0q , e 1 " c p1, 0q , e 2 " c p1, 1q , e 3 " c p0, 1q , e 4 " c p´1, 1q , e 5 " c p´1, 0q e 6 " c p´1, ´1q , e 7 " c p0, ´1q , e 8 " c p1, ´1q . (3) The equilibrium distribution for two-dimensional nine-speed (D2Q9) model can be given as: where the weighting factors are: ( The fluid density and momentum can be obtained from the discrete distribution function, and the pressure (p) is given by: p " c 2 ρ{3 (8)

Discrete Element Method (DEM)
Discrete element method (DEM) was proposed in 1979 [8] and has been widely used in geomechanics.In our model, the classic linear spring-dashpot contact model is used [8].The force-displacement relationship is given by: F n " ´kn α ´Cn ν n F s " ´ks η ´Cs ν s (10) where F n is the normal contact force, F s the tangential contact force, k n the normal spring stiffness, k s the tangential spring force, C n the normal damping, and C s the tangential damping.The shear slider will work at the contact point, when F n and F s satisfy the following inequality: where µ is the coefficient of friction.The motion of the particle is described by: where m is particle mass, F the total force acting on it, I the rotation inertia, and M is the total torque.In this work the explicit Verlet algorithm is applied, which is often used in DEM and MD simulation [42].LBM-DEM coupling at each time step is realized by computing the fluid field by LBM firstly and then the hydrodynamic force and contact force between solids are calculated.Finally, different external forces are added together to update the particle position by integrating the above differential equation.The computation of LBM-DEM coupling is shown in Figure 1.

Discrete Element Method (DEM)
Discrete element method (DEM) was proposed in 1979 [8] and has been widely used in geomechanics.In our model, the classic linear spring-dashpot contact model is used [8].The force-displacement relationship is given by: where n F is the normal contact force, s F the tangential contact force, n k the normal spring stiffness, s k the tangential spring force, n C the normal damping, and s C the tangential damping.
The shear slider will work at the contact point, when n F and s F satisfy the following inequality: where μ is the coefficient of friction.The motion of the particle is described by: where m is particle mass, F the total force acting on it, I the rotation inertia, and M is the total torque.In this work the explicit Verlet algorithm is applied, which is often used in DEM and MD simulation [42].LBM-DEM coupling at each time step is realized by computing the fluid field by LBM firstly and then the hydrodynamic force and contact force between solids are calculated.Finally, different external forces are added together to update the particle position by integrating the above differential equation.The computation of LBM-DEM coupling is shown in Figure 1.In order to consider the bonding effect in cemented granular media, the contact bond model is incorporated into the DEM, which is simple but can efficiently capture bonding effect between particles [43].In this model, the bond behavior is decided by the normal bonding strength n R and the tangential bonding strength t R , and t R is set equal to n R .For two contacting particles, the failure envelope can be summarized by: In order to consider the bonding effect in cemented granular media, the contact bond model is incorporated into the DEM, which is simple but can efficiently capture bonding effect between particles [43].In this model, the bond behavior is decided by the normal bonding strength R n and the tangential bonding strength R t , and R t is set equal to R n .For two contacting particles, the failure envelope can be summarized by: Computation 2016, 4, 15 where N s is the normal strength, and T s is the tangential strength.The detail about this model can be found in [43].

Fluid and Solid Interaction
Midway bounce-back boundary condition is applied in this work, which is easy to be implemented and widely used in particle suspension model [19].For a moving wall, an additional term is added into the classic bounce-back procedure to keep nonslip boundary condition: where u w is the velocity on the midpoint of the boundary link, x f is the fluid node next to the solid, e α the direction from fluid node to solid node, e α its opposite direction, and r f is the post-collision distribution function.
The calculation of hydrodynamic force is a major topic in particle suspension simulation.Ladd proposed a shell model basing on momentum exchange, where the hydrodynamic force between fluid and solid is calculated on each boundary link [19]: The total hydrodynamic force acting on the particle is the sum of δF w px w , e α q on each boundary link of the particle surface: The corresponding total torque can be given by: where R is the mass center of the particle.In Ladd's model [19], the fluid exists inside and outside of the particle boundary, so it is actually a shell model and the momentum exchange is calculated for both interior and exterior fluid.As a result, the interior fluid may have an influence on the particle motion [38].Thus, Aidun et al. [20] proposed a non-shell model by excluding the interior fluid, and the momentum exchange is only considered for exterior fluid.As Aidun suggested, for non-shell model a named "impulse force" should be added to the conventional momentum exchange method.This "impulse force" is only calculated when the particle moves to cover or uncover a lattice node, which serves as a correction for the hydrodynamic force and the necessity of this correction for non-shell model has been confirmed by [44].Thus in our simulation the non-shell model is applied and the "impulse force" is used as a correction.When the particle moves across a fluid node the "impulse force" is written as [20]: where CN and UN are respectively fluid nodes being covered and uncovered during a time step.

Benchmarks
To validate the algorithm and the code, two benchmark cases are considered: single particle sedimentation, and two-particle sedimentation.For single particle sedimentation, we replicate the results computed by Feng et al. [45] using the finite element method (FEM), which is often used to validate the fluid-solid coupling methods [20,46,47].For two-particle sedimentation, our results capture the important phenomena named "drafting, kissing and tumbling" or DKT motion [46,47], and agree well with their results quantitatively.For the consideration of stability and accuracy, the dimensionless relaxation time (τ) in LBM is taken as 1, and a relatively high grid resolution (each particle consisting of 64 LBM grids) is applied to reduce the influence of no smoothness brought by midway bounce-back model.In a 2D system, a high grid resolution can be afforded, while it may be computationally expensive for 3D.As an alternative, some smoothing models can be used such as the immersed moving boundary method [41].

Single Particle Sedimentation
Single particle sedimentation geometry is shown in Figure 2a, where d is the particle diameter, and l " 1.5d " 4 cm is the width of the narrow channel.Gravity is along the X axis in the positive direction.The density and kinematic viscosity of fluid is 1000 kg/m 3 and 10 ´6 m 2 /s, respectively.The particle with density of 1010 kg/m 3 is initially placed off the centerline of channel, and then released with a zero velocity.Because the particle is heavier than the fluid, it will settle down along the direction of gravity, and oscillate around the centerline owing to the influence of channel walls as in Figure 2b.Eventually a steady state can be achieved, where the particle reaches a terminal velocity with no lateral motion, and the terminal particle Reynolds number can be calculated: where u is the terminal velocity of the particle, d is the particle diameter, and ν is the kinematic viscosity.
In the case of Re = 6.2, Figure 2b shows the comparison of particle settling trajectories obtained using LBM-DEM and by Feng et al. using FEM [45], and they agree well with each other.

Benchmarks
To validate the algorithm and the code, two benchmark cases are considered: single particle sedimentation, and two-particle sedimentation.For single particle sedimentation, we replicate the results computed by Feng et al. [45] using the finite element method (FEM), which is often used to validate the fluid-solid coupling methods [20,46,47].For two-particle sedimentation, our results capture the important phenomena named "drafting, kissing and tumbling" or DKT motion [46,47], and agree well with their results quantitatively.For the consideration of stability and accuracy, the dimensionless relaxation time (τ) in LBM is taken as 1, and a relatively high grid resolution (each particle consisting of 64 LBM grids) is applied to reduce the influence of no smoothness brought by midway bounce-back model.In a 2D system, a high grid resolution can be afforded, while it may be computationally expensive for 3D.As an alternative, some smoothing models can be used such as the immersed moving boundary method [41].

Single Particle Sedimentation
Single particle sedimentation geometry is shown in Figure 2a, where d is the particle diameter, and is the width of the narrow channel.Gravity is along the X axis in the positive direction.The density and kinematic viscosity of fluid is 1000 kg/m 3 and 10 −6 m 2 /s, respectively.The particle with density of 1010 kg/m 3 is initially placed off the centerline of channel, and then released with a zero velocity.Because the particle is heavier than the fluid, it will settle down along the direction of gravity, and oscillate around the centerline owing to the influence of channel walls as in Figure 2b.Eventually a steady state can be achieved, where the particle reaches a terminal velocity with no lateral motion, and the terminal particle Reynolds number can be calculated: where u is the terminal velocity of the particle, d is the particle diameter, and  is the kinematic viscosity.In the case of Re = 6.2, Figure 2b shows the comparison of particle settling trajectories obtained using LBM-DEM and by Feng et al. using FEM [45], and they agree well with each other.

Two-Particle Sedimentation
For further validating the program, two-particle sedimentation is simulated.The channel geometry is shown in Figure 3a, where the width and length of channel are 2 cm and 8 cm, respectively.The fluid has the properties of water with kinematic viscosity 10 −6 m 2 /s and density

Two-Particle Sedimentation
For further validating the program, two-particle sedimentation is simulated.The channel geometry is shown in Figure 3a, where the width and length of channel are 2 cm and 8 cm, respectively.The fluid has the properties of water with kinematic viscosity 10 ´6 m 2 /s and density 1000 kg/m 3 .The diameter of particle is 0.2 cm, and particles density is 1010 kg/m 3 .The first particle is at 0.999 cm and 0.8 cm in Y and X direction, respectively.The second particle is at 1.0 cm and 1.2 cm in Y and X direction, respectively.For comparison, all these parameters are the same as those in [46,47].Initially, the two particles are released with a zero velocity, and the upper one will settle in the wake of the other, which leads to the rearrangement mechanism called "drafting, kissing, and tumbling" or DKT motion [48].Drafting means the particle in the wake will accelerate with an increasing acceleration owing to the low pressure in the wake.Then two particles contact with each other, which is called kissing.Due to the instability of contacting particles aligned in the direction parallel to the motion, they tend to "tumble" to another position.As a result, the center line of particles rotates an angle, and eventually the two particles depart from each other due to the unsymmetrical wake [48].
1000 kg/m 3 .The diameter of particle is 0.2 cm, and particles density is 1010 kg/m 3 .The first particle is at 0.999 cm and 0.8 cm in Y and X direction, respectively.The second particle is at 1.0 cm and 1.2 cm in Y and X direction, respectively.For comparison, all these parameters are the same as those in [46,47].Initially, the two particles are released with a zero velocity, and the upper one will settle in the wake of the other, which leads to the rearrangement mechanism called "drafting, kissing, and tumbling" or DKT motion [48].Drafting means the particle in the wake will accelerate with an increasing acceleration owing to the low pressure in the wake.Then two particles contact with each other, which is called kissing.Due to the instability of contacting particles aligned in the direction parallel to the motion, they tend to "tumble" to another position.As a result, the center line of particles rotates an angle, and eventually the two particles depart from each other due to the unsymmetrical wake [48].
As show in Figure 3b, our result captures the rearrangement mechanism of "drafting, kissing and tumbling".Furthermore, the quantitative comparison of particle trajectories in Y direction is presented in Figure 4, which shows that before particles collision our result agree well with that in [46,47].After kissing stage, exact agreement may not be expected due to different discrete methods are used [47].As show in Figure 3b, our result captures the rearrangement mechanism of "drafting, kissing and tumbling".Furthermore, the quantitative comparison of particle trajectories in Y direction is presented in Figure 4, which shows that before particles collision our result agree well with that in [46,47].After kissing stage, exact agreement may not be expected due to different discrete methods are used [47].
Computation 2016, 4, 15 7 of 18 1000 kg/m 3 .The diameter of particle is 0.2 cm, and particles density is 1010 kg/m 3 .The first particle is at 0.999 cm and 0.8 cm in Y and X direction, respectively.The second particle is at 1.0 cm and 1.2 cm in Y and X direction, respectively.For comparison, all these parameters are the same as those in [46,47].Initially, the two particles are released with a zero velocity, and the upper one will settle in the wake of the other, which leads to the rearrangement mechanism called "drafting, kissing, and tumbling" or DKT motion [48].Drafting means the particle in the wake will accelerate with an increasing acceleration owing to the low pressure in the wake.Then two particles contact with each other, which is called kissing.Due to the instability of contacting particles aligned in the direction parallel to the motion, they tend to "tumble" to another position.As a result, the center line of particles rotates an angle, and eventually the two particles depart from each other due to the unsymmetrical wake [48].
As show in Figure 3b, our result captures the rearrangement mechanism of "drafting, kissing and tumbling".Furthermore, the quantitative comparison of particle trajectories in Y direction is presented in Figure 4, which shows that before particles collision our result agree well with that in [46,47].After kissing stage, exact agreement may not be expected due to different discrete methods are used [47].

Numerical Results and Discussion
This section is the main part of this paper.Firstly, the biaxial compression test is simulated to validate the contact bond model and show the influence of bonding strength on macroscopic mechanical behavior of cemented granular media.Then the sand production of cemented granular media is simulated, and the influences of bonding strength and flow rate on sand production are discussed.
The main aim of this work is to investigate the importance of often ignored bonding effect on sand production in previous work.For simplicity and clarity, we treat the porous structure as a hexagonal packing with monodisperse spheres to highlight the bonding effect between particles without the influence of packing way and particle size distribution.In future work, we will refine the porous structure with more realistic sand grains, including the effect of the particle size distribution [25] and non-spherical particles such as polyhedron [49].

Biaxial Compression Simulation
Although the bonding effect on macroscopic mechanical behavior of cemented granular media is complex, some general features can still be observed by compression test.These effects include that in cemented materials the macroscopic strength is enhanced and shear band is more likely to be observed [50].In order to validate the contact bond model and show the influence of bonding strength on macroscopic mechanical behavior, biaxial compression is simulated to capture these effects.The main steps in biaxial compression simulation are sample generation, installing specified isotropic stress and compression.
Firstly, the rectangular DEM samples are bounded by four frictionless rigid walls as in Figure 5a.For capturing the bonding effect in cemented granular media, the contact bond model is added to the particles contacting with each other.In the simulation, three samples are considered, which have the same micro parameters except bonding strength.Due to the limitations of pressure differences that can be applied when using LBM, the samples used in the biaxial compression simulation and sand production simulation have a higher deformability and smaller strength than the real materials, which is a common consideration on simulating sand production by LBM-DEM model [23].Based on this principal, the main DEM parameters used in biaxial compression and sand production simulation are chosen empirically and summarized in Table 1. Figure 6 shows the deviatoric stress-strain responses of cemented granular media with different bonding strength.It is observed that the synthetic cemented granular media presents a brittle behavior, and the higher the bonding strength is, the more the macroscopic strength is enhanced, which is compatible with the experimental finding in [50].In addition, the shear band is captured in the biaxial compression simulation as in Figure 7. Therefore our model successfully captures the main features of cemented granular media, and it can be used for a further study.
After sample generation, the specimen is isotropically compacted by four walls to reach the isotropic state of stress.Then the biaxial compression is carried out, where the top and bottom walls are moved inward as strain loading condition, and the left and right walls are allowed to displace to retain the confining stress as 111 Pa.In the simulation, the velocity of loading wall is set to a small enough value to ensure a quasi-static process.
Computation 2016, 4, 15 9 of 18 captured in the biaxial compression simulation as in Figure 7. Therefore our model successfully captures the main features of cemented granular media, and it can be used for a further study.In cemented granular media, bond is a primary property determining its macroscopic behavior, and the failure of granular media is generally owing to the bond breakage [24].Thus the number of cumulative bond breakage is monitored during the biaxial compression simulation to show the relationship between bond breakage and macroscopic mechanical behavior.Figure 8 presents the number of accumulative bond during the biaxial compression simulation.It is observed that when the failure of sample occurs, the number of bond breakage increases greatly, just as in [51].Thus, in sand production simulation, the number of accumulative bond breakage will be monitored as an indicator of the damage in the sample, and the larger the number is, the closer the sample is to macroscopic failure.In cemented granular media, bond is a primary property determining its macroscopic behavior, and the failure of granular media is generally owing to the bond breakage [24].Thus the number of cumulative bond breakage is monitored during the biaxial compression simulation to show the relationship between bond breakage and macroscopic mechanical behavior.Figure 8 presents the number of accumulative bond during the biaxial compression simulation.It is observed that when the failure of sample occurs, the number of cumulative bond breakage increases greatly, just as in [51].Thus, in sand production simulation, the number of accumulative bond breakage will be monitored as an indicator of the damage in the sample, and the larger the number is, the closer the sample is to macroscopic failure.

Sand Production Simulation
Sand production is a phenomenon occurring during the extrusion of hydrocarbons from reservoirs, which makes a lot of trouble on the oil or gas production.Just as discussed above, the behavior of natural sands evidently is distinct from clean sands usually employed in laboratory owing to the influence of bond between sand grains [25].Therefore, sand production of cemented granular media is simulated in this part by LBM-DEM with contact bond model, and the influences of bonding strength and flow rate on sand production are discussed.

Physical Model
A cross section of a single perforation through a weakly cemented reservoir is shown in Figure 9a.Owing to the high computational costs, it is difficult to consider the whole cross section.Thus the right part of the cross section is isolated to consider the flow erosion and stress loading on matrix, just as in Figure 9b.

Sand Production Simulation
Sand production is a phenomenon occurring during the extrusion of hydrocarbons from reservoirs, which makes a lot of trouble on the oil or gas production.Just as discussed above, the behavior of natural sands evidently is distinct from clean sands usually employed in laboratory owing to the influence of bond between sand grains [25].Therefore, sand production of cemented granular media is simulated in this part by LBM-DEM with contact bond model, and the influences of bonding strength and flow rate on sand production are discussed.

Physical Model
A cross section of a single perforation through a weakly cemented reservoir is shown in Figure 9a.Owing to the high computational costs, it is difficult to consider the whole cross section.Thus the right part of the cross section is isolated to consider the flow erosion and stress loading on matrix, just as in Figure 9b.The sample is generated by hexagonal packing as in biaxial compression simulation, and for ensuring fluid connectivity in two dimensions, the hydraulic radius is used, which is a common method to deal with flow in two dimensions [7].In three-dimensionally hexagonal packing of monodisperse sphere, the pore space is one third of the particle radius.In order to ensure the pore space in 2D close to that value, the hydraulic radius of the solid particle is take as 85% the mechanical radius.The cavity in the left part of cemented sample represents the tunnel perforation, which is 1.45 cm, and the rectangular container is created to collect the sand grains eroded by flow.In the following simulation, the dimensionless relaxation time (τ) is equal to 1, and there are 1396 lattices in Y direction leading to the grid size at 1.116 × 10 −5 m.From Equation (2) the time step can be The sample is generated by hexagonal packing as in biaxial compression simulation, and for ensuring fluid connectivity in two dimensions, the hydraulic radius is used, which is a common method to deal with flow in two dimensions [7].In three-dimensionally hexagonal packing of monodisperse sphere, the pore space is one third of the particle radius.In order to ensure the pore space in 2D close to that value, the hydraulic radius of the solid particle is take as 85% the mechanical radius.The cavity in the left part of cemented sample represents the tunnel perforation, which is 1.45 cm, and the rectangular container is created to collect the sand grains eroded by flow.In the following simulation, the dimensionless relaxation time (τ) is equal to 1, and there are 1396 lattices in Y direction leading to the grid size at 1.116 ˆ10 ´5 m.From Equation (2) the time step can be calculated (dt = 1.6 ˆ10 ´4 s) which is a safe value for a stable solution in DEM integration, because it is smaller than the DEM critical time step limit (9 ˆ10 ´4 s) calculated by: where m is the mass of particle, K n is the normal spring stiffness.Moreover, under this grid system the largest Mach number in sand production simulation is O(10 ´3) which is much smaller than 0.1, so the deleterious compressibility effect can be annihilated [15].In order to couple the two mechanisms in sand production: (i) mechanical instabilities and (ii) hydro-mechanical instabilities [6], an external stress is loaded in Y direction and a fluid flow is applied in X direction.Thus the top and bottom wall are moved inwards as a stress loading in Y direction with the left and right wall fixed.During the simulation the strain in y direction is monitored to reflect loading condition.Meanwhile, the fluid with the property of water is set a pressure drop in X direction to capture the flow erosion, and in Y direction the top and bottom wall are treated as no-slip boundary condition.

Damage Evolution
Before quantitative discussion, the phenomenological observation of damage evolution in cemented granular media is presented firstly to understand the process of sand production, which is important to study the instability mechanism on sand production [5].Although visualization is quite cumbersome in experiment, it can be easily achieved using simulations (Figure 10).where m is the mass of particle, Kn is the normal spring stiffness.Moreover, under this grid system the largest Mach number in sand production simulation is O(10 −3 ) which is much smaller than 0.1, so the deleterious compressibility effect can be annihilated [15].
In order to couple the two mechanisms in sand production: (i) mechanical instabilities and (ii) hydro-mechanical instabilities [6], an external stress is loaded in Y direction and a fluid flow is applied in X direction.Thus the top and bottom wall are moved inwards as a stress loading in Y direction with the left and right wall fixed.During the simulation the strain in y direction is monitored to reflect loading condition.Meanwhile, the fluid with the property of water is set a pressure drop in X direction to capture the flow erosion, and in Y direction the top and bottom wall are treated as no-slip boundary condition.

Damage Evolution
Before quantitative discussion, the phenomenological observation of damage evolution in cemented granular media is presented firstly to understand the process of sand production, which is important to study the instability mechanism on sand production [5].Although visualization is quite cumbersome in experiment, it can be easily achieved using simulations (Figure 10). Figure 10 shows the damage evolution of the cemented sample with bonding strength of 0.078 N under the effective stress and flow erosion, where the color of the sand grain represents the number of bond breakage, and the warmer the color is, the more bonds on the sand grain are broken.Initially, the sample is intact with no bonding breakage (Figure 10a).With the increase of the strain, the bond near the cavity begins to be broken owing to the concentration of effective stress, and some small cracks forms in this zone (Figure 10b).Then, these small cracks propagate into the sample, which leads to the localized failure near the cavity (Figure 10c).Eventually, the flow erodes the localized failure zone and transfers the loose sand grains into the cavity (Figure 10d).It is observed that the experimental observation "initial sand production is evoked by localized failure" in [5] is examined by our simulation.What is more, although the localized failure is serious near the cavity, the zone far from the cavity is still intact, which means that sand production occurs well before macroscopic failure, just as the experimental observation by X-ray CT photographs in [5]. Figure 10 shows the damage evolution of the cemented sample with bonding strength of 0.078 N under the effective stress and flow erosion, where the color of the sand grain represents the number of bond breakage, and the warmer the color is, the more bonds on the sand grain are broken.Initially, the sample is intact with no bonding breakage (Figure 10a).With the increase of the strain, the bond near the cavity begins to be broken owing to the concentration of effective stress, and some small cracks forms in this zone (Figure 10b).Then, these small cracks propagate into the sample, which leads to the localized failure near the cavity (Figure 10c).Eventually, the flow erodes the localized failure zone and transfers the loose sand grains into the cavity (Figure 10d).It is observed that the experimental observation "initial sand production is evoked by localized failure" in [5] is examined by our simulation.What is more, although the localized failure is serious near the cavity, the zone far from the cavity is still intact, which means that sand production occurs well before macroscopic failure, just as the experimental observation by X-ray CT photographs in [5].

Effect of Bonding Strength
In cemented granular media, bonding strength is a primary property that influences its mechanical behavior.It is necessary to study the effect of bonding strength on sand production, so three cemented samples with different bonding strength are considered, which are 0.0078 N, 0.039 N and 0.078 N, respectively.Just as illustrated above, during sand production simulation the accumulative bond breakage is monitored as an indicator of the damage caused by effective stress and flow erosion.
Figure 11 shows the number of accumulative bond breakage in three samples under the same flow condition (dp/dx = 2.86 Pa/m).First, a staged growth of accumulative bond breakage is observed for all samples, which means a great increase of bond breakage is always followed by a very slow growth.This is owing to the fact that under effective stress and flow erosion, instability and temporary stability take place alternately in cemented sample.Before macroscopic failure, the cemented granular media is in a temporarily stable state until the effective stress increases to a value strong enough to break this temporarily stable state.Because the macroscopic failure has not occurred, another temporarily stable state can be rebuilt.As a result, on the temporarily stable state, a few bonds are broken in the sample, but when the temporarily stable state is broken, a great increase of bond breakage can be observed.These two states occurring alternately lead to the staged growth of accumulative bond breakage (Figure 11).In cemented granular media, bonding strength is a primary property that influences its mechanical behavior.It is necessary to study the effect of bonding strength on sand production, so three cemented samples with different bonding strength are considered, which are 0.0078 N, 0.039 N and 0.078 N, respectively.Just as illustrated above, during sand production simulation the accumulative bond breakage is monitored as an indicator of the damage caused by effective stress and flow erosion.
Figure 11 shows the number of accumulative bond breakage in three samples under the same flow condition (dp/dx = 2.86 Pa/m).First, a staged growth of accumulative bond breakage is observed for all samples, which means a great increase of bond breakage is always followed by a very slow growth.This is owing to the fact that under effective stress and flow erosion, instability and temporary stability take place alternately in cemented sample.Before macroscopic failure, the cemented granular media is in a temporarily stable state until the effective stress increases to a value strong enough to break this temporarily stable state.Because the macroscopic failure has not occurred, another temporarily stable state can be rebuilt.As a result, on the temporarily stable state, a few bonds are broken in the sample, but when the temporarily stable state is broken, a great increase of bond breakage can be observed.These two states occurring alternately lead to the staged growth of accumulative bond breakage (Figure 11).Although the tendencies of three curves are same, some obvious differences are still observed in Figure 11.First, the bond breakage occurs until the strain reaches a critical point, and higher bonding strength sample corresponds to a higher critical strain.This is because, for intact sample, the seepage force is too weak to break the bond, and the initial bond breakage is evoked by the concentration of effective stress near the cavity.Thus only beyond the critical strain, the effective stress is strong enough to break the bond, and the higher the bonding strength is, the higher effective stress is needed to break it.As a result, higher bonding strength leads to a higher critical strain to evoke the initial bond breakage.Another significant difference in Figure 11 is that under the same strain and flow condition, the accumulative bond breakage in low bonding strength sample is much greater than that in high bonding strength sample.Therefore, the cemented material with lower bonding strength is more likely to produce sand grains, which is confirmed in Figure 12. Figure 12 shows the states of two cemented samples with different bonding strength (0.0078 N and 0.078 N) under the same strain (0.024) and flow condition (dp/dx = 2.86 Pa/m).For higher bonding strength sample, as in Figure 12b, the cavity is stable, and no sand grains are produced.However, for lower bonding strength sample, as in Figure 12a, the stability of cavity is broken, and sand production continues up Although the tendencies of three curves are same, some obvious differences are still observed in Figure 11.First, the bond breakage occurs until the strain reaches a critical point, and higher bonding strength sample corresponds to a higher critical strain.This is because, for intact sample, the seepage force is too weak to break the bond, and the initial bond breakage is evoked by the concentration of effective stress near the cavity.Thus only beyond the critical strain, the effective stress is strong enough to break the bond, and the higher the bonding strength is, the higher effective stress is needed to break it.As a result, higher bonding strength leads to a higher critical strain to evoke the initial bond breakage.Another significant difference in Figure 11 is that under the same strain and flow condition, the accumulative bond breakage in low bonding strength sample is much greater than that in high bonding strength sample.Therefore, the cemented material with lower bonding strength is more likely to produce sand grains, which is confirmed in Figure 12. Figure 12 shows the states of two cemented samples with different bonding strength (0.0078 N and 0.078 N) under the same strain (0.024) and flow condition (dp/dx = 2.86 Pa/m).For higher bonding strength sample, as in Figure 12b, the cavity is stable, and no sand grains are produced.However, for lower bonding strength sample, as in Figure 12a, the stability of cavity is broken, and sand production continues up to end of the simulation.Therefore, for samples under the same condition, lower bonding strength leads to more sand production than higher bonding strength.to end of the simulation.Therefore, for samples under the same condition, lower bonding strength leads to more sand production than higher bonding strength.

The Effect of Flow Rate
Flow rate is the other important factor that influences the process of sand production in cemented reservoir, which leads to the hydro-mechanical instability.In the simulation, three kinds of pressure gradient are considered to study the influence of flow rate on sand production, which are 0.286 Pa/m, 0.572 Pa/m and 2.86 Pa/m, respectively.
Figure 13 shows the influence of flow rate on accumulative bond breakage in cemented granular media with bond strength of 0.0078 N and 0.078 N. It is observed that the influence of flow rate becomes significant until the strain reach some critical value, which confirms the conclusion that the role of flow in sand production is to erode the localized failure zone [4,5].Before the critical strain, the cemented granular media is almost intact, and the hydrodynamic force provided by the highest flow rate cannot erode the sample.Therefore, there is little difference between various pressure gradient conditions.After the critical point, the localized failure near the cavity occurs, and the flow erosion starts to take effects, and the coupling effect of mechanical instability and hydro-mechanical instability becomes significant, so the curves of different flow rates begin to deviate from each other.
Comparing Figure 13a,b, in high bonding strength sample (0.078 N), the accumulative bond breakage under different flow conditions is almost same, but in low bonding strength sample (0.0078 N), the gap between different flow conditions is obvious and higher flow rate leads to more bond breakage.Therefore the influence of flow rate on sand production depends on the bonding strength, and it is more significant in lower bonding strength sample than in higher bonding strength sample.The same result was also obtained in [5], that the effect of flow rate on weak rock is in contrast to that on more competent rock.In competent rock, the localized failure zone is hard to be eroded by flow owing to the residual bond with high strength, so the coupling effect between mechanical instability and hydro-mechanical instability is not significant.However, in weak rock, the localized failure zone is easily eroded, which leads to a significant coupling effect between mechanical instability and hydro-chemical instability, so the influence of flow rate on sand production is more obvious in lower bonding strength sample.

The Effect of Flow Rate
Flow rate is the other important factor that influences the process of sand production in cemented reservoir, which leads to the hydro-mechanical instability.In the simulation, three kinds of pressure gradient are considered to study the influence of flow rate on sand production, which are 0.286 Pa/m, 0.572 Pa/m and 2.86 Pa/m, respectively.
Figure 13 shows the influence of flow rate on accumulative bond breakage in cemented granular media with bond strength of 0.0078 N and 0.078 N. It is observed that the influence of flow rate becomes significant until the strain reach some critical value, which confirms the conclusion that the role of flow in sand production is to erode the localized failure zone [4,5].Before the critical strain, the cemented granular media is almost intact, and the hydrodynamic force provided by the highest flow rate cannot erode the sample.Therefore, there is little difference between various pressure gradient conditions.After the critical point, the localized failure near the cavity occurs, and the flow erosion starts to take effects, and the coupling effect of mechanical instability and hydro-mechanical instability becomes significant, so the curves of different flow rates begin to deviate from each other.
Comparing Figure 13a,b, in high bonding strength sample (0.078 N), the accumulative bond breakage under different flow conditions is almost same, but in low bonding strength sample (0.0078 N), the gap between different flow conditions is obvious and higher flow rate leads to more bond breakage.Therefore the influence of flow rate on sand production depends on the bonding strength, and it is more significant in lower bonding strength sample than in higher bonding strength sample.The same result was also obtained in [5], that the effect of flow rate on weak rock is in contrast to that on more competent rock.In competent rock, the localized failure zone is hard to be eroded by flow owing to the residual bond with high strength, so the coupling effect between mechanical instability and hydro-mechanical instability is not significant.However, in weak rock, the localized failure zone is easily eroded, which leads to a significant coupling effect between mechanical instability and hydro-chemical instability, so the influence of flow rate on sand production is more obvious in lower bonding strength sample.

Conclusions
In this study, the sand production in cemented reservoir is simulated at pore scale by LBM-DEM.In order to account for the bonding effect in cemented granular media, a contact bond model is introduced into the LBM-DEM framework, which is validated by benchmark cases.By sand production simulation, the experimental observation that "initial sand production is evoked by localized failure" is confirmed, and it is also found that instability and temporary stability take place alternately in cemented granular media, until the macroscopic failure occurs.Then the influence of bonding strength on sand production is discussed.Under the same condition, the cemented granular media with low bonding strength is easier to produce sand than that with high bonding strength.Finally, the effect of flow rate is invested, and it illustrates that the influence of flow rate on sand production depends on the bonding strength.In higher bonding strength sample, the influence of flow rate on sand production is not significant.However, in low bonding strength sample, the higher the flow rate is, the more severe the erosion in localized failure zone becomes.Therefore, this simulation captures the main feature of sand production in cemented reservoir and discusses the influence of the micro parameters, which helps us understand the underlying physics of this problem at pore scale.Since only 2D simulations are performed in this paper, in the future we will consider 3D simulations and enlarge our computational scale to make it more suitable for actual situation.

Conclusions
In this study, the sand production in cemented reservoir is simulated at pore scale by LBM-DEM.In order to account for the bonding effect in cemented granular media, a contact bond model is introduced into the LBM-DEM framework, which is validated by benchmark cases.By sand production simulation, the experimental observation that "initial sand production is evoked by localized failure" is confirmed, and it is also found that instability and temporary stability take place alternately in cemented granular media, until the macroscopic failure occurs.Then the influence of bonding strength on sand production is discussed.Under the same condition, the cemented granular media with low bonding strength is easier to produce sand than that with high bonding strength.Finally, the effect of flow rate is invested, and it illustrates that the influence of flow rate on sand production depends on the bonding strength.In higher bonding strength sample, the influence of flow rate on sand production is not significant.However, in low bonding strength sample, the higher the flow rate is, the more severe the erosion in localized failure zone becomes.Therefore, this simulation captures the main feature of sand production in cemented reservoir and discusses the influence of the micro parameters, which helps us understand the underlying physics of this problem at pore scale.Since only 2D simulations are performed in this paper, in the future we will consider 3D simulations and enlarge our computational scale to make it more suitable for actual situation.

Figure 1 .
Figure 1.The computation of lattice Boltzmann method (LBM) and discrete element method (DEM) coupling.

Figure 1 .
Figure 1.The computation of lattice Boltzmann method (LBM) and discrete element method (DEM) coupling.

Figure 2 .
Figure 2. The benchmark case of single particle sedimentation.(a) The geometry of the simulation; (b) the comparison of particle trajectories.

Figure 2 .
Figure 2. The benchmark case of single particle sedimentation.(a) The geometry of the simulation; (b) the comparison of particle trajectories.

Figure 3 .
Figure 3.The benchmark case of two particles sedimentation.(a) The geometry of the simulation; (b) the settling trajectory of two-particle sedimentation.

Figure 4 .
Figure 4. Transverse coordinates of the centers of the two particles.

Figure 3 .
Figure 3.The benchmark case of two particles sedimentation.(a) The geometry of the simulation; (b) the settling trajectory of two-particle sedimentation.

Figure 3 .
Figure 3.The benchmark case of two particles sedimentation.(a) The geometry of the simulation; (b) the settling trajectory of two-particle sedimentation.

Figure 4 .
Figure 4. Transverse coordinates of the centers of the two particles.Figure 4. Transverse coordinates of the centers of the two particles.

Figure 4 .
Figure 4. Transverse coordinates of the centers of the two particles.Figure 4. Transverse coordinates of the centers of the two particles.

Figure 5 .
Figure 5.The DEM sample for biaxial compression simulation.(a) The hexagonal packing of the granular media; (b) the boundary condition for biaxial compression simulation.

Figure 6 .
Figure 6.The deviatoric stress-strain responses of cemented granular media in biaxial compression simulation.

Figure 7 .
Figure 7. Shear band after macroscopic failure of cemented granular media in the biaxial compression simulation.

Figure 5 .
Figure 5.The DEM sample for biaxial compression simulation.(a) The hexagonal packing of the granular media; (b) the boundary condition for biaxial compression simulation.

Figure 5 .
Figure 5.The DEM sample for biaxial compression simulation.(a) The hexagonal packing of the granular media; (b) the boundary condition for biaxial compression simulation.

Figure 6 .
Figure 6.The deviatoric stress-strain responses of cemented granular media in biaxial compression simulation.

Figure 7 .
Figure 7. Shear band after macroscopic failure of cemented granular media in the biaxial compression simulation.

Figure 6 .
Figure 6.The deviatoric stress-strain responses of cemented granular media in biaxial compression simulation.

Figure 5 .
Figure 5.The DEM sample for biaxial compression simulation.(a) The hexagonal packing of the granular media; (b) the boundary condition for biaxial compression simulation.

Figure 6 .
Figure 6.The deviatoric stress-strain responses of cemented granular media in biaxial compression simulation.

Figure 7 .
Figure 7. Shear band after macroscopic failure of cemented granular media in the biaxial compression simulation.

Figure 7 .
Figure 7. Shear band after macroscopic failure of cemented granular media in the biaxial compression simulation.

Figure 8 .
Figure 8.The deviatoric stress-strain response and number of accumulative bond breakage during biaxial compression: (a) bonding strength is 0.0078 N; (b) bonding strength is 0.039 N; and (c) bonding strength is 0.078 N.

Figure 8 .
Figure 8.The deviatoric stress-strain response and number of accumulative bond breakage during biaxial compression: (a) bonding strength is 0.0078 N; (b) bonding strength is 0.039 N; and (c) bonding strength is 0.078 N.

Figure 9 .
Figure 9.The cemented granular media used in sand production simulation.(a) The cross section of a single perforation in cemented reservoir; (b) the geometry of sample used in sand production simulation.

Figure 9 .
Figure 9.The cemented granular media used in sand production simulation.(a) The cross section of a single perforation in cemented reservoir; (b) the geometry of sample used in sand production simulation.

Figure 10 .
Figure 10.The process of sand production in cemented granular media.(a) Initially, the sample is intact with no bonding breakage; (b) the bond near the cavity begins to be broken; (c) small cracks propagate into the sample leading to the localized failure near the cavity; (d) the flow erodes the localized failure zone.

Figure 10 .
Figure 10.The process of sand production in cemented granular media.(a) Initially, the sample is intact with no bonding breakage; (b) the bond near the cavity begins to be broken; (c) small cracks propagate into the sample leading to the localized failure near the cavity; (d) the flow erodes the localized failure zone.

Figure 11 .
Figure 11.The number of accumulative bond breakage for three cemented samples under the pressure gradient of 2.86 Pa/m.

Figure 11 .
Figure 11.The number of accumulative bond breakage for three cemented samples under the pressure gradient of 2.86 Pa/m.

Figure 12 .
Figure 12.The state of cemented granular media with different bond strengths under the same effective stress and flow erosion: (a) continuous sand production for low bonding strength 0.0078 N; and (b) the stable cavity for high bonding strength 0.078 N.

Figure 12 .
Figure 12.The state of cemented granular media with different bond strengths under the same effective stress and flow erosion: (a) continuous sand production for low bonding strength 0.0078 N; and (b) the stable cavity for high bonding strength 0.078 N.

Figure 13 .
Figure 13.The number of accumulative bond breakage for synthetic cemented materials: (a) with high bonding strength 0.078 N; and (b) with low bonding strength 0.0078 N under different flow rate.

Figure 13 .
Figure 13.The number of accumulative bond breakage for synthetic cemented materials: (a) with high bonding strength 0.078 N; and (b) with low bonding strength 0.0078 N under different flow rate.

Table 1 .
The main parameters used in discrete element method (DEM).
(dt = 1.6 × 10 −4 s) which is a safe value for a stable solution in DEM integration, because it is smaller than the DEM critical time step limit (9 × 10 −4 s) calculated by: calculated