Computational study of reservoir sand production mechanisms

A

cross-sectional area to the flow

c

concentration of fluidised solids transported

$Ǥ¯$

flow potential

g

acceleration due to gravity

ɧ

invariant related to the third invariant of the deviatoric stress tensor

I₁

first invariant of the stress tensor

J₂

second invariant of the deviatoric stress tensor

J₃

third invariant of the deviatoric stress tensor

k

rock permeability

L

length between the regions of pressure drop

$m˙$

rate of solid mass eroded

P

pore fluid pressure

P_p

pore pressure

q_i

fluid flux

$qˆuc$

uniaxial compressive stress

$R=$

ratio of the yield stress in triaxial tension to the yield stress in triaxial compression

S₁

maximum principal stress

S₂

intermediate principal stress

S₃

minimum principal stress

$Sijd,Sjid$

components of the deviator stress tensor

$Sjie$

equivalent stress

S_H

maximum horizontal principal stress

S_h

minimum horizontal principal stress

S_m

mean stress

S_o

cohesion

S^vm

von Mises stress

S_z

total vertical stress

$sˆ$

degree of saturation

$Ɏf$

fluid specific weight

z

height being considered

z_p

height of the phreatic surface

z_s

position of the surface of the rock mass

$zs′$

height of the surface that distinguishes between the dry and partially saturated zones

γ^p

plastic shear strain

ΔP

drawdown (pressure drop between the reservoir pressure and the wellbore pressure)

δ_ij

Kronecker delta

ϵ

void ratio

$ϵ¯pl$

equivalent plastic strain

$ϵ¯˙pl$

equivalent plastic strain rate

ϑ

dilation angle as measured on the meridional (S _m t) plane

λ

sand production coefficient

μ

fluid viscosity

ρ_d

dry density of the rock

ρ_f

fluid density

ρ_s

density of solids

σ_hyd

hydrostatic stress

σ_z

effective vertical stress

Φ

porosity

Φ_i

initial porosity

φ

angle of internal friction

Sand production is the generation of solid particles as a result of the disintegration of reservoir rocks which are weakened in the course of extraction of hydrocarbons. When it is not circumvented, this phenomenon has the potential to cause substantial losses in the quantity and quality of hydrocarbons produced. It also increases the cost of the operation (Rahmati et al., 2012; Wang et al., 2016), attributed to erosion of pipelines and other surface facilities such as valves, plugging of separators and the production liner by sand deposits, additional and more expensive intervention measures and a decrease in the rate of production (Rahmati et al., 2012; Tronvoll and Fjaer, 1994; Wang et al., 2016). The deterioration of equipment and well integrity has cumulative and culminating adverse implications, including wellbore failure and increase in disposal costs and downtime (Eshiet, 2012; Penberthy and Shaughnessy, 1992). The environmental consequences of sand production are severe, even in the short term. For instance, wear and tear of downhole and surface equipment due to the highly erosive properties of sand-laden fluids can be severe enough to cause total failure of downhole and surface facilities, leading to serious environmental, health and safety issues (Penberthy and Shaughnessy, 1992). Sand production may also collapse the formation because of significant reductions in pore pressure and compaction of the reservoir rock; the combination of both processes is likely to cause surface subsidence (Penberthy and Shaughnessy, 1992).

Sanding is generally instigated as a posterior event to failure of reservoir rock, which happens in response to in situ stress regimes and induced changes in stress conditions when hydrocarbons are being extracted from wellbores. This phenomenon is described by Wu et al. (2006a) as comprising three major processes occurring consecutively, starting with rock material failure at the vicinity of the wellbore, when local stresses surpass the tensile, compressive or shear strength of the rock. This is followed by disaggregation or dislodgement of particles at the failure zone and finally the movement of loose particles as fluidised materials. The same phenomenon can be depicted as just two processes. Morita and Boyd (1991) describe these as, first, the inducement of stresses and weakening of the rock at the surroundings of the wellbore by virtue of drilling activities, depletion of reservoir pressure, drawdown and so on and, secondly, the erosion of disintegrated materials. The latter categorisation by Morita and Boyd (1991) supports the recognition of the two key mechanisms that drive the sanding process. The initial phase is the mechanical instability and failure of the rock surrounding the wellbore, and the second and final phase is the hydrodynamic erosion of the failed material around the well cavity (Eshiet and Sheng, 2013; Li et al., 2018; Papamichos et al., 2001; Wang et al., 2016), which is then transported towards the wellbore in response to fluid flow and flow-induced pressure gradients (Wang et al., 2016). The second phase is also associated with hydromechanical instability (Eshiet and Sheng, 2013; Rahmati et al., 2013; Wang et al., 2016), since the effects of fluid flow, drawdown and pore pressure depletion alter localised stress concentrations, and the removal of solid particles augments the rock porosity while facilitating a realignment of the remaining solids and interparticle forces. Both phases are, hence, interrelated. Hydrodynamic actions enhance rock damage, thereby increasing the propensity for mechanical failure.

Several mechanisms impact the sanding process. Some of these factors are discussed by, for example, Younessi et al. (2013), Papamichos et al. (2010), Chin and Ramos (2002) and Nouri et al. (2002). These include water sensitivity of rock strength, capillary cohesion in the failed zone, multiphase flow and pore pressure gradient (Papamichos et al., 2010); triaxial stress conditions (Younessi et al., 2013); seepage, depletion, erosion, water-cut and material weakening (Nouri et al., 2002; Wu et al., 2006a) and fluid viscosity and producing time (Chin and Ramos, 2002). The geometry of the failure region near the wellbore is determined by Younessi et al. (2013) to be a function of lateral stress anisotropy and its magnitude. The width of this failure region is inversely proportional to lateral stress anisotropy; after yield, a critical drawdown (CDD) pressure is then necessary to induce sanding. Water breakthrough reduces capillary cohesion and weakens the strength of water-sensitive sandstones, thereby causing earlier onset of sanding. However, in two-phase flow conditions, the sanding rate is reduced and there is a rise in ‘sand production excess stress’ (stress magnitudes above the critical value for onset of sand production) due to capillary cohesion enabled by irreducible water saturation (Papamichos et al., 2010). The critical bottom-hole pressure (CBHP) increases with cohesion and reservoir pressure, whereas CDD decreases with increase in cohesion (Nouri et al., 2002). On the other hand, CBHP decreases with friction angle and elastic modulus, while CDD increases with friction angle and elastic modulus (Nouri et al., 2002).

Although a few analytical models exist, such as that proposed by Gholami et al. (2016), most non-experimental studies of sand production are realised through numerical computations. Numerical simulation of the sanding process is more frequently implemented by applying the continuum approach. This method relies on the underlying principles of continuum mechanics (e.g., Wu, 2005) where materials are treated, on the macroscopic scale, as a single and consistent entity using a fixed combination of parameters. With this approach, discontinuities and deformations, if required, must be explicitly incorporated (e.g. Nouri et al., 2009). Where microscale phenomena are of interest, discontinuum approaches are more suitable because of their inherent ability to capture local occurrences. The discrete element method (DEM), as a notable example of a discontinuum approach, has been applied by Cui et al. (2016) and Climent et al. (2013, 2014a, 2014b) to examine the impact of fluid flow and in situ stresses on the sanding process. A prime condition for a successful implementation of DEM is accurate couplings with computational fluid dynamics techniques.

Early studies on sand production were limited only to the mechanical failure of reservoir rocks (Rahmati et al., 2013; Vardoulakis et al., 1996), where the onset of sanding was attributed to mechanical instabilities of rock and localised failure. Examples of such works are presented by Morita et al. (1989), Vaziri (1988) and Perkins and Weingarten (1988). However, the concept of hydromechanical instabilities governed by both surface and internal erosion was later introduced by Vardoulakis et al. (1996). The contributions of seepage forces and other fluid-flow- and pressure-gradient-related parameters are now commonly integrated in sand production models. Currently, the mechanical behaviour of the rock material is generally extended to include its performance in the post-yield hardening/softening stage. This recognises the unique nature of the rock, particularly in the plastic range and how stress–strain relationships at this region contribute to the yield surface and failure envelope. The main rock failure models available include the Mohr–Coulomb, Drucker–Prager, Hoek–Brown, Mogi, Lade and Weibols–Cook models. A comparison of their performance is reported by Mehranpour and Kulatilake (2016). The Hoek–Brown and Mohr–Coulomb failure models are more popularly applied mostly because of their simplicity. This is evident in sand production studies, where various versions of Mohr–Coulomb failure criteria are preferred (e.g. Gholami et al., 2016; Hayavi and Abdideh, 2017a; Li et al., 2018; Nouri et al., 2006, 2009; Wang and Sharma, 2016; Wang et al., 2016). Interestingly, in the paper of Gholami et al. (2016), Mohr–Coulomb, Hoek–Brown and Mogi–Coulomb material failure criteria are compared when applied in a developed analytical elliptical model that predicts the volume of sanding based on the evolving shape of a borehole in the course of drilling and production. Both Mogi–Coulomb and Hoek–Brown failure criteria produce results that were somewhat closer to reality. However, the outcome using the Mohr–Coulomb criterion was unrealistic.

This study explores the mechanisms of key phenomena and their contributions to the reservoir erosion process by implementing the combination of a material failure model and an erosion model for solid mass generation. Factors such as pressure gradient, drawdown, pore fluid velocity, material strain and mud weight are examined. The investigation forms part of a research project on fracturing and reservoir erosion presented by Eshiet (2012).

A crucial aspect of this work is the departure from the more popularly adopted material failure model – that is, Mohr–Coulomb – and the Hoek–Brown model, in favour of the Drucker–Prager model. The Mohr–Coulomb criterion has a wider range of application in geo-engineering problems and is relatively simpler; however, it has two major shortcomings. First, it neglects the influence of the intermediate principal stress, and, with respect to this, the major principal stress is assumed to act independently. This results in an underestimation of the material yield strength (Jiang and Xie, 2011). Secondly, the shape of the yield surface on the deviatoric plane is an uneven hexagon with six acute angles that restrict the convergence inflow theory (Jiang and Xie, 2011). The yield function as defined by the Drucker–Prager model smoothens the surface defined by the Mohr–Coulomb criterion as well as modifies the von Mises criterion (Öztekin et al., 2016). Moreover, the Drucker–Prager model accounts for contributions from the intermediate principal stress, which, unlike counterpart Mohr–Coulomb models, underscores the fact that, for geomaterials, the uniaxial compressive strength is lower than the biaxial compressive strength.

Another distinctive component of this work is the post-yield behaviour of the rock. The assumption of an elastic–perfectly plastic behaviour for rocks is now regarded as misrepresentative, since they enter a new set of regimes once the yield point is reached. Hardening/softening models are therefore requisite representations of the post-yield plastic behaviour of rocks. The mobilisation of certain properties serves as a determinant and indicator of strain hardening and softening. During strain hardening, the friction angle is mobilised as it increases due to the accumulation of plastic strain and becomes fully mobilised and thereafter remains constant only after the peak strength value is attained (Nouri et al., 2007). Strain softening is the succeeding regime. At this the stage (when the maximum strength is reached), the rock cohesion, which is, hitherto, constant, is mobilised by softening and the tensile strength decreases, instigating the formation of microcracks (Nouri et al., 2007; Papanastasiou and Vardoulakis, 1992; Vaziri et al., 2008). This allows for a formulation denoting an expansion and a contraction of the yield surface during strain hardening and softening, respectively (Vermeer and de Borst, 1984). Post-yield behaviour is a function of the triaxial compressive state of the rock: the strain hardening regime is predominant at high effective confining stress conditions, while a more extended strain softening regime is more likely in low effective confining stress conditions (Nouri et al., 2007). In this study, the strain softening regime is modelled as being constricted at very high confining stress states; the tension cut-off is a function of the plastic shear strain, reducing rapidly down to zero at the softening phase. Post-yield isotropic hardening is predetermined and explicitly prescribed according to Eshiet and Sheng (2013).

The other characteristic facet is the use of the equivalent plastic strain (⁠$ϵ¯pl$⁠) instead of the plastic shear strain. The equivalent plastic strain allows an arbitrary three-dimensional (3D) strain state to be represented by a sole parameter.

Static stress/displacement analyses were conducted using finite-element method (FEM) techniques. These were implemented through Abaqus, a software suite for finite-element analysis (Dassault Systèmes, 2018). Finite-element formulations are conventionally anchored on concepts of continuum mechanics (solid and fluid mechanics), where stress and displacement fields are considered to be continuous within the element of interest (Pian and Tong, 1972). An implicit integration scheme (a characteristic feature of ‘Abaqus/Standard’) is adopted, where the relevant sets of equations are solved involving both the current and later states of the model system, at every solution increment. This numerical scheme is ideal for static and slow dynamic activities, which typically have high stiffness.

Well completion is the process whereby a newly drilled well is prepared for production. This includes a sequence of activities such as the installation of the various types of casing, cementing, perforating and gravel packing and the installation of completion components – for example, the Christmas tree, tubing hanger and downhole gauges and valves. A well may be completed either by allowing it to remain open (top sets or barefoot completion) or by casing it (case completion). It is imperative for casing completions to be accompanied with perforations to permit inflow of fluids. Thus, cased wells consist of two key components: the main wellbore and a series of perforation tunnels perpendicular to the axial axis of the wellbore, and equidistanced vertically/horizontally and azimuthally from each other. A 3D representation of a cased well is depicted in this study (Figure 1). One direct derivative of the principle of symmetry is the possibility of using only a quarter section of the model instead of constructing the entire domain. The perforator phasing is set at 90°; therefore, this design aligns four perforation tunnels spaced at right angles apart in each longitudinal well segment. This implies that only one perforation tunnel is included in a quarter of the segment.

The geometric dimensions of the reservoir domain, wellbore and perforation tunnel are given in Table 1. These dimensions, although arbitrary, were chosen in accordance with typical well geometries in oil/gas fields. Typically, the length of perforation tunnels may exceed 2 m, and its actual depth is influenced by several factors, including the following (Seibi et al., 2008): initial speed of the discharged explosives, effective surface area, casing strength, cement strength and type of formation. The penetration depth is inversely proportional to the effective surface area and directly proportional to the initial speed (Seibi et al., 2008). Two kinds of elements were used to build the finite-element model (Eshiet, 2012): the rock (main domain) was created with continuum, 3D, eight-node, trilinear displacement and pore pressure elements (C3D8P), while the cement sheath was built with 3D, four-node, full integration membrane elements (M3D4). These solid and membrane elements are defined by first-order (linear) interpolation. First-order elements were preferred to high-order elements because of their ability to generate accurate results without significant oscillations at a lower computational cost. In contrast to higher-order elements, they are fundamentally constant-strain elements; nonetheless, this limitation is not appreciable in this case because of the relatively high rate of displacement involved in the failure and erosion process. Non-constant and linear strain responses provided by individual higher-order elements are essential, for instance, in problems relating to bending and compression. Generally, in this study, the piecewise linear solutions provided by first-order elements yield similar resolutions in comparison to second-order elements.

A schematic diagram of the formation domain is sketched as in Figure 1(c). This is the segment of the rock used for the analysis. There are seven outer boundaries, labelled as ‘Face-1’, ‘Face-2’, ‘Outer radial’, ‘Wellbore face’, ‘Face-top’, ‘Face-bottom’ and ‘Perforation face’ (Figure 1(c)). Initial and boundary conditions are defined in terms of stress, displacement and pore pressure. Initial in situ stress conditions are specified in Table 5, comprising horizontal and overburden effective stresses given as 34·47 and 51·71 MPa, respectively.

Face-1 and Face-2 are boundary surfaces created as a result of the separation of the quarter section from the remainder of the full model (Figure 1(c)). Displacement at the boundaries is controlled by constraining these surfaces in the second degree of freedom – that is, in the tangential direction (θ = 0). Also, the outer radial and wellbore surfaces, ‘Outer radial’ and ‘Wellbore face’, are constrained in the first degree of freedom (i.e. radial direction, r = 0), while the bottom surface, ‘Face-bottom’, is constrained vertically. The entire rock domain, as well as the outer radial surface and the perforation tunnel surface (Perforation face), is assigned an initial pore pressure of 37·92 MPa.

The horizontal and vertical stresses presented in Table 5 are in situ and reflect the initial effective stress regime at the reservoir. The initial reservoir depth considered is approximately 3450 m (i.e. 3448 m), corresponding to the total vertical stress, S_z (i.e. $σz+Pp$⁠), of 89·63 MPa. A saturated bulk density of 2650 kg/m³ is assumed for sandstone (Manger, 1963; Papamichos and Stavropoulou, 1998) and an initial void ratio (ϵ) of 0·351 is associated with a porosity (Φ) of 0·26 (Equation 11).

Brittle rocks generally exhibit high friction and dilation angles (e.g. Jafarpour et al., 2012; Zhang et al., 2016). This is more so in rock masses due to contributions from discontinuities. At the subsurface, factors such as confining stress affect rock dilation. Dilation increases with decreasing confining stress; however, irrespective of the level of confining stresses, the dilation angle and friction angle are sometimes mobilised to very high peak values during hardening. Dilation and friction angles of 40 and 45°, respectively, were adopted here as given, for instance, by Jafarpour et al. (2012) for sandstones. Friction angles of brittle rocks can be mobilised up to approximately 80° (e.g. Renani and Martin, 2018).

The failure behaviour of the rock (sandstone) material is described using a linear Drucker–Prager model with hardening, and the cement sheath considered to be linearly elastic. Material properties used for the model include parameters for the Drucker–Prager model, permeability, void ratio, specific weight, elastic modulus and Poisson’s ratio. These and others are given in Tables 2 and 3.

Models describing material behaviour are generally selected or formulated based on factors including the kind of material, the analysis to be conducted, the availability of experimental data, loading conditions and the magnitude of stresses expected to be encountered (Eshiet, 2012). The Drucker–Prager model was chosen for reasons that take into account its ability to describe the following: the behaviour of frictional and brittle materials including pressure-dependent yield; post-yield isotropic hardening and softening; the behaviour of materials with higher compressive strength in comparison with tensile strength; and the volume change associated with material behaviour in the inelastic (plastic) range. In addition, the preference for this failure model was encouraged by its ability to incorporate elastic material models, particularly where small strains are anticipated.

The yield criterion in Drucker–Prager models is dependent on the shape of the yield surface, categorised into three forms: linear, hyperbolic and exponential. On the meridional plane, the differences between these forms are mostly pronounced at low confining stresses. The linear form was deemed adequate since the modelling conditions entailed perpetually high confining stresses. It assumes that the deviatoric stress and shear stress are linearly dependent on the equivalent pressure stress (representing the mean stress), which is similar to the relationship described in the hyperbolic and exponential forms at high confining stresses. Since the material was expected to show a post-yield hardening behaviour, an extension to the Drucker–Prager model was made to capture the increase in yield strength during plastic deformation and increase in plastic strain or strain rate if rate dependency is considered. Being a quasi-static event, rate dependency was precluded and the hardening process was made a function of the extent of plastic straining. The hardening behaviour was monitored with respect to changes in the equivalent stress, defined in this case as the uniaxial compressive yield stress, expressed as

where $qˆuc$⁠, $ϵ¯pl$ and $ϵ¯˙pl$ are the uniaxial compressive stress, the equivalent plastic strain and the equivalent plastic strain rate, respectively. The yield surface can be described using invariants.

The hydrostatic stress, σ _hyd, is the mean of the normal stress components of the Cauchy stress tensor. The hydrostatic stress is negative when expressed in terms of the equivalent pressure, which is compressive. The equivalent pressure stress tensor, S _m, represents the mean stress given in terms of the first invariant of the stress tensor.

where I ₁ is the first invariant of the stress tensor written in terms of the principal stresses as

The deviatoric stress is represented using the equivalent stress or von Mises stress, related to the second invariant of the deviatoric stress tensor as

where J ₂ is the second invariant of the deviatoric stress tensor, expressed as (e.g. Malvern, 1969)

where $Sijd$ and $Sjid$ are components of the deviator stress tensor; $S1d$⁠, $S2d$ and $S3d$ are the maximum, intermediate and minimum principal deviatoric stresses, respectively; and S ₁, S ₂ and S ₃ are the maximum, intermediate and minimum principal stresses respectively. J ₂ can be defined as a function of either the deviatoric stress tensor or the associated principal stress values. The deviatoric stress is calculated through Equations 4a and 4b in terms of the equivalent pressure and the hydrostatic stress, respectively

where S _m is the mean stress as earlier defined. Another invariant used in the linear model is related to the third invariant of the deviatoric stress tensor J ₃, used in the form

where

The yield criterion for the linear Drucker–Prager model is

where φ and S _o are the angle of internal friction and the cohesion, respectively, and t is determined with respect to the invariants of the deviatoric stress tensor. t is a measure of the deviatoric stress. It can be plotted as a substitute for the deviatoric stress in the deviatoric–mean stress (S _m t) plane to define the failure surface and is given as

t is a parameter representative of the deviatoric component of the stress tensor (consisting of shear stresses) on the S _m–t plane. It is defined in Equation 6b with respect to the second and third invariants of the deviatoric stress tensor (S ^vm, ɧ). $R=$ is defined as the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression. The plastic flow of the material is modelled using the flow potential (⁠$Ǥ¯$⁠)

where ϑ denotes the dilation angle as measured on the meridional (S _m t) plane. During hardening, flow occurs at an angle ϑ to the t axis.

Parameter values for the Drucker–Prager model and material properties (Tables 2 and 3) were selected to lie within the range of values typically representing sandstones (Eshiet and Sheng, 2013). Material properties for the cement sheath and the domain, the well dimensions and the initial conditions are given in Tables 1, 4 and 5. Also, the pore pressure gradient (Figure 2) is positive when referenced from the wellbore and the highest values of pore pressure occur at the extreme boundaries. This is necessary to instigate flow in the direction of the well as would normally occur during petroleum production.

Contrary to models adopted by some others, the onset of sanding is hereby defined by a distinct criterion in variance with the more general failure criterion that describes rock behaviour. Some erosion prediction models attempt to synchronise the initial shear or compressive failure of the rock material with initiation of sanding; however, this assumption has been established as overly conservative. Two criteria are, hence, depicted in this model: a material failure criterion effectively described by the Drucker–Prager model and a sanding criterion given in an eroded solid mass generation equation formulated by Papamichos and Stavropoulou (1998) 

where $m˙$ is the rate of solid mass eroded, ρ _s is the density of solids, λ denotes the sand production coefficient, Φ is the porosity, c is the concentration of fluidised solids transported and q _i is the fluid flux. The left-hand term of the equation, $m˙/ρs$⁠, is denoted as the erosion velocity, v _e, and the term (q _i q _i)^1/2 is given as the pore fluid velocity, v _fp.

As shown by experiments, initiation of sanding occurs when a critical external stress value is exceeded, which is incorporated in Equation 9 and makes the sand production coefficient (λ) dependent on the plastic shear strain (γ ^p) (Papamichos and Stavropoulou, 1998; Papamichos et al., 2001). This implies that erosion can take place in the rock material only when its maximum strength is surpassed and the failure regime is in the plastic softening stage. The erosion capacity increases with plastic strain, since λ depends on the plastic strain. This occurs at a rate λ ₁, increasing until a maximum value λ ₂ is attained. Papamichos and Stavropoulou (1998) and Papamichos et al. (2001) assume the following function for λ = λ(γ ^p)

In this analysis, the plastic shear strain is represented by the equivalent plastic strain (⁠$ϵ¯pl$⁠). The equivalent plastic strain is the integral of the equivalent plastic strain rate, given as

whereas the equivalent plastic strain rate is given in terms of the plastic strain rate, $γ˙p$

For known values of void ratio (ϵ), the porosity, Φ, is determined using the relation between the two parameters, which is expressed as (Knappett and Craig, 2012: p. 23)

rewritten as

Although the formulae presented here are appropriate for volumetric erosion, the premise on which the equations are derived has been carefully adapted by applying the surface velocity in order to fit the processes encountered in surface erosion (Eshiet, 2012). Surface erosion is caused by the hydrodynamic action of fluid when it interacts with the rock surface; this excludes the disintegration and flow of loose materials which can take place inside the material (Eshiet, 2012). Emphasis is given to the surface erosion process, since there is a high tendency for the occurrence of larger stresses at the surface of the well, which causes a predominance of surface erosion at the vicinity of wellbores.

To capture the erosion process, a mesh adaptivity technique referred to as arbitrary Lagrangian–Eulerian (ALE) adaptive meshing was employed that entails periodic mesh smoothing using a combination of Lagrangian and Eulerian analyses. This method is particularly essential in cases where there are excessive deformations, making it necessary to check the distortion of elements. To apply the ALE adaptive meshing, an adaptive mesh domain was defined after the drilling operation. This domain consists of regions surrounding the perforation tunnel (Figure 1(b)). The mesh density of this region was made finer than the outer domain due to its sensitivity. The ALE mesh domain consists of interior and boundary regions. The boundary regions comprise nodes subdivided into free surface nodes, corner nodes, edge nodes and constrained nodes, which are all subject to the relevant mesh constraints applied (Table 6).

Interior nodes refer to nodes completely confined by elements within the adaptive mesh domain and have an unconstrained direction of motion. Surface nodes of boundary regions are constrained from movement in the normal direction, but are free to slide tangentially. Edge nodes of boundary regions are only able to slide on the edge with their position determined by the locations of adjacent edge nodes. Corner nodes are positioned at vertices, and their movement due to mesh-smoothing operations is constrained. Nodal displacement is controlled by applying adaptive mesh constraints in preferred directions.

Nodal motion is either defined explicitly or controlled by a mesh-smoothing algorithm with constraints determined by the boundary. To control node motion explicitly, it is necessary to apply spatial mesh constraints, which may be applied on any node apart from those in meshes with Lagrangian constraints. In this work, both Lagrangian and spatial type constraints were used. A spatial mesh constraint was applied using a subroutine, where the constraints rely on given information from nodes or material points. There are mainly two types of spatial mesh constraints: displacement/rotation and velocity/angular velocity spatial constraints. Velocity spatial constraints were applied because of the need to control the rate of movement of nodes at certain regions (perforation tunnel surface). This was accomplished by repeated adjustment of adaptive mesh velocities using calculated values of erosion velocity, carried out during mesh-smoothing operations. Magnitudes of the adaptive mesh velocity were determined by a mesh-smoothing algorithm.

The erosion equation used to compute the erosion velocity is solution dependent and therefore frequently requires updated values of output variables including plastic strain, fluid velocity and porosity. Porosity was calculated using values of the changing void ratio (Equations 11a and 11b). Values of these output variables were obtained from material integration points of elements and the actual nodal values used inferred from them. Lagrangian type constraints were applied at the interface between the adaptive domain and the rest of the rock mass (Figure 1) to prevent the occurrence of mesh smoothing and to indicate that the material within the area should follow the underlying material in line with nearby elements of non-adaptive regions.

Thus, in line with the research objective, the ALE adaptive meshing technique was adopted to account for the wearing and loss of rock materials due to the erosion process. A flow diagram illustrating the sequence of procedures for controlling the behaviour of the adaptive mesh domain and computing rock ablation is demonstrated in Figure 3.

For nodes at the outer sides of the ALE adaptive mesh domain (top, bottom, left and right surfaces), adaptive meshing is not carried out in the direction in which boundary conditions are applied but in other directions (i.e. the radial and vertical directions). Thus, the outer boundary regions of the ALE adaptive mesh domain consist of constrained nodes, contrary to nodes (interior nodes) inside the interior region (Figure 4(a)).

The motion of the nodes in the mesh of the ALE adaptive domain is controlled by a mesh-smoothing algorithm with constraints enforced by the boundary conditions (at the boundary regions). In other parts of the ALE domain, the motion of the nodes is explicitly defined. In this regard, two types of ALE adaptive mesh constraints are imposed: Lagrangian and spatial mesh constraints. The former is applied at the inner ALE adaptive domain boundary (Figure 4(b)) and the latter at the surface of the perforation tunnel (Figure 4(b)). The motion of nodes subjected to spatial mesh constraints are defined by periodically updated mesh velocities; these are, in essence, velocity spatial constraints.

The model analysis was completed in five steps, described as follows.

Within this step, the initial geostatic stress field was defined and equilibrium was established in order to represent an undisturbed rock material under steady-state equilibrium and subjected to geostatic loading. This was effectuated after prescribing initial values of void ratio, pore pressure, effective stresses and distributed loading serving as the overburden rock layer. The loads and initial stresses equilibrated, thereby preventing deformations. During the geostatic procedure, atmospheric pressure was ignored, and, since there was no initial flow, the pore fluid was assumed to be in hydrostatic equilibrium. Pore fluid pressure changes linearly with depth. It is zero at the water table/phreatic surface, increases linearly with depth below this level and becomes negative above the water table/phreatic surface. For an incompressible fluid with specific weight independent of depth, the pore fluid pressure is defined thus

where $Ɏf$⁠, ρ _f and g are the fluid specific weight, fluid density and acceleration due to gravity, respectively. Also, z _p denotes the height of the phreatic surface, while z is the height in consideration. Equation 12 ensures that the pore fluid pressure is negative above the phreatic surface. Shear components of stresses are considered negligible, which allows the assumption of horizontal stresses that vary with depth, but remain constant horizontally. The term ‘horizontal stresses’ is used as a collective word for both radial and circumferential stresses. If shear stresses are ignored, vertical equilibrium is achieved under the following condition (Dassault Systèmes, 2018)

Equation 13 is rewritten in terms of effective stresses as follows

where S _z and σ _z are the total and effective vertical stresses, respectively; ρ _d is the dry density of the rock; $sˆ$ is the degree of saturation; and Φ_i is the initial porosity. ρ _d and Φ_i are taken to be constant. Vertical equilibrium is achieved when conditions for Equations 13 and 14 are met. Equations 13 and 14 denote the differences in the total and effective stresses, respectively, with respect to the height of the point (z) being considered. These equations are expressed as functions of the degree of saturation and porosity. The first and second terms at the right-hand side of Equation 13 are the effective stress and pore pressure components, respectively. If z is designated as the height in consideration and $zs′$ is the height of the surface that distinguishes between the dry and partially saturated zones, Equation 14 is solved to become

where z _s is the position of the surface of the rock mass. The horizontal stresses (S _H, S _h) are assumed to remain constant horizontally, but they are dependent on vertical variations in vertical stresses.

In this step, the wellbore and perforation tunnel was removed by a contact deactivation procedure. In order to maintain the stability of the new configuration, pressure was applied on the perforation face.

New boundary conditions were set in this step to apply pore pressure on the perforation tunnel face. At this stage, the pore pressure at the perforation face is the same as that in the far-reach region (outer radius).

Drawdown was instigated by reducing both the pore pressure at the perforation face (changing the boundary condition) and the applied pressure (mud pressure) on the same face.

Erosion was simulated at a constant drawdown pressure using ‘adaptive meshing’. To enable this process, an adaptive mesh domain was defined at the radial region close to the wellbore and perforation tunnel. Erosion was described through a spatial adaptive mesh constraint applied to all the nodes located at the surface of the perforation tunnel. A velocity spatial adaptive mesh constraint was specified at the perforation surface.

FEM analyses were used to study the mechanisms of subsurface erosion due to fluid flow, particularly at the macroscopic scale. This section comprises two parts: sensitivity analysis and validation. Within the first part, sensitivity analyses are carried out to determine the influence of mesh density (element size) on results obtained from the FEM simulations. Validation of numerical results by comparisons with other results forms the basis of the second part.

The mesh dimensions were mainly adjusted in the circumferential direction. The pattern of meshing was skewed so as to become finer towards the well and perforation regions. The finest meshes were placed at this region due to the sensitivity and importance of events occurring there. The region was also delineated and defined as an adaptive domain (described in the section headed ‘Modelling the erosion process’). The mesh density is defined in terms of the number of elements generated, and, for this analysis, it is equivalent to the number of elements within a specified region. Results of various models built using meshes of different densities are presented in terms of three parameters: the von Mises stress, the deformation and the erosion rate. The first two parameters are presented as contour plots (Figures 5–7) so as to illustrate the spatial distribution of the parameter quantity as well as to emphasise the differences in mesh density. The third parameter (erosion rate) is presented in a convergence plot (Figures 7 and 8).

The initial, boundary and operating conditions are similar to those described by Eshiet and Sheng (2013). A constant drawdown of 3·72 MPa was maintained throughout the simulation period. For all cases, the edge seeds at the outer region (away from the adaptive zone) were set to have a bias ratio of 5 to enable consistency in the element aspect ratio at this region. The outer region in this context refers to the far-reach zone beyond the adaptive mesh domain (Figure 4). Edge seeds are markers placed along the edge of an area of the model domain to prescribe its mesh density. The bias ratio is the distance between nodes of elements (i.e. element length) at the extreme end of the coarse mesh section divided by the distance between nodes of elements at the extreme end of the dense mesh section of the edge.

The spatial distribution of both von Mises stress and deformation is consistent for the different mesh densities, except for the higher values of deformation noticed when the adaptive region is discretised into 900 elements. Likewise, the rates of erosion (sand production) match, as depicted in Figure 8.

Additional tests were conducted whereby the changes in mesh density were restricted only to the adaptive zone, keeping the mesh density at the outer region constant. In all cases, the maximum element aspect ratio at the ALE adaptive mesh domain was less than 5. The results show variances in the sand production rate, with higher values recorded as the mesh density of the adaptive zone is reduced (Figure 9). Convergence was observed at a mesh density of 12 480. Below this number of elements, there was no apparent change in the rate of sand production. The divergence in the results when higher mesh densities (>12 480) were used at the adaptive zone is attributed to the widening disparity in element sizes, particularly at the boundary between elements in the adaptive zone and the other domain. The disparity increases with higher mesh densities. The fixed mesh pattern and density at the outer region prevent smooth transitions between boundary elements, with the transfer of information from the nodal and integrated points becoming less accurate with decreasing element sizes at the adaptive zone. The selection of appropriate mesh densities was therefore based on the initial convergence test portrayed in Figure 8. The primary region where the erosion process occurs is the area at and around the well face and perforation tunnel, denoted as the adaptive zone. More attention was therefore given to this area due to its importance. Based on the results of the mesh sensitivity analysis, an adaptive mesh density of 9507 was deemed to be appropriate as it falls within the limits of convergence. Also, a mesh density of 21 253 for the whole model was used.

Results from the FEM (erosion) model show trends comparable to those reported in previous research (Papamichos and Stavropoulou, 1998; Papamichos and Vardoulakis, 2005; Papamichos et al., 2001). The parameters examined include sand production, pore pressure, pore fluid velocity and plastic strain. These were referenced from the wellbore region. Patterns of cumulative sand production as functions of time were compared with results presented by Papamichos et al. (2001) and Papamichos and Stavropoulou (1998) (Figures 10 and 11, respectively). Although the quantity of eroded material increases with time, the rate of erosion is transient and may decrease following phenomena that could lower the susceptibility of the formation material to wellbore instability. Wellbore instability is an antecedent of wellbore collapse, fracturing, wellbore closure or wellbore enlargement/washout. This is noticed in Figure 11, where the cumulative amount of sand produced approaches a constant value at later stages. The two patterns (Figures 10 and 11) are typical of the transient nature of the erosion process. The effect of external stress conditions on the magnitude and intensity of sand production is highlighted in Figure 12. The extent and intensity of sand production increases with external stresses, and the similarity of this trend is depicted by the developed FEM model as well as the model by Papamichos and Stavropoulou (1998). This is a qualitative comparison, demonstrating a similarity between patterns due to increasing external stresses/vertical stresses. Likewise, the magnitude and intensity of sand production increases with the flow rate. Note that vertical pressures, as denoted by the FEM model, are equivalent to vertical stresses, which are forms of external stresses in this context.

The similarity in trend between the two models is illustrated in Figure 13. Drawdown is directly associated with flow rate because it controls the pore pressure gradient, which is, in turn, related to the fluid flow rate. In addition to the trend, the magnitudes of sand produced, as shown in Figures 11 and 13, are of the same order. The quantity of sand produced was deliberately converted from cubic metres to grams to match units adopted in the relevant literature. A sandstone density of 2650 kg/m³ (Papamichos and Stavropoulou, 1998) was used.

In Figure 13(a), drawdown is directly linked to flow rate since it governs the pore pressure differential between the external and internal boundaries. Changes in flow rate occur according to prescribed drawdown conditions, and they are directly proportional to each other. Thus, flow rate increases with drawdown and vice versa.

The strong link between these two parameters is the premise for the comparison between Figures 10(a) and 10(b). The difference between the reservoir pressure and the bottom-hole (wellbore) pressure is the reservoir drawdown. Drawdown is the driving force that propels and steers flow into the wellbore. The relationship between drawdown and fluid flow rate is best represented by Darcy’s law, which is commonly used to describe hydrocarbon reservoir fluid flow. In its simplest form – applied in single-phase oil flow – the flow rate in a linear flow system (q) is given by

k is the rock permeability, μ is the fluid viscosity, A is the cross-sectional area to the flow, L is the length between the regions of pressure drop and ΔP is the drawdown (i.e. pressure drop between the reservoir pressure and the wellbore pressure). The proportionality constant is defined by the combination of two independent properties, k and μ. Ideally, this constant linearly relates reservoir fluid flow rate to drawdown. The relationship is also positive. Several other experimental observations and more sophisticated expressions exist which support this relationship – for instance, in the study by Song et al. (2015), where a positive and almost linear correlation between drawdown and flow rate is illustrated.

The resemblance in pore pressure distribution between the model by Papamichos and Stavropoulou (1998), the model by Papamichos and Vardoulakis (2005) and the developed FEM model is portrayed in Figure 14. The figure shows a pore pressure gradient that becomes steeper at the wellbore region. The values of the pore pressure in Figure 14 are functions of the applied external pore pressure and wellbore pressure in each case. The applied external pore pressures in Figures 14(a)–14(c) are 1·5, 0·4 and 37·92 MPa, respectively. Figure 15 depicts the spread of pore fluid velocity, which generally follows a power-law distribution, decreasing away from the wellbore. The plots (Figure 15) for both Papamichos and Vardoulakis (2005) and FEM models have similar features.

The exact magnitudes of the vertical pressures for Figures 14(a) and 14(b) and 15(a) are not known, but the trends are analogous to the FEM model results. Pore pressure and pore fluid velocity distributions are not affected by changes in vertical pressures (Figures 18 and 23).

Figures 14(a)–14(c) are qualitative comparisons of the trend in pore pressure distributions in the radial direction. They juxtapose results derived from this study with those by Papamichos and Vardoulakis (2005) and Papamichos and Stavropoulou (1998). The unit for pore pressure is megapascal, and, for radial distance, metre or centimetre. The order of scale for pore pressure is approximately similar for the three illustrations, bearing in mind that, unlike others, the bottom-hole (wellbore) pressure for the model used in this study is 34 MPa (Figure 2). The actual range of the active region of the pore pressure scale in Figure 14(c) is 4·0 MPa. The scale for the radial distance presented in Figure 14(c) is reflective only of the size of the reservoir domain considered. A reduction in this scale will not impinge on the pattern of pore pressure distribution.

The influence of certain key parameters was examined to ascertain their effect on the failure and subsequent erosion of the rock material. These include the following: drawdown, depth of wellbore (perforation depth) and erosion criterion. Also, the optimal mud pressure for given operational and reservoir conditions was determined.

The lithostatic or overburden pressure (vertical stress) at any point within the subsurface is representative of the depth of that location referenced from the ground or sea surface. A direct positive linear correlation can, hence, be drawn between the depth of a level at the subsurface and the vertical stress at the same position. Thus, the depth of the bottom-hole section of the wellbore is regarded, in this instance, to be the location of the perforation tunnel and is given by the magnitude of the vertically downward pressure at that level. As depicted in Figure 16, sand production increases with depth, with the highest amount occurring at a vertical pressure of 103 MPa. At and below a vertical pressure of 90 MPa, the quantity of sand produced is relatively very low and declines to a much reduced rate earlier (shortly after 1 d). A further increase in vertical pressure to 97 MPa and beyond results in large increments in both the quantity and rate of sand production (Figure 16). These boosts in the sanding phenomenon occur in tandem with vertical pressure, with the highest values ensuing, in this case, at a vertical pressure of 103 MPa (Figure 16). These observations are analogous to test results from other studies (e.g. Papamichos and Stavropoulou, 1998; Papamichos et al., 2001), where an increase in external radial stress (Papamichos and Stavropoulou, 1998) and confining stress (Papamichos et al., 2001) is shown to cause a corresponding rise in sand production rate. This parallel analogy is supported by the strong linkages between in situ stresses. The relationship between the vertical stress and horizontal stresses (i.e. circumferential and radial stresses) is complex but can be simplified for the homogeneous subsurface – for example, using proportionality terms involving Poisson’s ratio or friction angle (Kaiser et al., 2016; Mayne and Kulhawy, 1982). Generally, changes in vertical stress induce likewise effects on horizontal stresses at the same location.

Several mechanisms are responsible for augmentations in sand production with depth and confining stresses. The first mechanism is the higher plastic strains which are likely to occur at deeper levels as in situ and induced stresses surrounding the wellbore become greater. The sanding phenomenon is largely dependent on a range of factors. One of the crucial contributing components is the plastic shear strain; any increment in this parameter increases the propensity for the onset of sanding as well as the severity of the erosion process. For this study, the rock plastic strain is directly linked with the sand production coefficient (λ) (Equations 9a–9c), which is an experimentally derived parameter. In this regard, a pair of sand production coefficient parameters is required, λ₁ and λ₂, with the potential for erosion growing with plastic strain until the peak value of λ(λ₂) is attained. Due to elevated activities such as high fluid flow rates at the vicinity of the perforation tunnel and wellbore, rock deformation is most likely to take place at this region (Figure 17), escalating the chances of the material drifting into the plastic regime.

The second mechanism is guided by the yield and post-yield paths and the mode and rate at which material failure occurs. The post-yield behaviour adopted here is defined by strain hardening denoted by an increase in the rock yield strength in the course of plastic deformation or plastic straining or increases in strain rate, where rate dependency is taken into account. At deeper subsurface levels, the wellbore is more prone to instabilities and there is a higher tendency for rock material failure as a result of the interplay between in situ and induced stresses, pore pressure and fluid flow. The occurrence of shear, compressive or tensile failures depends on the kind of prevailing forces near the wellbore and any mitigating measures that are applied. The course of rock material failure is a crucial stage that serves as a precursor for subsequent processes of reservoir sand erosion.

The third mechanism is influenced by rock compaction and/or consolidation, which is more likely to happen at greater depths as a direct consequence of overburden and lateral confining pressures. Compaction increases with depth (Bjorlykke, 1978, 2015), which encourages a similar pattern of reductions in both porosity and void ratio. A denser rock material implies that more material would be available to be eroded.

The fourth mechanism is instigated by the drawdown which not only is prevalent at greater depths, but also depends on the upstream pressure condition of the reservoir. Drawdown is the differential pressure that causes fluid flow near the borehole. Fluid flow rate at this location increases with drawdown; the high pressure gradient renders it more vulnerable to erosion (Alshmakhy and Maini, 2012).

There is a threshold depth beyond which the sand production rate is rapidly accelerated (Figure 16). In this work, this exists between vertical pressures of 90 and 97 MPa. There is a considerable rise in the extent of erosion when the vertical pressure is increased from 90 to 97 MPa, which indicates the existence of a critical vertical pressure (and depth) that signifies the onset of the quickening of the sanding process. The point of a sudden rise in sand production rate is, herein, termed the ‘critical depth’. Determining the critical depth should be a crucial aspect of the design process of well operations. This would inform decisions concerning the ideal location of perforation tunnels and how this relates to other design parameters.

High velocities are noticed at the wellbore/perforation region which tend to decrease with increasing depth (Figure 18). The reason for the sudden increment in pore fluid (interstitial) velocities around the wellbore is attributed to, among other factors, the pore pressure distribution and changes in drawdown conditions, which will be discussed later in this section. In addition, Figures 19(a) and 19(b) show a clear increase in the surrounding pore fluid velocity when the vertical pressure (depth) is reduced. This is easily noticed at the tip of the perforation tunnel.

As earlier mentioned, drawdown is the differential pressure simply denoted as the difference between the reservoir pressure and the wellbore flow pressure. It is an important operational parameter because of its direct relationship with reservoir flow towards the wellbore. The wellbore skin is the immediate and surrounding area around the wellbore, which may be in different states of permeability. Ideally, this zone should be sufficiently permeable to enhance productivity. A positive skin condition suggests that rock material around the wellbore is somewhat impaired and well productivity is low. Contrarily, a negative or very low positive skin condition indicates that the rock material at the vicinity of the wellbore encourages high productivity (e.g. Jianchun et al., 2014). In cases of positive wellbore skin conditions, drawdown manipulation serves as an effective way of improving the productivity of the well and has been shown to expedite the recovery of depleted heavy oil reservoirs (Alshmakhy and Maini, 2012). As demonstrated, for instance, through Equation 16, reservoir pressure drawdown is linked with and has a direct impact on fluid flow rate. A high drawdown increases the pressure gradient near the wellbore, causing corresponding increments in the fluid flow rate in the same region. On the other hand, a lower drawdown reduces the pressure gradient and fluid flow rate. During well completions, one of the objectives is to increase inflow to the well, which can be achieved by instituting a sufficient pressure drawdown environment; however, increasing drawdown may aggravate sand production and wellbore collapse (Tronvoll and Sønstebø, 1997). Figure 20 shows the cumulative sand production with time at various drawdown conditions. Sand production increases with drawdown, and a drastic increase in eroded sand is noticed when a constant drawdown of 10·34 MPa is applied.

The link between drawdown and fluid flow is such that, for a given proportionality constant, an increase in drawdown accelerates fluid flow rate and vice versa. This applies to Darcy and non-Darcy flows. The proportionality constant primarily consists of both rock permeability and fluid viscosity. It is safe to assume that these properties are constant and uniform in homogeneous single-fluid-phase reservoirs. High drawdown engenders large pressure gradients close to the wellbore, creating conducive conditions for the generation of viscous forces that dominate existing capillary forces at the interparticle level (Alshmakhy and Maini, 2012). Gravity forces also compete with viscous and capillary forces (Løvoll et al., 2005; Rosado-Vazquez et al., 2007). Therefore, for the sustainability of improved fluid flow in this region, it is necessary that the viscous forces do not suppress contributions from gravity forces. Generally, gravity forces enhance oil recovery in areas that are not dominated by viscous forces (Rosado-Vazquez et al., 2007).

Drawdown was instigated by changing the pore pressure boundary conditions at the perforation region. Progressive reductions in pore pressures at the perforation region cause pore pressure gradients which are correspondingly larger within the vicinity of the wellbore. The greater pressure gradient results in considerable increases in the pore fluid velocity at the wellbore region. Figure 21 shows that, apart from an obvious rise at the near-borehole region, pore fluid velocity also increases with drawdown. It distinctly shows higher pore fluid velocities at the immediate surroundings of the perforation tunnel. This is further emphasised in Figures 22(a) and 22(b), where the pore fluid velocity, particularly around the perforation, increases when the drawdown is increased from 3·72 MPa (Figure 22(a)) to 10·34 MPa (Figure 22(b)). The rise in pore fluid velocity near the wellbore is attributed to high and disproportionate pore pressure gradients at the same zone. This facilitates the development of drag (viscous) forces which impact on the evolution of stresses within a close range of the wellbore. At this area, tensile radial stresses are induced as a result of fluid flow, and, when these surpass the rock tensile strength, there is tensile failure. It is also possible for compressive (shear) failure to occur if the tangential stress near the wellbore exceeds the shear strength of the rock; this is generally instigated by reservoir depletion or drawdown (Nauroy, 2011; Vaziri et al., 2002). Tensile and shear forces provoke material degradation and rock failure, which are essential requirements for sand production. An effective sand production management strategy can be achieved by manipulating drawdown, flow rate or pressure gradient.

The strain state of the rock is a central component of the sanding process. Like most engineering materials, rocks display a linear stress–strain relationship below the proportionality limit, the stress state beyond which the stress–strain relationship becomes non-linear. In certain instances, the material may still be within the elastic region if the yield point has not been attained and the strain is still recoverable. This implies that the proportionality limit and yield point are not necessarily synchronised, as there could be a slight difference between the two conditions. In this work, both stress–strain levels (i.e. proportionality limit and yield point) are in sync, giving way to the onset of plastic strain immediately after the elastic limit of the material is exceeded. The criterion for the onset of sanding is a function of the plastic strain. Similar to its stress counterparts, there are three basic types of strains: compressive, tensile and shear. The plastic regimes associated with any of these strains are usually co-opted for predictions of the onset of sand production. The equivalent plastic shear strain has been used to establish the onset of sand production (e.g. Mohamad-Hussein and Ni, 2018; Papamichos et al., 2001; Wang et al., 2005). This presupposes a material failure and sanding criterion based on the shear stress state, negating any tensile condition that may be imposed by the interplay of activities. Although the rock material failure may occur due to shear stresses, this does not necessarily signify the occurrence of sand production. It is possible for the radial stress near the wellbore to become negative due to imbalances caused by high pore pressure gradients which are, in comparison, greater in magnitude at this region (Hayavi and Abdideh, 2017a). Conditions for tensile failure could be satisfied if the negative radial stress is sufficiently high. Where this happens, the sanding process is dependent on the plastic shearing of the rock caused by compressive and shear stresses, and the tensile state induced by pore pressure gradient, fluid flow and drawdown.

An alternative to the equivalent plastic shear strain is hereby given in the form of the equivalent plastic strain (⁠$ϵ¯pl$⁠) (Equations 10 and 10b). It is the strain equivalent of the von Mises stress and enables an arbitrary 3D strain state to be denoted by a single value. As defined by its composition, it accounts for all categories of plastic strains. Plastic strain is shown to increase with drawdown (Figure 23). This is particularly pronounced within the vicinity of the perforation tunnel, extending outwards to a region of almost 1·5 m before tapering off to zero. The equivalent plastic strain $ϵ¯pl$ is a scalar variable that is used to indicate inelastic deformation and yield. A value of $ϵ¯pl$ greater than zero indicates material yield, and its magnitude shows the extent of plastic deformation. Figures 23(a) and 23(c) indicate a rise in $ϵ¯p$ as drawdown is increased from 3·72 to 10·34 MPa.

The sanding process is self-driven. The eroded area which is initiated at the vicinity of the wellbore and perforation tunnel enlarges due to the erosion process. This weakens the material, causing a redistribution of stresses to more intact areas situated further away from the wellbore, resulting in its softening, hence susceptibility to erosion. The extent of erosion therefore reduces as it progresses away from the wellbore vicinity; the maximum magnitude of erosion is at the perforation and wellbore. Figure 24 illustrates the start of plastic strains at the surface of the perforation tunnel at a constant drawdown of 3·72 MPa.

The distribution of pore pressure does not show significant differences with time. However, there is an apparent steepness of the pore gradient near the wellbore region, which is consistent for all simulation runs. It is worth noting that the larger pore fluid (interstitial) velocity that occurs close to the wellbore is attributed to the higher pore gradient, which also contributes to increases in stresses around this region, leading to a weakening and eventual erosion of the material. At a constant drawdown of 3·72 MPa, changes in vertical pressure do not cause significant changes in pore pressure distributions, as seen in Figure 25. In fact, the pore pressure variations match, following a similar trend. Further tests still need to be conducted to ascertain the influence of changes in drawdown conditions.

The minimum distance from the perforation tip beyond where there is zero equivalent plastic strain for a constant drawdown of 3·72 MPa is presented in Figure 26(a), which shows a flattening at the top, indicating that, after 1·2 d, plastic strains do not spread much further. An extension of Figure 26(a) is found in Figure 26(b), where a comparison is made at various drawdown conditions. The results are quite intriguing, displaying a progressive increase in the outer radial regions that are not plastic-strained with increasing constant drawdown. Despite the effect of drawdown, the occurrence and general trend is consistent with the literature, which asserts that the stress–strain conditions around wells cause the development of plastic regions at both the wellbore vicinity and a radius further away from the wellbore. The extent of the plastic radius is a function of the prevailing stress–strain conditions, among other factors. A similar pattern is observed (Figure 26(d)) when the vertical pressure is varied.

In contrast to expectations, the minimum radial distance away from the perforation tip after which no plastic strain occurs increases with reduced vertical pressure. An intriguing phenomenon was noticed when the vertical pressure was reduced to 55 MPa. From Figure 26(c), the minimum distance for a vertical pressure of 55 MPa reduced at 4·5 d to 15·7165 mm from a previous value of 39·2671 mm. This is an unusual occurrence because the distance is expected to recede with a lowering of the vertical pressure and not the contrary. To interpret this observation, it should be noted that, generally, the radial distance from the wellbore after which there is zero plastic strain increases with time. After a time period of 1·2 d, the maximum distance is attained, before becoming fairly constant (Figure 26(a)).

Apparently, there is a tendency for localisation of plastic strain around the wellbore. Plastic strain is mainly concentrated at the vicinity of the wellbore, which is indicative of the elevated level of activities in this area. Generally, after 4 d and as shown in Figure 26, regions at a distance beyond a 1 m radius from the perforation tunnel are still not strained above the proportionality limit/yield point.

On application of 55 MPa vertical pressure, pockets of areas with signs of no strains were noticed after only 1·2 d, becoming even more pronounced when the vertical pressure was reduced to 34·47 MPa (not shown in this paper). At this juncture, the authors can only speculate the phenomenon to be in agreement with the possibility of failure taking place at the far reach regions of the domain, leading to the creation of erosion channels known as ‘wormholes’ that eventually propagate to the wellbore region. Obviously, further work is required to confirm the authenticity of this claim, and, if proven, it will buttress the previous notion that the erosion process starts away from the wellbore vicinity and works its way towards the wellbore/perforation surface through wormholes.

It should be noted that the graphical plots only display the locations of zero plastic strain, and, although the differences in values of distance (due to differences in drawdown and vertical pressures) are relatively small, it is expected to become significant when the duration is increased.

As stated earlier, the equivalent plastic strain parameter was adopted as the criterion guiding the onset of erosion, whereby ablation of the material will occur if the value of the equivalent plastic strain, $ϵ¯pl$⁠, exceeds a cut-off value above zero. Above this, the material is assumed to have yielded and the actual sanding process signified by the detachment of rock particles can take place only if $ϵ¯pl$ reaches or exceeds a predetermined value referred hereafter as the ‘erosion or sanding criterion’ or ‘cut-off equivalent plastic strain’. For the preliminary study, values for $ϵ¯pl$ were chosen arbitrarily so as to enable a more pronounced evidence of erosion. Vicissitudes in erosion criterion reveal an inverse relationship, with lower values resulting in greater sand production. Thus, in Figure 27, significant changes occur when the cut-off is reduced to 0·01. The importance of determining an accurate and more realistic criterion value is not underplayed here. This will be a subject for future work as it invariably entails laboratory experimentation. The sand production coefficient (λ(γ ^p)) is a function of the sanding criterion and evolves in a manner reliant on the prevailing state of plastic strain. Plastic strains that develop are constantly juxtaposed with the cut-off equivalent plastic strain, and the updated value of the sand production coefficient relies on whether the cut-off plastic strain is exceeded or not.

The efficacy of the pressure applied on the wellbore/perforation face in the erosion process was tested by varying the pressures, all other conditions remaining the same. The applied pressures are considered representative of the well operation ‘mud pressure’ normally applied to maintain the integrity of the wellbore during drilling and production phases. The values used were arbitrarily selected and the outcome, as presented in Figure 28, portrays a drastic reduction in the quantity of sand eroded with increasing mud pressure. Nevertheless, this effect becomes relatively less significant after the mud pressure is increased above 37·2 MPa. Beyond this value, the amount of pressure applied seems to have very little effect, leading to the conclusion that the optimal mud pressure under the prevailing well operation condition is approximately 37·2 MPa. At higher magnitudes of mud pressures, negligible reductions in the quantity of sand produced are observed.

Mud pressure (alternatively known as mud weight or mud density) is one of the key operational parameters that determine well stability during drilling and production (Eshiet and Sheng, 2018; Hayavi and Abdideh, 2017a, 2017b; Millikan, 1938). Oil well control measures usually include aspects of mud pressure management to prevent influx of ‘kick’ (formation fluid) into the well during drilling. Wellbore instability can be prevented if an appropriate magnitude of mud pressure is maintained such that the integrity of the well is sustained throughout drilling and operation stages. Bottom-hole pressure (BHP) is a derivative of mud pressure. At the well bottom, mud pressure is denoted by a function expressed as the BHP. It is imperative to apply the right magnitude of mud pressure to institute a balance against the stress concentration at the wellbore wall generated by in situ stress regimes, pore pressure and forces instigated by reservoir fluid flow. Mud pressure is instrumental to the prevention of wellbore failures (Hayavi and Abdideh, 2017b) as well as other more severe phenomena, including wellbore collapse and breakouts (Gholami et al., 2016). In overbalanced conditions, the penetration of mud pressure increases the pore pressure around the wellbore, which, in turn, reduces the effective support that would have been otherwise provided by the mud pressure. High pore pressure and pore pressure gradients encourage sand production (Wu et al., 2006b). A link can therefore be drawn between applied mud pressure and sand production. Wellbore stability is enhanced through the application of mud pressure, with an inverse consequence on the propensity for reservoir erosion. An optimum or, at least, a safe mud pressure window is one of the strategic prerequisites for minimising the sanding process.

An extended form of parametric study was conducted by examining the influence of certain reservoir and well production parameters on the failure and subsequent erosion of the formation material using FEM modelling techniques. The major factors considered include drawdown, wellbore/perforation depth and erosion criterion. In addition, the optimal mud pressure was determined. The adoption of methodologies based on FEM necessitated the conduct of sensitivity analyses and validations of predicted results. The mesh size distribution employed guaranteed high resolution and accuracy of results. A mesh density of 9507 elements in the adaptive zone (well face and perforation region) was used, which falls within the range of convergence. An increase in the mesh density of this region above a threshold value led to disparity (and non-convergence) in results when the mesh density of the outer rock domain was kept constant. Selection of appropriate mesh densities was based, hence, on the initial convergence test. There is a qualitative resemblance in results between the FEM model and models by other researchers (Papamichos and Stavropoulou, 1998; Papamichos and Vardoulakis, 2005; Papamichos et al., 2001). The discrepancy in the magnitude of the variables is attributed to differences in initial, boundary and operating conditions; geometry/configuration; and size. Results from literature are based on laboratory experiments, which are not only smaller in scale in comparison to the developed FEM model, but are also slightly different in configuration.

Variations in wellbore depth indicate an increase in sand production with wellbore/perforation depth and a corresponding decrease in values of higher fluid flow velocities typically observed at the wellbore/perforation zone. Changes in constant drawdown conditions show an increase in sand production with drawdown. Sand production was considerably high at the drawdown pressure of 10·34 MPa. There is also an increase in plastic strain with increasing drawdown. Similarly, variations in the magnitude of pressure applied on the wellbore/perforation face, which represents the applied mud pressure, indicate a significant reduction in the severity of sand production when mud pressure is raised, with an optimal value occurring when the mud pressure is increased to 37·2 MPa. Above this value, the changes are negligible. Measurements of the minimum distances from the perforation tip after which there are no plastic strains show that subsequent to the initial occurrence of plastic strain, the plastic zone increases at a steady rate prior to becoming approximately constant at the posterior of a given period. Furthermore, there is a progressive increase of the plastic region with decreasing constant drawdown. The same effect is observed when the wellbore/perforation depth is reduced. The results obtained show trends which enhance further understanding of the sanding phenomenon.

This research work was adopted from a section of the doctoral thesis on reservoir erosion presented by Eshiet (2012).