Numerical simulation of sanding using a coupled hydro-mechanical sand erosion model

Majid Fetrati, Ali Pak

Department of Civil Engineering, Sharif University of Technology, Tehran, Iran

ABSTRACT Mechanical failure of materials adjacent to the production cavity and material disaggregation caused by fluid drag are considered as the most important parameters that affect sand production.In light of such factors, the coupling of two mechanisms-mechanical instability and hydrodynamic erosion-is indispensable in order to model this phenomenon successfully. This paper examines the applicability of a coupled hydro-mechanical erosion criterion for simulating sand production using the finite element method. The porous medium was considered fully saturated. The onset of sanding and production of sand were predicted by coupling mechanical failure and subsequent erosion of the grain particles utilizing a sanding model. To consider the erosion process, the Papamichos and Stavropoulou (1998)’s sand erosion criterion was incorporated into the finite element code. Arbitrary Lagrangian-Eulerian(ALE) adaptive mesh approach was used to account for large amounts of erosive material loss. Besides, in order to address the problem of severe mesh distortion, the “mesh mapping technique” was employed. Sand production in a horizontal wellbore and in a field case was simulated to demonstrate capabilities of the proposed model. In addition, principal parameters affecting sand production,including in situ stresses, cohesion,perforation orientation,and drawdown were examined.The results indicated the efficiency of the model used in evaluation of sanding in the field. Parametric studies indicated that in situ stresses and formation cohesion could be considered as dominant factors affecting the amount of sand production.

Keywords:Sand production Finite element method Sanding criterion Hydro-mechanical coupling

1. Introduction

Sand production is a phenomenon of migrating considerable detached grain particles together with the formation fluid outflow.Sand production can cause a number of problems such as formation collapsing, clogging up of perforation channels, failure of sand control measures, and collapse of some parts of horizontal boreholes(Willson et al.,2002;Wang and Sharma,2016).Also,damage to well equipment,reduction in well productivity,and unfavorable environmental effects caused by sand production can increase the production costs noticeably (Dusseault and Santarelli, 1989;Rahmati et al., 2013). Robust numerical models capable of simulating the onset of sand production and sanding amount are essential for predicting the time and placement of sanding. These models facilitate selection of the sand control measures and effective utilization of the selected measures.In other words,sand control methods,such as gravel packs and chemical consolidation,can be optimized using these models.

In general, two principal mechanisms are involved in sanding phenomenon: (1) rock failure and localized damage around the cavity because of stress concentration; and (2)hydro-mechanical instabilities caused by internal and surface erosion owing to seepage forces leading to the release and transport of degraded rock (Morita and Boyd,1991; Vardoulakis et al.,1996). The most recent sand production models calculated sand production volume only based on mechanical failure, without capturing the hydro-mechanical erosion effect(Morita et al.,1989,1998;Burton et al.,1998; Vaziri et al., 2002; Shen, 2011). Several models tried to couple both the above-mentioned mechanisms (Vardoulakis et al., 1996; Skjaerstein et al., 1997; Papamichos and Malmanger, 1999; Chin and Ramos, 2002; Wang et al., 2005;Gravanis et al., 2015). However, most of the hydro-mechanical analyses did not consider the continuous change of the boundary condition, i.e. erosion of the wellbore. For addressing this issue,the numerical model with capability to model material erosion from the free surfaces of the wellbore is required. Furthermore,gradual updating of the boundary conditions along with the growth of the erosion process on the new surfaces should be considered.

Sand production phenomenon has been investigated and discussed for decades. Various approaches have been developed to help understanding the mechanism of sanding, e.g. empirical methods, laboratory experiments, analytical analyses, and numerical models. The empirical method is a technique based on field experience that establishes a relationship among sand producing well data and operational parameters (Veeken et al., 1991). Well monitoring and providing sand production field data for this method are time-consuming,and usually the extracted correlations are only applicable to a specific area. Another approach is laboratory experiments which have a long history in prediction of sand production (Isehunwa and Olanrewaju, 2010; Fattahpour et al.,2012). Some of the researchers used a real block of sandstone to simulate this phenomenon (Kooijman and Halleck, 1992), whilst some others conducted their tests on synthetic sandstones(Fattahpour et al., 2012). Although these kinds of experiments are suitable for validating numerical models with reliable sanding data,they are also time-consuming and expensive. To address this problem, some researchers proposed various analytical relationships (Isehunwa and Olanrewaju, 2010; Araujo Guerrero et al.,2014; Gholami et al., 2016; Tehrani et al., 2016). Although analytical methods are fast and straightforward to utilize,they are merely proper for predicting the onset of sand production. Another weakness of this method is that simplified geometrical and boundary conditions are used when predicting sanding phenomenon.

Numerical models are powerful tools for studying various aspects of sand production, and experimental results and/or field data can be utilized to calibrate and validate the numerical results(Rahmati et al., 2013). In general, numerical approaches can be classified into two categories - continuum and discontinuum.Discontinuum approach is an effective way to simulate the separation of solid particles from the rock formation, which is a useful tool to understand the mechanism of sanding (Rahmati et al., 2013). Following this, some researchers coupled discrete element method (DEM) with lattice Boltzmann approach to simulate sanding phenomenon(Ghassemi and Pak,2015;Han and Cundall, 2017). However, there are some limitations in this approach. Not only does the discontinuum approach need considerable computational effort compared to the continuum approach,but it also has major limitations in simulating full-scale problems.Thus,it can be applied merely for small-scale problems.For overcoming the limitation of the DEM, the continuum approach was used for simulation of sanding problem in this study.

Literature review indicates that few numerical investigations could verify their simulation results against experimental and/or field observations. For example, Eshiet and Sheng (2013) used Abaqus in their study,but they did not verify the numerical results against experimental observations or field data.

This paper examined the applicability of a coupled hydromechanical and erosion criterion for simulating sand production based on the finite element method. For simulating erosion, sand erosion criterion suggested by Papamichos and Stavropoulou(1998) was selected and incorporated into the finite element code. The erosion equation considered both plastic strain in the rock and fluid flow velocity, which are effective parameters in sanding phenomenon.

2. Model theory

2.1. Adopted material constitutive model

As mentioned previously, two important mechanisms play the major roles in the sand production process: (1) the mechanical instability of the rock adjacent to the borehole, (2) movement of sand particles due to the drag force of the flowing fluid. Thus, an elastoplastic constitutive model should be first considered when a continuum-based approach is concerned. The Mohr-Coulomb model including a hardening/softening behavior was employed in this investigation. For the general state of stress, the model is conveniently written in terms of the stress invariants as(Menetrey and Willam,1995):

whereRmcspecifies the shape of the yield surface in π plane,cis the material cohesion, φ is the friction angle,qis the Mises equivalent stress,pis the equivalent pressure stress, and θ is the deviatoric polar angle defined as

where σ is the stress tensor,ris the third invariant of deviatoric stress,Sis the deviatoric stress,Iis the identity matrix, and the symbol “:” is called the double contracted product, or double dot product, of two second-order tensors. This is defined, for two second-order tensors,AandB, by

It should be mentioned that the plastic strain is calculated using a non-associated flow rule.The flow potentialG,for yield surface of Mohr-Coulomb model is selected as a hyperbolic function and the smooth elliptic function in the meridional and deviatoric stress plane, respectively Menetrey and Willam (1995):

where

where ε is a parameter for the meridional eccentricity,c0is the initial cohesion, ψ is the dilation angle, anderepresents the deviatoric eccentricity.

It has been shown that for ameliorating the accuracy of sand production modeling, an elasto-perfectly plastic model is not satisfactory; therefore, a hardening/softening behavior of the rock matrix must be considered in numerical simulation (Nouri et al.,2009). Entering rock vicinity of the cavity into a strain-softening regime is a prerequisite for disaggregating of the material. To capture the softening behavior of rock,a strain dependent cohesion is considered (Vermeer and de Borst,1984). According to this equation, increasing the strain leads to cohesion degeneration:

Finally, intensity of the damage of the materials near the wellbore or perforation tunnels is assessed by considering the magnitude of the equivalent plastic strain:

where dεplis the plastic strain increment.

2.2. Sanding criterion

As discussed, the sanding process is related to rock failure and hydro-mechanical instability. In this paper, when considering the latter step, one sanding criterion was used. The sanding erosion shown by Eq. (14) was proposed by Papamichos and Stavropoulou(1998):

where ˙mis the erosion rate of solid mass, C is the transport concentration of the fluidized solids, ρsis the solid density, φ is the porosity of rock,qiis the fluid flux,the symbol‖‖ is the norm of the function that only allocates positive values ofqi, λ is the sand production coeffciient,is the critical plastic shear strain at peak strength, and λ1and λ2are the calibration constants which can be determined via calibrations with sand production tests in the laboratory.In this study,the sand erosion model described above was added into the finite element code by writing a subroutine in FORTRAN named U-mesh-motion. Then, the applicability of the model was investigated through comparison of the numerical results with the recorded data from a laboratory experiment,and also with the field data from a real sanding case.

3. Techniques used in numerical simulations

3.1. Arbitrary Lagrangian-Eulerian adaptive mesh technique

In the conventional Lagrangian analysis,each node is connected to the material particle, and material deformation results in element deformation because the material boundary matches the element boundary. Although following moving surfaces and applying boundary conditions are comfortable in this procedure,the mesh distortion may occur by large deformations. Consequently, the simulation results will no longer be reliable. On the contrary, in an Eulerian analysis, material flows through elements while nodes stay fixed in space, and Eulerian elements may be partially or entirely void.In this approach,mesh distortion will not occur,whereas applying boundary conditions would be challenging due to difficulty in tracking moving surfaces.

The adaptive meshing technique in combination with Lagrangian and Eulerian analyses was used to carry out mesh-refinement operations. This technique is generally known as Arbitrary Lagrangian-Eulerian (ALE) analysis. The mechanism of ALE adaptive meshingmoving independently of mesh from the material-makes it feasible to preserve the quality of mesh throughout the simulation, even when either excessive deformation or loss of materials occurs(Fig.1).

Given that sand production phenomenon leads to significant deformation near the wellbore and perforations, it is imperative that the mesh with high quality is secured. According to the importance of mesh quality, in this study, a high-quality mesh during the sand production procedure was achieved by using Arbitrary Lagrangian-Eulerian (ALE) adaptive meshing technique.

Fig.1. Different meshing techniques.

Fig. 2. Schematic illustration of the model geometry.

3.2. Mesh mapping

Another technique that may be needed for simulation of sanding is “mesh mapping”. As time elapses, more sand produces from the well, and it can cause large deformation near the wellbore during simulation. Consequently, there is a tendency for instability in the model.One solution employed in this study is mesh mapping technique. Mesh mapping technique maps the solution from an old,distorted mesh to a new one in order to the simulation can proceed.Using this technique,the instability can be eliminated by providing a new meshwith betterqualityfor the subsequent stepsof the analysis.

In this study,when distortion happened in some elements such that a proper discretization is not provided anymore for the problem domain, the mesh mapping technique was employed.

4. Numerical simulation of laboratory experiment and real field case

4.1. Sand production in horizontal hole inside an artificial sandstone block

This experiment was conducted on a weak artificial sandstone block (26.25 cm × 26.25 cm × 38 cm) with a 25.4 mm diameter horizontal hole along the longer axis of the block (see Fig. 2)(Kooijman et al.,1996). To overcome the difficulty of dealing with highly variable physical features and mechanical properties of a natural poorly consolidated rock, artificial sandstone with consistent properties was used herein. Sanding phenomenon was simulated by applying far-field stresses using flat jacks on different surfaces of the block with a ratio of vertical to horizontal effective stresses of σv/σh=2 during the test.Also,a differential pressure of 172 kPa between the outer sides of the sample and internal hole allowed fluid to flow from the external circumference into the borehole. The pressure inside the internal hole was constant and equal to atmospheric during the test. The characteristics of the synthetic sandstone block used represent a highly permeable sandstone sample which has been poorly consolidated (Table 1).The purpose of this experiment was to examine the erosion of the wall of the horizontal well, and measure the amount of the producing sand under various load conditions (Kooijman et al.,1996).

Table 1Input parameters (Kooijman et al.,1996).

Fig. 3. Finite element mesh.

Fig. 4. Vertical loading applied at the upper boundary and pore pressure during the experiment (Kooijman et al.,1996).

Fig. 5. Cumulative sand production (Nouri et al., 2007; Pak et al., 2015).

Due to symmetry in two perpendicular directions, only onequarter of the model was examined in the course of numerical modeling. In the analysis, the initial conditions, both stresses and pore pressures, were set to zero. Plane strain condition was assumed for simulation.The finite element mesh used is shown in Fig. 3, and the magnitude of vertical loading applied at the upper boundary and the variation of the flow rate applied to the right and top boundaries, are shown in Fig. 4. Throughout the numerical simulation, the pressure inside the central hole was set to atmospheric pressure following the experimental procedure. As can be seen from the input data,during the experiment,the pore pressure was rapidly reached to 172 kPa at maximum and was decreased to zero at two middle stages.Also,the proportion of vertical effective stress to the horizontal effective stress was equal to 2.In numerical modeling, similar conditions were simulated. Physico-mechanical properties used in the numerical simulations are shown in Table 1.

Two-dimensional (2D) 4-node elements with the primary unknowns of displacements(uandv)and pore pressure(p)were used in the numerical simulation. The ALE technique was adopted as described before. Eq. (14) was exercised for simulating erosion in this experiment.

Fig. 6. Variation of void ratio at different time: (a) 16,200 s; (b) 18,860 s, (c) 20,160 s, (d) 20,500 s, (e) 21,160 s, and (f) 21,700 s.

Fig. 7. The velocity of fluid flow on the borehole face vs. circumferential angle at various time (circumferential angle is measured clockwise around the borehole).

Due to the weak nature of the artificial sandstone used in this experiment (Kooijman et al., 1996), the deformation near the borehole was large. Also, in Fig. 4, both amounts of load and its variations are high at the end of the experiment,and these changes occur in a small amount of time,which makes the elements in the vicinity of the well be distorted excessively during the simulation.Consequently,in addition to ALE technique,map solution has been used for avoiding severe element distortion.

The obtained numerical results in this study are compared with the experimental data and the outcomes of other numerical simulations (Pak et al., 2015; Nouri et al., 2007) that used fast Lagrangian analysis of continua (FLAC) in Fig. 5. The obtained results match well with the experimental data. The amount of sanding is constant between 10,000 s and 15,000 s,and 17,500 s to 19,000 s in the experiment; this pattern also can be seen in the result of this study.Considering that both the stress magnitude and stress rate were increasing near the end of experiment and numerical modeling, it can be seen that a large amount of sand produced in a short period of time after 20,000 s. It should be mentioned that restarting of sand production after the first constant rate in the numerical modeling started sooner than the laboratory experiment between 13,500 s and 17,500 s.This occurrence can be explained by sand arching. It is probable that sand arching occurred in the experiment;however,this phenomenon cannot be seen in finite element modeling.

Fig. 6 shows the variation of void ratio with elapsed time surrounding the borehole.Besides,in this figure,the shape of the crosssection of the hole can be seen at different time.According to these contours, the amount of void ratio is low at the beginning of the sanding. However, the larger amount of sand is produced from the right and left sides of the borehole,higher void ratios can be seen in those areas in the contours.As this experiment has been carried out on weak sandstone, the void ratio has high values. The maximum void ratios for “a”and “f”are 0.5536 and 0.6468,respectively.

Fig. 8. Total sand production (SP) of North Sea sandstone field (Papamichos and Malmanger,1999).

Table 2The features of the three sand production time period(Papamichos and Malmanger,1999).

Fig. 7 depicts the velocity profile away from the surface of the borehole at various time.It can be inferred that the velocity of fluid flow nearby the borehole is almost uniform at the time of initial production(9600 s and 15,500 s).Nevertheless,as the eroded zone unevenly expands nearby the borehole,the fluctuation of fluid flux is escalated. According to Fig. 7, when the circumferential angle is around 70°,the maximum fluid velocity occurs.It can be perceived from erosion location that the maximum erosion occurs in places where the maximum fluid velocity takes place.

4.2. Sanding in a real well at the North Sea

In this section, the applicability of the introduced sand erosion model was examined for prediction of the sanding amount in a sandstone field in the North Sea. Volumetric sanding data were gathered from a wellbore located in North Sea reservoir(Papamichos and Malmanger,1999).Fig.8 illustrates the total sand production recorded for this reservoir. The entire period of time is divided into three segments for a better interpretation.The features of these three segments are shown in Table 2.The cumulative sand production for each time interval is shown in Fig. 9. Sand production data are approximated with a parabolic function of time(Papamichos and Malmanger,1999).

In the course of numerical modeling,a 2D finite element model was utilized to predict the amount of sanding from a cased and perforated well. The total number of perforations is 1400, and the number of perforations for each meter at depth is 20.Assuming that 10 of these perforations are on the right side and 10 on the left side of the well, a plane strain assumption is considered in numerical modeling.The 2D plane strain condition is selected for simulation.This premise has been formerly used by other researchers (e.g.Wang, 2003; Wang et al., 2004, 2005; Li et al., 2018).

Fig. 9. Cumulative sand production related to three sand production (SP) intervals(Papamichos and Malmanger,1999).

Fig.10. Perforated-casing completion model.

It was needed to take into account only one-quarter of the model in the numerical simulation because of the symmetry. The wellbore radius is 0.1 m, and the size of one-quarter of the model were chosen as 5 m×5 m to eliminate the boundary effect on the simulation results. The length and radius of the perforation are 0.286 m and 0.01 m,respectively.The 2D 4-node elements with the primary unknowns of displacements (uandv) and pore pressure(p) were used in the numerical simulation. Smaller finite element mesh was used in the vicinity of the perforation, and it gradually coarsened farther away from the perforation (Fig.10).

The simulation consisted of three principal steps.At first,initial pore pressure was applied to the whole domain, and stress equilibrium was obtained in the formation. There were no wellbore or perforation tunnel in this stage. In the second step, the elements which were used to symbolize both perforation and wellbore were deactivated in the interaction module.Moreover,a distributed load whose value was the same as the value of pore pressure of formation,was applied to the perforation surface.In the third step,in order to achieve the desired drawdown, the pore pressure magnitude was decreased to a specific amount on the perforation freesurface. Consequently, drawdown results in fluid flow, and sand starts to displace into the wellbore. The input parameters for simulating this reservoir are shown in Table 3. Eq. (14) was employed to predict the sanding amount in this simulation.

Table 3Input parameters.

As mentioned before, the total time is divided into three time periods for a better interpretation of the field data. At first, sand production phenomenon has been simulated individually for three time intervals. After that, the simulation is also carried out for the total time without changing the input data.

The features of three time periods are shown in Table 4. The simulation results for three time periods and for the total time can be seen in Figs.11 and 12,respectively.Good agreement is observed between the numerical results and field data. Therefore, the proposed sanding criterion has the ability to predict the sanding amount with acceptable accuracy, for a real well in a sandstone reservoir.

Fig.13 illustrates the longitudinal fluid flux along the perforation for three drawdowns corresponding to three time intervals. It can be concluded that the larger drawdown creates the higher fluid velocity in the perforation.Also,it is evident that oscillation of fluid flux on the tip of the perforation is higher than that in other sections.This can be elucidated by changing the perforation shape due to sand production at the tip of the perforation. Considering that the fluid flows into the perforation from all directions whose shape changes unevenly on the tip,the oscillation is elevated on the tip of the perforation compared to that in the other parts.

5. Parametric study

In this section, some principal parameters affecting sanding phenomenon were studied. For this, the effects of in situ stresses,cohesion, perforation orientation, and drawdown for sandstone field were examined.

Table 4The features of three-time periods used in numerical modeling.

Fig.11. Cumulative sanding amount for each time interval.

Fig.12. Total sanding amount.

Fig.13. Longitudinal fluid flux along the perforation channel.

Fig.14. The quantities of sand produced in different isotropic in situ stresses.

5.1. In situ stresses

In situ stresses in the reservoir can be considered as one of the most important factors affecting sand production. In general, with increasing depth, in situ stresses increase. In order to understand the effect of this parameter on sanding,different in situ stresses are applied to the model with the assumption that the stresses change isotropically, while other factors keep unchanged.

Fig.14 indicates that with increasing in situ stress,more sand is produced and vice versa.Also,it can be seen that when in situ stress equals 40 MPa, no sand is produced in the first 25 h because the sanding criterion is not satisfied in that period.

5.2. Cohesion

Fig.15 shows the influence of formation cohesion on sanding.For this,various values of cohesion for the reservoir are considered.It can be concluded that with lower cohesion, rock easily yields;consequently, more sand is produced. On the contrary, if the cohesion of the reservoir increases (e.g. by chemical stabilization),the sanding amount decreases. In this case, with increasing cohesion from 5 MPa to 8 MPa,no sand is produced.

5.3. Perforation orientation

Another parameter affecting sand production is perforation orientation. To investigate the effect of this parameter on sanding,the minimum horizontal stress and maximum horizontal stress are considered to be 41 MPa and 45 MPa,respectively.Fig.16 compares the impact of perforation orientation on sand production. The results show that higher sand is produced from the perforation parallel toSh(minimum horizontal stress) compared to that of the perforation tunnel parallel toSH(maximum horizontal stress).

5.4. Drawdown

Higher drawdown can result in higher fluid velocity. Consequently, drag forces exerted on sand particles will increase, which results in more sand production from the reservoir. Fig.17 shows the cumulative sand production for various drawdowns in this model,in which the base drawdown is 2.5 MPa.From the figure,it shows that when the amount of drawdown increases, the severity of the sanding also increases. In order to control this, either sand control methods should be used, or the drawdown should be adjusted to optimum level to prevent excessive sanding and exerting damage to the well.

Fig.15. Sanding amount for different reservoir’s cohesion values.

Fig.16. Total sanding amount with different perforation orientations.

Fig.17. The quantities of sand production for different drawdowns.

6. Conclusions

In this context, a hydro-mechanical sand erosion mechanism using a finite element platform along with a sanding criterion has been examined for simulation of sand production in the wellbore with or without perforations. All the aforementioned numerical simulations and analyses were based on the continuum-based approach. The numerical results were validated against laboratory experiment and field data. Following conclusions can be drawn:

(1) Sanding criterion that incorporates“plastic shear strain”and“erosion rate” can predict the onset of sanding and sanding amount with an acceptable level of accuracy.

(2) Producing sand from well and perforations causes large deformation that may induce severe mesh distortion in the model.ALE and mesh mapping techniques are suitable tools to achieve a high-quality mesh throughout the simulation.

(3) The effect of in situ stresses on sanding is high. Generally,smaller in situ stresses cause lower sand production, and in situ stresses increase with increasing depth. Therefore, by studying the effect of this parameter on sanding,the critical depth can be obtained.

(4) Stabilizing formation with resin or cement is one of the sand control methods. For the effect of this technique on sand production, sanding amount was obtained under various cohesion values. The result reveals that increasing cohesion could decrease sand production effectively. For example, for formation with cohesion of 8 MPa, no sand is produced.

(5) Drilling the perforations along the minimum horizontal stress direction could increase sand production.

(6) Effect of drawdown on sanding is also considerable. For instance,when drawdown increases from 2.5 MPa to 5 MPa,sand production increases by 120%.

Declaration of Competing Interest

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

The research was funded by the Iran National Science Foundation (INSF) (Grant No.96001589).

List of notations