Three -dimensional CFD -DEM -IBM simulations of sand production in a sandstone formation, using periodic boundaries, were performed using 10000 frictional elastic spheres bonded together and compressed at 1 MPa of overburden pressure. Sand production simulati on geometry and procedure are proposed, in which the cone penetration test (CPT) has been used to investigate a physical perforation penetration of the cemented sandstone material with the real - world grain size distribution from the Ustyurt -Buzachi Sedimen tary Basin . The Immersed Boundary Method (IBM) wa s adapted for the sand production simulation geometry to simulate the fluid flow near the well casing. Oil with low viscosity and density wa s used as an injection fluid (reservoir fluid). Erosion near the perforation tunnel due to the pressure drawdown was examined, where the production of sand particles was initiated during the first flow due to the drag force that lifted the sand particles from the perforation damage z one. At the beginning of the simulati on a sand arch was captured around the perforation tunnel and due to the fluid flow it collapsed and perforation cavity became larger. The amount and mass of produced sand particles were calculated. Keywords: DEM; CFD, IBM, sand production, cone penetratio n test, cemented sand , sandstone reservoir 1. INTRODUCTION The modelling of amount of sand production is still one of major issues in the oil and gas production. This problem is particular prominent in terrigenous/sandstone reservoirs, which is common to observe in majority of the oil and gas reservoir fields. Typical substances of the sandstone reservoir are quarts, mica, rock fragments and different type of mine ral grains with silica, cements or clay (Bjorlykke, 2010). In addition, such reservoir has a high porosity (0.15) and permeability (from 25 to 100 Darcies). When sand particles start to erode from weak sandstone formations for different reasons, the impact could lead to degradation of the equipment and deterioration of the production flow rate. The filling of the well with sand can cease the flow and result in collapse of the casing. The damage to the surface equipment by erosion or plugging can increase th e cost of production. In addition, the management of the eroded sand is vital to prevent the negative impact to t he environment (Carlson et al., 1992). This wide range phenomenon depends on various factors. Veeken et al. (1991) classified parameters and fa ctors influencing sand production into the three categories: formation, completion and production. The formation category is mostly associated with rock and reservoir properties, like strength, permeability and pressure distribution; the type and technique of well completion, such as 2 wellbore orientation, perforation, etc. lie within the completion category; and the production category includes flow rate, drawdown pressure and depletion of the reservoir pressure. Based on these parameters, the types of san d production (transient, continuous and catastrophic) and sanding mechanisms (shear and tensile failures, erosion and transportation of fine particles) will differentiate. Experimental and numerical models of sand production problem are the key approaches in finding a solution for the sand production problem, at the same time with increasing computing power the numerical modelling becomes one of the powerful tools for predicting sand production (Rahmati et al., 2013). This can be developed by incorporating the physics of fluid and solid interactions at the micro scale to improve modelling of the detachment of particles from the intact sandstone and their transport with the fluid flow. It requires the coupled modelling technique of Discrete -Element -Method (DE M) and Computational -Fluid - Dynamics (CFD). Cundall and Strack (1979) introduced the Discrete Element Method (DEM) which is a powerful tool to model sanding mechanism, where each sand particle is considered as an individual entity. This method can simulate disaggregation of particles from the rock matrix, i.e. degradation and dissemination rock due to the erosion to small sand particles. Using this discontinues approach we can calculate the mass and volume of produced sand and study its mechanism. In the las t twenty years, several researchers numerically studied the sand production phenomena using the CFD -DEM approach since it was first introduced by Tsuji et al. (1993). O'Connor et al. (1997) and Jensen and Preece (2000 ) were implemented the 2D continuity equation and Darcy’s law using the FEM into the 2D DEM code (MIMES) developed at the Massachusetts Institute of Technology and introduced the application of 2D CFD -DEM modelling to simulate the sanding problem in oil recovery, in particular to investigate th e role of flow rate in the matrix collapse and the magnitude of bond strength, respectively. O'Connor et al. (1997) used irregular particle shapes as sand gr ains and Jensen and Preece (2000 ) used n -sided polygon shapes; and in both studies the particles we re bonded together by cohesive strength (to represent cementation of sand particles), where only the tensile bond failure mode was considered. Jensen and Preece (2000 ) reported that as the bonds cohesion decreased, the number of particles breaking free fro m the matrix increased. Li and Holt (2002) implemented the 2D Darcy’s Law into the 2D DEM code (PFC2D), where fluid flow networks and simulating flow along the flow paths connecting the voids through the pipes. Further, Li et al. (2006) have modelled the hollow cylinder tests with fluid flow to study the sanding mechanism, where the 2D Darcy’s law using the FEM with a superimposed mesh implemented to the 2D code (PFC2D). Circular disks (scaled -up) were used as a particle shape and inserted bon ds between them, where bonds are a simple beam and have normal and shear stiffness and strength. In their simulations, Li et al. (2006) obtained three typical failure patterns, which is analogous to those observed in laboratory experiments. In case, if the material is prone to localized compressive failure due to grain crushing, the observed failure pattern is a slit -like breakout. The uniform failure of material around the borehole was observed, when the tensile strength was low and the material was weak. In the last case, if the material is relatively strong and does not form shear bands, the observed failure pattern is that of dog-eared breakouts. All the coupling techniques mentioned above used Darcy’s Law for the sand production simulation. The Darcy f luid flow is not the best description of the physics where particles are being transported by the fluid, it is insufficient to represent real 3D fluid flow effect during the oil recovery process and it is valid for slow and viscous flow. O'Connor et al. (1 997) and Jensen 3 and Preece (2001) suggested using a much larger number of particles to avoid the influence of boundary effects. The need to model the problem using 3D CFD -DEM (Rahmati et al., 2013) is conditioned by the three -dimensional complex geometry p ore flow networks and intergranular. Cheung (2010) implemented the simplified 1 -D fluid scheme calculating solely the continuity equation into the PFC3D to simulate the impact of fluid flow on sand production. The scaled -up DEM particles generated in a blo ck sample, where particles then bonded by parallel bond and compressed for the perforation simulations, the cylindrical specimen was cut from this block sample and the inner cylindrical wall was simply deleted to mimic the perforation hole. After, this spe cimen was overlain by fluid cells in cylindrical coordinate, where the fluid flow assumed to flow radially from the outer cylinder to inner perforation hole. This simplified 1 -D fluid scheme was unable to simulate the full effect of fluid flow as circumfer ential and vertical flow were considered to be significant, where the drag force was applied in one direction. Only the effect of the fluid flow on the particles was considered, but not the effect of the particles on the fluid flow. Also, the radial flow i s not representative at the perforation tip, where the flow is in multidirectional. Zhou et al. (2011) coupled 3D CFD based on the Navier -Stokes equations into the 3D DEM code (PFC3D) to simulate the sand erosion by fluid flow in a weakly bonded sand asse mbly. A portion of cylindrical domain was used as geometry, where spherical particles (without up - scaling) bonded together, and linear -elastic law was used to calculate the contact forces between particles. They simply removed the inner cylindrical wall, w hich mimic a perforation hole. Authors concluded that the fluid flow rate and bonding properties are significantly affects to sand erosion, erosion rate increased when the sand strength decreased and fluid rate increased. They simulated the simplified cond itions and concluded the positive applicability of microscale modelling to study sand erosion. In order to represent laboratory and field conditions better, model boundary conditions were suggested to be used by the authors. Climent (2016) coupled 3D DEM with 3D CFD based on the Navier -Stokes equations (PFC3D). Author used sandstone analogue prepared by Cheung (2010) and tested it under different fluid flow conditions. She reported that with increasing the fluid flow velocity, the plastic region and sand p roduction rate were also increase. However, this model constrained to consider influence of gradient of pressure on particles and was not able to represent the fluid pressure on each particle. Cui et al. (2016) developed particle -fluid model (PFC3D) using the Darcy’s Law equations and proposed a new methodology for grid block permeability and porosity distribution at in response both to deformation and sanding, also they used a new approach to sand grains seepage forces calculation on. Scaled -up spherical particles (bonded by parallel bonds) were generated in a rectangular DEM box, which was then overlain with fluid cells in cylindrical coordinate system. For the imitation of perforation hole, the cylindrical borehole was created by a slow grain stiffness reduction inside the borehole to zero and grains inside the borehole were removed. This block -shaped sample was tested under different pressure conditions and far -field stresses. Authors concluded that the size and failure mode is affected by fluid flow th e breakout zone and sand production as well as the boundary stress. The less sanding was observed at boundary fluid pressure despite higher boundary stress, whereas strong frictional interlocking caused by high tangential stresses around the borehole cause d by high confining stress resulted in moderate sand production. A large amount of sand was produced at high flow rates despite lower far -field stress. 4 All the coupling techniques using 3D DEM -fluid models and sand production simulations mentioned above were performed in PFC3D. All authors, except Zhou et al. (2011) used scaled - up DEM particles and bonded these particles by parallel bonds, and simplifies perforation technique and they suggested to conduct physical perforation tests; to use periodic boundaries to make calibration process more efficient, to use the DEM particles similar to the real grain sizes for the future work. Climent (2016) concluded th at the parallel bond model calibration was still simplistic and a development of the better calibration with more macro parameters was needed and upscaling of model in capturing behaviour of fine scale system was desirable to extend knowledge of the fluid -particle interaction in the multiscale system. This preprint paper presents the sample preparation procedure/method with physical perforation penetration and sand production modelling in periodic cell using the 3D CFD -DEM modelling and complex boundaries using the Immersed Boundary Method (IBM), and all the suggestions made by the above mentioned researchers were taking into account. As the preliminary results, one simulation test case has been presented to show the capability of this method to capture the damage zone (skin) due to the perforation penetration and the sanding mechanism (erosion near the perforation hole) due to the fluid flow. 2. NUMERICAL METHODS AND GOVERNING EQUATIONS Kafui et al. (2002) developed DEM fluidised bed simulation code at Bi rmingham University by coupling the DEM code provided by Cundall and Strack (1979) and the CFD code based on 2D and 3D Navier -Stokes solvers for compressible flows provided by Kuipers et al. (1993). Further, Guo et al. (2012) incorporated the IBM method into this CFD -DEM approach to model gas and particle two -phase flows with complex and moving boundaries. In this study, the 3D CFD -DEM -IBM code was used, where the IBM part was adapted for the geometry of sand production modelling. All the DEM sand particle s are considered to be elastic and frictional, where the particle -LSO (Large -Sized Objects) and particle -particle contact forces without adhesion are calculated using the theories of Hertz (1881) to determine the normal force, and Mindlin (1949) to calcula te the tangential force. A simple 3D bond contact model for cemented sandstone material developed by Rakhimzhanova et al. (2018, 2019 a) was used to calculate the interparticle normal contact force with adhesion, where authors simply modified the previous existing JKR model for auto -adhesive silt sized particles provided by Johnson et al. (1971) and Johnson (1976). Further, the 2D and 3D Navier -Stokes solvers for compressible flow were modified to an incompressible flow solver since the CFD -DEM modelling of sand production is not a large -scale simulation that can capture the density changes caused by pressure or temperature changes as in r eal field conditions. The flow can be approximated as incompressible if flow -induced pressure changes do not cause significant density changes. 2.1. FLUID -PARTICLE INTERACTION FORCE We consider the single -phase fluid and large numbers of fine solid particles form numerical modelling. The total force 𝐹𝑖 acting on particle i incorporated into the CFD -DEM -IBM approach has several components: fci is the sum of the solid particle -particle contact force and the particle - 5 LSO contact force, f fpi is the fluid -particl e interaction force, a gravitational force 𝑚i𝑔. The translational and rotational motions of each particle are governed by Newton’s second law: 𝑚i𝑑𝜈i 𝑑𝑡=Fi=fci+ffpi+mig (1) 𝐼i𝑑𝜔i 𝑑𝑡=Ti (2) where 𝐼i, 𝜈iare the moment um inertia, angular and linear velocities of the particle i, respectively; and T i is the torque arising from the tangential components of the contact force. Following Anderson and Jackson (1967), the fluid -particle interaction force, f fpi, can be expressed a s: ffpi=−𝜈pi∇𝑝+𝜈pi∇∙τf+𝜀fdi (3) where 𝜈pi is the volume of particles, 𝑝 is the local fluid pressure, 𝜀 is the local void fraction, fdi is the drag force, τf is the local viscous stress tensor. Di Felice’s (1994) proposed an empirical correlation to calculate the drag force fdi on particle i 𝑓𝑑𝑖=1 2𝐶𝐷𝑖𝜌𝑓𝑗𝜋𝑑𝑝𝑖2 4𝜀𝑗2|𝑢𝑓𝑗−𝑣𝑖|(𝑢𝑓𝑗−𝑣𝑖𝜒)𝜀−𝜒 (5) where 𝜌fj, ufj and 𝜀j are the fluid density, fluid velocity and void fraction of the fluid cell 𝑗, where particle 𝑖 of diameter 𝑑pi resides; 𝜀−𝜒 is a correction for the presence of other particles and the correlation equation for 𝜒 obtained by Di Felice (1994) which tak es account of the variation of this exponent in the intermediate flow regime as well as the near constant values in the low and high Reynolds number flow regimes: 𝜒=3.7−.065 𝑒𝑥𝑝 [−(1.5−𝑙𝑜𝑔 10𝑅𝑒𝑝𝑖)2 2] (7) 𝐶Di is the fluid drag coefficient for a single unhindered particle: 𝐶Di=[0.63+4.8 √Repi]2 (6) The particle Reynolds number 𝑅𝑒𝑝𝑖, which is based on the superficial slip velocity between fluid and particle, is Repi=𝜌f𝑑pi𝜀j|u−v| 𝜇s (8) We expect that the viscous stress tensor depends only on the fluid motion (Bird et al., 1960) τf=[(𝜇b−2 3𝜇s)∇∙uf]𝛿+𝜇s[(∇uf)+(∇uf)T] (4) where 𝜇b and 𝜇s are the bulk and shear viscosities, uf is velocity of the fluid, respectively. 6 In this study, the incompressible flow was considered; therefor e the fluid 𝜌f is always a constant value. 2.2. PARTICLE -FLUID INTERACTION FORCE The continuity and momentum equations for the hydrodynamics of incompressible fluid in a granular system are 𝑑(𝜌f𝜀) 𝑑𝑡+∇∙𝜌f(𝜀uf)=0 (9) 𝑑(𝜌f𝜀uf) 𝑑𝑡+∇∙𝜌f(𝜀𝑒ufuf)=−∇𝑝+∇∙𝜏f−Ffp+𝜀𝜌fg (10) where 𝐹𝑓𝑝is the fluid -particle interaction force per unit bed volume and it is obtained as Ffp=1 𝑉cell∑ ffpi=−(1−𝜀)𝑛c 𝑖=1∇𝑝+(1−𝜀)∇∙𝜏f+1 𝑉cell∑ 𝜀ffpi𝑛c 𝑖=1 (11) where 𝑛c is the number of particles in a fluid cell of volume 𝑉cell. The momentum equation can be expressed by substituting Eq.11 into Eq.10. 𝑑(𝜌f𝜀uf) 𝑑𝑡+∇∙(𝜌f𝜀𝑒ufuf)=−𝜀∇𝑝+𝜀∇∙𝜏f−Ffp∗+𝜀𝜌fg (12) where Ffp∗ expressed as Ffp∗=1 𝑉cell∑ 𝜀ffpi𝑛c 𝑖=1 (13) 2.3. PARTICLE -FLUID INTERACTION WITH THE IMMERSED BOUNDARY METHOD Peskin (1972) introduced the Immersed Boundary Method (IBM) to simulate arbitrarily complex geometries with sophisticated and mobile boundaries, which do not conform to the CFD -DEM Eulerian rectangular grids. Guo et al. (2013) incorporated the IBM method into an existing CFD -DEM code to model gas -particle two -phase flows with complex and moving boundaries based on the previous IBM development for modelling of particulate systems consisting of gas and solid particles proposed by Kajishima et al. (2001). This method is used when the multiphase system involves a fluid, many fine solid particles and one or more LSOs and it deals with the interaction between fluid flow and LSO. The volume of a single fluid cell 𝑉cell, can be expressed as the sum of the volume of the fine particles, 𝑉p, th volume of the LSO, 𝑉obj, and the volume of the fluid in this cell, 𝑉f 𝑉cell=𝑉p+𝑉obj+𝑉f (13) The volume fraction of the fluid cell excluding the fine particles is defined as the effective void fraction: 7 𝜀𝑒=1−𝑉p 𝑉cell=𝑉obj+𝑉f 𝑉cell (14) and the volume fraction of the LSO in a fluid cell is defined as 𝛼=𝑉obj 𝑉cell (15) Yuki et al (2007) proposed a simple algorithm to evaluate the value of 𝛼 at the collocation point for each cell 𝛼=0.5[1−tanh (∆ 𝜎𝜆𝑑𝑐𝑒𝑙𝑙)] (16) 𝜆=|𝑛𝑥|+|𝑛𝑦|+|𝑛𝑧| (17) 𝜎=0.05(1−𝜆2)+0.3 (18) where 𝑑𝑐𝑒𝑙𝑙 is the cell size, 𝑛=(𝑛𝑥,𝑛𝑦,𝑛𝑧) is a normal outward unit vector at a surface element and 𝜎 is a signed distance from the cell centre to the surface element. Kajishima et al. (2001) established a unified velocity, 𝑢, to describe th e fluid -LSO system by volume averaging the local fluid velocity, uf, and the local solid object velocity, Uobj, for each fluid cell: u=𝑉f 𝑉obj+𝑉fuf+𝑉obj 𝑉obj+𝑉fUobj=𝜀𝑒−𝛼 𝜀𝑒uf+𝛼 𝜀𝑒Uobj (19) where the local LSO velocity, Uobj, is computed by taking average of the local velocities at solid points over the volume of the LSO falling in a fluid cell, 𝑉obj: Uobj=1 𝑉obj∫ (vLSO+𝜔LSO+𝑟)𝑑𝑉=vLSO+1 𝑉obj∫ 𝜔LSO+𝑟𝑑𝑉𝑉obj 𝑉obj (20) in which, vLSO and 𝜔LSO are the translational and rotati onal velocities of the LSO and 𝑟 is the relative position from the centre of mass of the object to a local point on the object. The continuity and momentum equations for the hydrodynamics of incompressible fluid in a granular system and in presence of the LSO are (based on the governed equations for compressible flow by Guo et al., 2012) 𝑑(𝜀𝑒) 𝑑𝑡+∇∙(𝜀𝑒u)=0 (21) 𝜌f𝑑(𝜀𝑒u) 𝑑𝑡+𝜌f ∇∙(𝜀𝑒uu)=−𝜀𝑒∇𝑝+𝜀𝑒∇∙𝜏−Ffp∗+𝜌f𝜀𝑒g+f′ (22) where, the fluid -LSO is treated as a single c ontinuum with a unified velocity field, u , because the no -slip and impermeable conditions at fluid -solid interface with uf=Uobj was considered; and the term f′ is the virtual body force used to correct the velocity field, u, at the fluid -solid interface a nd inside the LSO. Subsequently, the viscous stress tensor 𝜏 is given by Eq.4 is used 8 with u replacing uf; the effective void fraction, 𝜀𝑒, is used in both governing equations to take into account the effect of the presence of fine particles. However, t he real void fraction is used in the Eqs. 5 and 11 to calculate the drag force fdi and the fluid -particle interaction force Ffp∗, respectively and obtained from 𝜀=1−𝑉p 𝑉cell−𝑉obj=𝜀𝑒−𝛼 1−𝛼 (23) The momentum Eq.22 is discretised in time using the first -order finite difference algorithm: 𝜌𝑓(𝜀𝑒𝑢)𝑛+1=𝜌𝑓(𝜀𝑒𝑢)𝑛+∆𝑡∙H𝑛−∆𝑡∙(𝜀𝑒∇𝑝)𝑛+1+∆𝑡∙(f′)𝑛+1 (24) where H=−∇∙𝜌f(𝜀𝑒uu)+𝜀𝑒∇∙𝜏f−Ffp∗+𝜌f𝜀𝑒g (25) and ∆𝑡 is the time step and superscripts implies the number of the time step. Based on the location of the LSO, the main three variables 𝛼, u and f′ of computational cells will be different: 1. For the computational cells outside the LSO, 𝛼=0, u=uf and there is no virtual force (f′)𝑛+1=0, therefore the momentum Eq. 22 reduces to Eq. 12. 2. For the computational cells inside the LSO, the momentum Eq. 22 is used with 𝛼=1, 𝑢𝑛+1=Uobj𝑛+1 and the virtual body force, (f′)𝑛+1, can be derived from Eq. 24: (f′)𝑛+1=𝜌𝑓(𝜀𝑒Uobj)𝑛+1−𝜌𝑓(𝜀𝑒𝑢)𝑛 ∆𝑡−H𝑛+(𝜀𝑒∇𝑝)𝑛+1 (26) 3. For the computational cells at the fluid -solid interface that are partially occupied by the LSO, the momentum Eq. 22 is used with 0<𝛼<1, 𝑢𝑛+1=Uobj𝑛+1 and the virtual body force, (f′)𝑛+1, can be defined by a first -order linear interpolation using 𝛼 defined in Eqs. 16 -18. (f′)𝑛+1=𝛼[𝜌𝑓(𝜀𝑒Uobj)𝑛+1−𝜌𝑓(𝜀𝑒𝑢)𝑛 ∆𝑡−H𝑛+(𝜀𝑒∇𝑝)𝑛+1] (27) Numerical Algorithm For the convenience of the reader we reiterate the previous algorithm ( Guo et al., 2012) with modifications for incompressible flow case. 1. Compute the new velocities and new positions from Eqns. 1 and 2. 2. Compute the new velocity and new position of LSO. 3. For each cell compute 𝜀𝑒, 𝛼 and 𝜀. 9 4. Compute the virtual body force 5. Based on the pressure -correction method we compute 𝑝′. First, 𝑢𝑛+1 is computed from the continuity system (Eq. 28). Second, 𝑢∗ is computed in Eq . 32 with provided 𝑝∗. Next, from Eqns. 33 we obtain 𝑝′. 6. The piecewise constant in cell pressure at (n+1)th time step can be found using 𝑝𝑛+1=𝑝∗+ 𝑝′. 7. Update (f′)𝑛+1 and 𝑢𝑛+1 using 𝑝𝑛+1. 8. Testing the convergence of continuity problem. If not, we repeat procedure from step 4 with 𝑝∗:=𝑝𝑛+1. This iteration continues until the convergence criterion is satisfied. 9. From Eq. 19 we estimate 𝑢𝑓𝑛+1. 2.4.PRESSURE -CORRECTION EQUATIONS A semi -implicit finite difference numerical technique employing a staggered grid is used in solving the incompressible Navier -Stokes equations on an equidistant 3D Cartesian grid. The adapted version of SIMPLE methodology (Semi -Implicit Method for Pressure -Linked Equations) proposed by Patanker (1980) is used as the algorithm. The pressures are defined at the center of each fluid cell, and the fluid velocity components are defined at the cell faces (Figure 1). Figure 1 . 3D control volume (cell ) For the ce ll center O, the continuity equation (21) discretized as (𝜀𝑒)0𝑛+1−(𝜀𝑒)0𝑛+Δ𝑡 Δ𝑥[(𝜀𝑒𝑢)2𝑛+1−(𝜀𝑒𝑢)1𝑛+1]+Δ𝑡 Δ𝑦 [(𝜀𝑒𝑢)4𝑛+1−(𝜀𝑒𝑢)3𝑛+1] +Δ𝑡 Δ𝑧[(𝜀𝑒𝑢)6𝑛+1−(𝜀𝑒𝑢)5𝑛+1]= 0 (28) The discretized momentum equation with virtual force 10 𝜌𝑓(𝜀𝑒𝑢)𝑛+1=(1−𝛼)𝜌𝑓(𝜀𝑒𝑢)𝑛+(1−𝛼)∆𝑡𝐻𝑛−(1−𝛼)∆𝑡(𝜀𝑒∇𝑝)𝑛+1+𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗)𝑛+1 (29) The momentum equation at n+1 at 1 - 6 faces: 𝜌𝑓(𝜀𝑒𝑢)1𝑛+1=(1−𝛼1)𝜌𝑓(𝜀𝑒𝑢)1𝑛+(1−𝛼1)∆𝑡(𝐻𝑥)1𝑛−(1−𝛼1)∆𝑡 ∆𝑥(𝜀𝑒)1𝑛+1(𝑃𝑂𝑛+1−𝑃𝑤𝑛+1)+ 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑥)1𝑛+1 (30.1) 𝜌𝑓(𝜀𝑒𝑢)2𝑛+1=(1−𝛼2)𝜌𝑓(𝜀𝑒𝑢)2𝑛+(1−𝛼2)∆𝑡(𝐻𝑥)2𝑛−(1−𝛼2)∆𝑡 ∆𝑥(𝜀𝑒)2𝑛+1(𝑃𝐸𝑛+1−𝑃𝑂𝑛+1)+ 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑥)2𝑛+1 (30.2) 𝜌𝑓(𝜀𝑒𝑢)3𝑛+1=(1−𝛼3)𝜌𝑓(𝜀𝑒𝑢)3𝑛+(1−𝛼3)∆𝑡(𝐻𝑦)3𝑛−(1−𝛼3)∆𝑡 ∆𝑦(𝜀𝑒)3𝑛+1(𝑃𝑂𝑛+1−𝑃𝑆𝑛+1)+ 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑦) 3𝑛+1 (30.3) 𝜌𝑓(𝜀𝑒𝑢)4𝑛+1=(1−𝛼4)𝜌𝑓(𝜀𝑒𝑢)4𝑛+(1−𝛼4)∆𝑡(𝐻𝑦)4𝑛−(1−𝛼4)∆𝑡 ∆𝑦(𝜀𝑒)4𝑛+1(𝑃𝑁𝑛+1−𝑃𝑂𝑛+1)+ 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑦) 4𝑛+1 (30.4) 𝜌𝑓(𝜀𝑒𝑢)5𝑛+1=(1−𝛼1)𝜌𝑓(𝜀𝑒𝑢)1𝑛+(1−𝛼1)∆𝑡(𝐻𝑧)1𝑛−(1−𝛼1)∆𝑡 ∆𝑧(𝜀𝑒)5𝑛+1(𝑃𝑂𝑛+1−𝑃𝐵𝑛+1)+ 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑧)5𝑛+1 (30.5) 𝜌𝑓(𝜀𝑒𝑢)6𝑛+1=(1−𝛼6)𝜌𝑓(𝜀𝑒𝑢)6𝑛+(1−𝛼6)∆𝑡(𝐻𝑧)6𝑛−(1−𝛼6)∆𝑡 ∆𝑧(𝜀𝑒)6𝑛+1(𝑃𝐹𝑛+1−𝑃𝑂𝑛+1)+ 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑧)6𝑛+1 (30.6) In this method, the pressure at (n+1)th time step can be computed as a sum of the predicted, 𝑝∗, and the corrected 𝑝′: 𝑝𝑛+1=𝑝∗+𝑝′ (31) The predicted mass flow field, 𝜌𝑓𝜀𝑒𝑢, is calculated by the predicted pressure field 𝑝∗ (the momentum equation at 𝑝∗ at 1- 6 faces) 𝜌𝑓(𝜀𝑒𝑢)1∗=(1−𝛼1)𝜌𝑓(𝜀𝑒𝑢)1𝑛+(1−𝛼1)∆𝑡(𝐻𝑥)1𝑛−(1−𝛼1)∆𝑡 ∆𝑥(𝜀𝑒)1𝑛+1(𝑃𝑂∗−𝑃𝑤∗) + 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑥)1𝑛+1 (32.1) 𝜌𝑓(𝜀𝑒𝑢)2∗=(1−𝛼2)𝜌𝑓(𝜀𝑒𝑢)2𝑛+(1−𝛼2)∆𝑡(𝐻𝑥)2𝑛−(1−𝛼2)∆𝑡 ∆𝑥(𝜀𝑒)2𝑛+1(𝑃𝐸∗−𝑃𝑂∗) + 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑥)2𝑛+1 (32.2) 11 𝜌𝑓(𝜀𝑒𝑢)3∗=(1−𝛼3)𝜌𝑓(𝜀𝑒𝑢)3𝑛+(1−𝛼3)∆𝑡(𝐻𝑦)3𝑛−(1−𝛼3)∆𝑡 ∆𝑦(𝜀𝑒)3𝑛+1(𝑃𝑂∗−𝑃𝑆∗) + 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑦) 3𝑛+1 (32.3) 𝜌𝑓(𝜀𝑒𝑢)4∗=(1−𝛼4)𝜌𝑓(𝜀𝑒𝑢)4𝑛+(1−𝛼4)∆𝑡(𝐻𝑦)4𝑛−(1−𝛼4)∆𝑡 ∆𝑦(𝜀𝑒)4𝑛+1(𝑃𝑁∗−𝑃𝑂∗) + 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑦) 4𝑛+1 (32.4) 𝜌𝑓(𝜀𝑒𝑢)5∗=(1−𝛼5)𝜌𝑓(𝜀𝑒𝑢)5𝑛+(1−𝛼5)∆𝑡(𝐻𝑧)5𝑛−(1−𝛼5)∆𝑡 ∆𝑧(𝜀𝑒)5𝑛+1(𝑃𝑂∗−𝑃𝐵∗) + 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑧)5𝑛+1 (32.5) 𝜌𝑓(𝜀𝑒𝑢)6∗=(1−𝛼6)𝜌𝑓(𝜀𝑒𝑢)6𝑛+(1−𝛼6)∆𝑡(𝐻𝑧)6𝑛−(1−𝛼6)∆𝑡 ∆𝑧(𝜀𝑒)6𝑛+1(𝑃𝐹∗−𝑃𝑂∗) + 𝛼𝜌𝑓(𝜀𝑒𝑈𝑜𝑏𝑗𝑧)6𝑛+1 (32.6) The momentum equation at n+1 as functions of momentum at * and p’ (after 32 -30 equations for all faces): 𝜌𝑓(𝜀𝑒𝑢)1𝑛+1=𝜌𝑓(𝜀𝑒𝑢)1∗−(1−𝛼1)∆𝑡 ∆𝑥(𝜀𝑒)1𝑛+1(𝑃𝑂′−𝑃𝑊′) (33.1) 𝜌𝑓(𝜀𝑒𝑢)2𝑛+1=𝜌𝑓(𝜀𝑒𝑢)2∗−(1−𝛼2)∆𝑡 ∆𝑥(𝜀𝑒)2𝑛+1(𝑃𝐸′−𝑃𝑂′) (33.2) 𝜌𝑓(𝜀𝑒𝑢)3𝑛+1=𝜌𝑓(𝜀𝑒𝑢)3∗−(1−𝛼3)∆𝑡 ∆𝑦(𝜀𝑒)3𝑛+1(𝑃𝑂′−𝑃𝑆′) (33.3) 𝜌𝑓(𝜀𝑒𝑢)4𝑛+1=𝜌𝑓(𝜀𝑒𝑢)4∗−(1−𝛼4)∆𝑡 ∆𝑦(𝜀𝑒)4𝑛+1(𝑃𝑁′−𝑃𝑂′) (33.4) 𝜌𝑓(𝜀𝑒𝑢)5𝑛+1=𝜌𝑓(𝜀𝑒𝑢)5∗−(1−𝛼5)∆𝑡 ∆𝑧(𝜀𝑒)5𝑛+1(𝑃𝑂′−𝑃𝐵′) (33.5) 𝜌𝑓(𝜀𝑒𝑢)6𝑛+1=𝜌𝑓(𝜀𝑒𝑢)6∗−(1−𝛼6)∆𝑡 ∆𝑧(𝜀𝑒)6𝑛+1(𝑃𝐹′−𝑃𝑂′) (33.6) For the cell O, the pressure -correction equation can be obtained by substituting Equations A6 into Equation 1 and reorganizing the equation: 𝑎𝑂𝑝𝑂′+𝑎𝑊′𝑝𝑊′+𝑎𝐸′𝑝𝐸′+𝑎𝑆′𝑝𝑆′+𝑎𝑁′𝑝𝑁′+𝑎𝐵′𝑝𝐵′+𝑎𝐹′𝑝𝐹′=𝑏 (34) where 𝑎𝑊′=(1−𝛼1)(∆𝑡 ∆𝑥)2 (𝜀𝑒)1𝑛+1 (35) 𝑎𝐸′=(1−𝛼2)(∆𝑡 ∆𝑥)2 (𝜀𝑒)2𝑛+1 (36) 𝑎𝑆′=(1−𝛼3)(∆𝑡 ∆𝑦)2 (𝜀𝑒)3𝑛+1 (37) 12 𝑎𝑁′=(1−𝛼4)(∆𝑡 ∆𝑦)2 (𝜀𝑒)4𝑛+1 (38) 𝑎𝑆′=(1−𝛼5)(∆𝑡 ∆𝑧)2 (𝜀𝑒)5𝑛+1 (39) 𝑎𝑁′=(1−𝛼6)(∆𝑡 ∆𝑧)2 (𝜀𝑒)6𝑛+1 (40) 𝑎𝑂=−(𝑎𝑊′+𝑎𝐸′+𝑎𝑆′+𝑎𝑁′+𝑎𝐵′+𝑎𝐹′) (41) and 𝑏=𝜌𝑓(𝜀𝑒)0𝑛+1−𝜌𝑓(𝜀𝑒)0𝑛+Δ𝑡 Δ𝑥[𝜌𝑓(𝜀𝑒𝑢)2∗−𝜌𝑓(𝜀𝑒𝑢)1∗]+Δ𝑡 Δ𝑦[𝜌𝑓(𝜀𝑒𝑢)4∗−𝜌𝑓(𝜀𝑒𝑢)3∗]+ Δ𝑡 Δ𝑧[𝜌𝑓(𝜀𝑒𝑢)6∗−𝜌𝑓(𝜀𝑒𝑢)5∗] (42) 3. NUMERICAL SAMPLE PREPARATION METHOD AND SIMULATION PROCEDURE S Prediction of sand production is a complex and difficult task because it covers multidisciplinary domains like geology, drilling, geophysics, geomechanics, fluid mechanics, reservoir engineering, petrol eum engineering, production engineering, etc. Additionally, a hydrocarbon’s field life has to be fully studied and take into account, where it has three major phases: Exploration, Development, and Production. All of each phase has specific objectives and a huge data for sand production problems can be collected. Capturing all the parameters influencing sand production classified by Veeken et al . (1991) and its mechanisms in the numerical models is difficult and the models have many limitations. Therefore, the aim of this study was not to simulate the real downhole conditions but to simplify them, reducing the number of parameters. To achieve this simplif ication and at the same time to mimic the real field conditions, for the well completion category one small perforation hole of the cased horizontal well was chosen as a simulation geometry (Figure 2). All hydrocarbon reservoirs themselves are complicated geological heterogeneous bodies with rock properties (porosity, permeability, saturation, etc.) varying within the reservoir and they are not the homogeneous porous media that are often presented on paper s and used in calculations (Dawe, 2004). Since the a im of this research is dedicated to studying the behavior of sandstone material and its sanding mechanism, and the scale of the area under consideration is small (around one small perforation hole) one lithology in terms of cemented sandstone has been chos en and considered as a homogeneous media . For the fluid characterization part, a single liquid phase - oil with low viscosity and density is considered. For the rock characterization part, the previous research conducted by Rakhimzhanova et al. (2018, 2019 a) is used. Authors developed a simple 3D bond contact model for cemented sandstone material by modifying the previous existing JKR model for auto -adhesive silt size sand particles and reported the results obtained from 3D DEM simulations of triaxial compression tests of cemented sandstone of different density from loose to medium dense. The material properties of the cemented analogue samples from the Ustyurt -Buzachi Sedimentary 13 Basin studied by Shabdirova et al. (2016) were re produced for the numerical samples . Authors did not up -scale the spherical DEM particles sizes, i.e. they were equal to the real grain sizes, and these particles were connected together by interface energy. All simulations were performed in a periodic cell , which implies that the ‘sample size’ is infinite, there are no ‘loading boundaries’ and there is no finite size effect that may affect the results. Authors obtained the Mohr -Coulomb strength criterion parameters c’ - inherent shear strength and φ’ - fricti on angle for both systems were obtained and identified the best match with the experimental results strength parameters. The stress -strain curves obtained for the numerical samples showed good agreement with the experimental results of the cemented sandsto ne in terms of the in terms of the shear strength (Figure 3 ). Consequently, the microscopic response of experimental cemented sandstone material was examined by numerical results. I t should be noted that the experimentally studied cemented analogue samples are not the core samples of overburden formations. The core samples are almost never available for comprehensive laboratory testing (Zoback, 2017). Finally , the sand production simulations were divided into five stages: particle generation, pluvial deposi tion, compression, perforation penetration and sand production. This research extends the recent 3D DEM simulations of cone penetration tests in cemented sandstone performed by Rakhimzhanova et al. (2019b) , where authors used the first four stages for cone penetration simulations (Subsections 3.1 -3.4). Authors conducted series of of numerical CPT tests at different bond strength values and presented the result of numerical CPT tests in cemented sandstone which was found to be in good agreement with the Soil Behaviour Type (SBT) classification system from CPT data. In this study, these penetrated samples are used for the three-dimensional simulations of sand production tests . a) Horizontal well (field scale) b) Perforation hole Figure 2. Well completion of horizontal cased well 14 Figure 3. Comparison of numerical and experiment al stress -strain curves for cemented sandstone samples (Rakhimzhanova et al., 2019) 3.1. PARTICLE GENERATION The same eight different particle sizes from the particle size distribution (PSD) data of the cemented analogue samples from the Ustyurt -Buzachi Sedimentary Basin (Figure 4) measured by Qicpic dynamics image analyzer (Shabdirova et al. 2016) were selected. Within a thin rectangular parallelepiped workspace of dimension 8.4 mm in width, 21 mm in height and 1.4 mm in thickness 10000 particles were randomly generated (Figure 5a): 0.15 mm ( 1179 ), 0.18 mm ( 1750), 0.2 mm ( 1505 ), 0.22 mm ( 1558 ), 0.25 mm ( 1973 ), 0.275 mm ( 1119 ), 0.3 mm ( 610), 0.355 mm ( 306). The workspace was divided into small 720 DEM boxes (box dimension of 0.0007 mm) : 12 in width, 3 0 in height and 2 in thickness. The following material properties were used for all particles: Young’s modulus E=70 GPa, Poisson ratio 𝜈=0.3 and experimental particle density 𝜌=2605 kg/m3. 3.2. PLUVIAL DEPOSITION To form the initial homogeneous bed a vertical gravity field was applied to the system, where all particles were pluvially deposited (Figure 5b). One finite wall was created at the bottom of workspace in order to avoid /prevent the periodic boundary condition at this side . The following material properties of iron were used for all walls: Young’s modulus E= 210 GPa, Poisson ratio 𝜈=0.29 and experimental particle density 𝜌=7 900 kg/m3. 3.3. COMPRESS ION Three top finite walls were cr eated at the top of the deposited specimen and moved down with a wall velocity of 0.1 m/s in order to vertically compress the sample; these three finite walls imitate the well casing (Figure 5 c). When a vertical stress level ( Sv overburden pressure) reached Sv = 500 kPa, t he interparticle friction and wall friction coefficients were set to 𝜇𝑝 = 0.3 and 𝜇𝑤 = 0.3 , respectively; the i dentified bond strength values for cemented sandstone by 15 Rakhimzhanova et al. (2019) Г=20 J/m2 was set in order to make a cementation of medium dense specimen (Figure 4 ), any bond breakage event was started to count from this moment . When the overburden pressure reached Sv = 1 MPa, the top three wall velocities were reset to 0. The Overburden pressure is defined as 𝑆𝑣=∫𝜌(𝑧)𝑔dz≈𝜌𝑟𝑔𝑧𝑧 0 (43) where, 𝜌𝑟 is rock density, 𝑔 is gravitational acceleration, and 𝑧 is depth of interest. Using the Eq.43, the depth of interest was found and it is equal to 39.1 m. Figure 4. Grain size distribution (Rakhimzhanova et al., 2019b) 3.4. PERFORATION PENETRATION After the compression of DEM cemented sandstone, a standard cone penetrometer was first created above the bed and then pushed down at a high constant rate of 0.1 m/s in order to avoid computational cost. The penetr ometer had a “radius” R of 0. 8 mm and an apex angle of 60o. It was defined as 2 rigid frictional walls inclined to X -axis at 60o to simulate the cone penetrometer and two rigid vertical frictional walls to simulate the penetrometer sleeve (Figure 5d). The cone penetrometer was removed from the penetrated specimen once the cone tip had reached the depth y/R=3.5. 16 a) Particle generation b) Pluvial deposition c) Compression d) Cone perforation Figure 5. Perforation penetration simulation stages (Rakhimzhanova et al., 2019b) 3.5. SAND PRODUCTION The following fluid properties were used for sand production simulation : oil density 𝜌𝑜𝑖𝑙 = 850 kg/m3 and shear viscosity ηoil = 0.52·10-3 Pa·s-1 (Agada et al., 2014). The following fluid boundary conditions were set for numerical samples: periodic boundary for left and right sides, impermeable free lip for back and front sides, and prescribed pressure for bottom and top side s. Top solid walls represent the perforated casing of the horizontal well , where the two lowest finite horizontal walls (at the end of compression) are not necessarily located to coincide 17 with any fluid cell boundary. Consequently, they treated as immersed boundaries d uring the sand production simulation with the no -slip and impermeable conditions . The upper two finite horizontal walls and two finite vertical walls that connected these upper and lower finite horizontal walls are coincide with fluid cell boundary. The di fferent pressures were specified for the bottom and top sides of the specimen to mimic the real field pressure drawdown (the difference between the reservoir pressure 𝑃𝑟 and the flowing wellbore pressure 𝑃𝑤) that drives fluids from the reservoir into t he wellbore. The reservoir pressure assumed to be equal to the pore pressure at depth, where the pore pressure is defined as a scalar hydraulic potential acting within an interconnected pore space at depth: 𝑃𝑝≈𝑃𝑟=∫𝜌𝑜𝑖𝑙(𝑧)𝑔dz≈𝜌𝑜𝑖𝑙𝑔𝑧𝑜𝑖𝑙𝑧 0 (44) where, 𝑧𝑜𝑖𝑙 = z = 39.1 m and it is a depth of interest (depth, where oil locates) . The reservoir pressure at 39.1 m was calculated using Eq. 44 and it is equal to 𝑃𝑟 = 326 296 Pa (pressure at the bottom side of specimen). The wellbore pressure assumed to be equal 𝑃𝑤 = 101 325 Pa (pressure at the top side of specimen) . 4. NUMERICAL RESULTS Figure 6 shows the sand production at 0.0295 sec for bond strength Г=20 J/m2, and at this state the quantity of produced sand grains is 266 (Figure 7 a) and the mass is 4.33 mg (Figure 7 b). Figure 6 . Sand production 18 a) Quantity b) Mass Figure 7. Data of produced sand The figure 8 shows the produced sand particle positions at initial (0 sec) and current (0.0295 sec) states, where the data of initial positions (data points located below the IBM walls, near the perforation hole) have been taken just after the perforation stage. It has shown that the all produced sand particles (data points located above the IBM walls) came from perforation damage zone (skin), and at the current state the migration of fine particles from remote zone through channels h ad not detected. 050100150200250300 0 0.005 0.01 0.015 0.02 0.025 0.03Quantity of produced sand grains Time, s0.00.51.01.52.02.53.03.54.04.55.0 0 0.005 0.01 0.015 0.02 0.025 0.03Mass of produced sand grains, mg Time, s 19 Figure 8 . Erosion near the perforation hole 00.00070.00140.00210.00280.00350.00420.00490.00560.00630.0070.00770.00840.00910.00980.01050.01120.01190.01260.01330.0140.01470.01540.01610.01680.01750.01820.01890.01960.02030.021 0 0.0007 0.0014 0.0021 0.0028 0.0035 0.0042 0.0049 0.0056 0.0063 0.007 0.0077 0.0084ini_D1=0.355 mm SP_D1=0.355 mm ini_D2=0.3 mm SP_D2=0.3 mm ini_D3=0.275 mm SP_D3=0.275 mm ini_D4=0.25 mm SP_D4=0.25 mm ini_D5=0.22 mm SP_D5=0.22 mm ini_D6=0.2 mm SP_D6=0.2 mm ini_D7=0.18 mm SP_D7=0.18 mm ini_D8=0.15 mm SP_D8=0.15 mm 20 At the very beginning of sand production simulation (0.00003 sec) a sand arch was captured ( Figure 9 ). Figure 9 . Sand arch at 0.00003 sec (the different colours indicate the cluster size: white particles have 1 particle in one cluster, light cyan 2 and dark cyan 3) Figure 10. Illustration of sand arch (Carlson et al., 1992) Carlson et al (1992) stated that sanding onset is caused by the fluid flow initiation and the drag forces associated with it. The sand grains detach from the matrix and transported to perforation hole. The effect is directly proportional to the fluid viscosity and its velocity and with high-pressure differentials d uring drawdown, which results in stable arch around the entrance to be formed aroun d a perforation cavity (Figure 10). The stability of the arch depends on the constant flow velocity and drawdown. After the change in the flow velocity and the drawdown, the arch collapses and a new one is generated. The figure 11 shows the perforation damage zone in term of compressive (orange) and tensile (mild green) contact forces at 0.0295 sec and it demonstrates that the size of arch becomes larger. Figure 1 2 shoes the fluid velocity vectors when the IBM method is used. The sample confi guration for sand production problem is utilized, where the cemented sandstone material has removed in order to see the fluid behaviour near the IBM walls and perforation hole. 21 Figure 11 . The perforation damage zone at 0.0295 sec Figure 12 . Fluid flow near the IBM walls 5. CONCLUSION This preprint paper is the extension of the work conducted by Rakhimzhanova et al. (2019b) , where authors performed 3D DEM simulations of cone penetration in cemented sandstone material with the realistic grain sizes . One of the previously numerically penetrated sandstone sample with bond strength Г=20 J/m2 is used to perform the numerical simulation of sand production in horizontal cased wells using 3D CFD -DEM -IBM approach , where t he IBM method was optimized for sample geometry. Sanding mechanism due to the pressure drawdown was examined, where the fluid flows (oil ) have been drive n due to the difference between the base and top boundary pressures of the workspace and lifted the eroded sand particles . The results show that all produced sand particles came from the perforation damage zone; and the migration of fine particles had not detected from remote zone through channels. Sand arch was numerically captured, and the perforation cavity changed the size and becomes larger due to the sand arch collapse ; the mass and amount of produced sand particles have been ca lculated. The aim of this paper is to demonstrate the capability of the numerical approach to study the sanding mechanism at macro and micro scales. This study will be continued ; number of improved numerical simulations will be con ducted and full results we ll be published in the near future. ACKNOWLEDGMENTS This research was supported by the Nazarbayev University Faculty Development Competitive Research Grant No 110119FD4502 and by the EU Project GEO -RAMP (H2020 -MSCA -RISE - 2014 GA 645665). The authors gratefully acknowledge Prof. Stefano Utili (Newcastle Unive rsity) and Dr. Helen Cheng (University College London) for their help during the first 22 author’s visits to the UK; Prof. Yu Guo (Zhejiang University) for his help and suggestions in the IBM part ; and Sand production research group at Nazarbayev Unives ity.