Erosion of material (sand production) in an oil wellbore

The process of optimizing the production value of an oil well is complex but can be simplified as a balance between the oil recovery rate (as measured by the volumetric flow rate of oil), the sustainability of the recovery (as measured by the amount of oil ultimately recovered), and the cost of operating the well. In practice, achieving this balance requires consideration of the erosion of rock in the wellbore, a phenomenon that occurs when oil is extracted under a sufficiently high pressure gradient. This erosion phenomenon is generally referred to as “sand production.” Depending on the flow velocities the sand may accumulate in the well, affecting the sustainability of recovery, or it may get carried to the surface along with the oil. The sand in the oil causes erosion in the piping system and its components such as chokes and pipe bends, increasing the costs of operating the well. Excessive sand production is, therefore, undesirable and will limit oil recovery rates. A typical measure of these limits on recovery rates is the so-called “sandfree rate.” A sandfree rate might be based on direct damage caused by sand production or might be based on the cost of sand management systems that limit the damage of the sand. The former measure is called the maximum sandfree rate, or MSR. The latter measure is called the acceptable sandfree rate, or ASR. The ASR concept has become possible with the availability of many commercial sand management systems as well as new designs of piping components. The ASR concept has also engendered a need for predicting the sand production rates to properly choose and size the sand management systems and piping components. In this example we focus on measures that can be obtained from Abaqus to predict these rates.

The geometry of an oil well has two main components. The first is the wellbore drilled through the rock. The second component is a series of perforation tunnels that project perpendicular to the wellbore axis. These tunnels, which are formed by explosive shape charges, effectively increase the surface area of the wellbore for oil extraction. The perforation tunnels are typically arranged to fan out in a helical fashion around the wellbore, uniformly offset in both vertical spacing and azimuthal angle.

A linear Drucker-Prager model with hardening is chosen for the rock, and the casing in the three-dimensional model is linear elastic.

The loading sequence for a wellbore analysis generally includes

• establishing geostatic equilibrium, based on the overburden loading;

• simulation of material removal operations, including drilling the wellbore and forming the perforation tunnel; and

• applying a drawdown pressure in the wellbore to simulate pumping.

This sequence is modeled slightly differently in the three-dimensional and planar models.

There are two main sources of eroded material in a wellbore. One of the sources is volumetric and is due to the material that is broken up due to high stresses and transported by the fluid through the pores. The other source is surface based and is due to the material that is broken up by the hydrodynamic action of the flow on the surface. Depending on the properties of the oil-bearing strata and the flow velocities, one or the other may be the dominant source of eroded material. Development of equations describing the erosion in a wellbore is an active research field. In this example we consider surface erosion only but choose a form for the erosion equation that has dependencies that are similar to those used by Papamichos and Stavropoulou for volumetric erosion. This approximation is reasonable because high stresses exist only in a very thin layer surrounding the wellbore. The erosion equation is

where ${\mathrm{V}}_e$ is the erosion velocity, $v_w$ is the pore fluid velocity, c is the transport concentration, n is the porosity, and $\lambda$ is the so-called sand production coefficient. $\lambda$ depends on the equivalent plastic strain (${\overline{\epsilon }}^{pl}$). It is zero below a cutoff equivalent plastic strain ($\overline{\epsilon }^{pl}{}_{\mathrm{cutoff}}$), is equal to ${\lambda }_1\left({\overline{\epsilon }}^{pl}-\overline{\epsilon }^{pl}{}_{\mathrm{cutoff}}\right)$, and is limited by a constant ${\lambda }_2$.

Both ${\lambda }_1$ and ${\lambda }_2$ must be determined experimentally. In this example $\overline{\epsilon }^{pl}{}_{\mathrm{cutoff}}=0.028$, ${\lambda }_1=4$, and ${\lambda }_2=0.01$. We choose $c=0.001$ as recommended by Papamichos and Stavropoulo. These values are chosen to show visible erosion in a reasonably short analysis time.

Erosion is modeled during the final step of each analysis. The erosion equation describes the velocity of material recession as a local function of solution quantities. Abaqus/Standard provides functionality through adaptive meshing for imposing this surface velocity, maintaining its progression normal to the surface as the surface moves, and adjusting subsurface nodes to account for large amounts of erosive material loss.

Erosion itself is described through a spatial adaptive mesh constraint, which is applied to all the nodes on the surface of the perforation tunnel. Adaptive mesh constraints can be applied only on adaptive mesh domains; in this example a sufficiently large extent of the finite element mesh near the wellbore and perforation tunnel surfaces is declared as the adaptive mesh domain. A cut section of the adaptive mesh domain for the three-dimensional model, including the perforation tunnel, is shown in Figure 2. The adaptive mesh domain for the planar model is the regular mesh near the wellbore (refer to Figure 6 and Figure 8). Identification of the adaptive mesh domain will result in smoothing of the near-surface mesh that is necessary to enable erosion to progress to arbitrary depths. All the nodes on the boundary of the adaptive mesh domain where it meets the regular mesh must be considered as Lagrangian to respect the adjacent nonadaptive elements.

A velocity adaptive mesh constraint is defined. The generality and solution dependence of the erosion equation are handled by describing the erosion equation in user subroutine UMESHMOTION (see Defining ALE adaptive mesh domains in Abaqus/Standard). User subroutine UMESHMOTION is called at a given node for every mesh smoothing sweep. Mesh velocities computed by the Abaqus/Standard meshing algorithm for that node are passed into UMESHMOTION, which modifies them to account for the erosion velocities computed at that node. The modified velocities are determined according to the equation for ${\mathrm{V}}_e$, where local results are needed for ${\overline{\epsilon }}^{pl}$, n, and $v_w$. To obtain these values, we request results for output variables PEEQ, VOIDR, and FLVEL respectively, noting that void ratio is related to porosity by $\left(1-n\right)=1/\left(1+\mathrm{VOIDR}\right)$. Since these output variables are all available at element material points, the utility routine GETVRMAVGATNODE is used to obtain results extrapolated to the surface nodes (see Obtaining material point information averaged at a node).

Rezoning, the process of creating a new mesh in the model's deformed configuration, is a useful technique in addressing element distortion in erosion problems.

The erosion model in this example is solution dependent in the sense that it defines an erosion velocity as a function of local values of equivalent plastic strain, PEEQ, and fluid velocity, FLVEL. Since these node-located results values are extrapolated from adjacent elements' material point results, using utility routine GETVRMAVGATNODE, the quality of the results used in the erosion model is dependent on local element quality. In practice, as elements near the surface deform, there is a tendency for instability in the erosion model.

This instability is mitigated by rezoning the model and providing a more regular mesh to continue the analysis. The rezoned model is created in Abaqus/CAE, and transfer of state variables occurs using mesh-to-mesh solution mapping in Abaqus/Standard. Rezoning occurs at 55,000 seconds (15.3 hours) into the 32-hour erosion analysis of the planar model.