#### Abstract

Investigating the shear failure caused by the concentration of compressive stress around noncircular boreholes is important both in the field and in the laboratory. This article deals with the numerical analysis of elliptical boreholes under a nonisotropic in situ stress field using the Mogi–Coulomb nonlinear failure criterion. The purpose of the presented numerical model is to simulate the progressive shear failure (breakout) around the borehole and investigate the impact of the eccentricity of the borehole on the stability and depth and width of the failure area. According to the obtained results, the breakout is V-shaped and is formed along the minimum principal stress. As the eccentricity of the borehole increases, the final dimension of the breakout becomes smaller; in other words, the increase in ellipticity strengthens the borehole against shear failure. However, as the eccentricity increases, the stress concentration at the breakout tip increases. Another finding of the study conducted in this article is the significant relationship between the width and the depth of the breakout failure, which makes the idea of estimating both horizontal in situ stresses using breakout dimensions seriously doubtful. Also, the interesting result obtained is that the stress concentration factor at the breakout tip for boreholes with different eccentricities is the same at the end of the breakout.

#### 1. Introduction

Often 10–20% of the total drilling cost is related to borehole instability, this causes a loss of about 500 to 1000 million dollars in the oil industry worldwide. Therefore, maintaining well stability is one of the main concerns of drilling engineers in the oil and gas industry [1–4]. One of the instabilities of the borehole is related to the shear instability of the rock in the borehole wall, which is known as the borehole breakout.

When the shear stress of the material reaches its shear strength in the borehole wall, the rock is crushed and falls into the borehole. The main driver of this type of failure is the concentration of compressive stress that occurs due to the drilling of the borehole. This phenomenon was first observed and reported by Cox [5] in Alberta wells and later confirmed by Babcock [6]. As the degree of in situ stress anisotropy increases, the intensity of stress concentration increases, and as a result, the failure zone becomes wider and deeper [7–14].

Breakout is formed symmetrically on both sides of the borehole along the minor principal in situ stress where the compressive stress is most concentrated. For this reason, in the last two decades, breakout has been used as an indicator to determine the direction and magnitude of in situ stresses [15–17]. Also, several laboratory studies conducted on predrilled rock samples show that the size of the borehole (hole radii) has a significant effect on the breakout initiation stress [9, 18–21].

In addition to physical models and laboratory studies on borehole breakouts, various theoretical and numerical models have also been conducted to investigate breakouts. Although the first theoretical models presented for the breakout had simple assumptions, the episodic progression of the breakout was not investigated in them [7, 22], but in the following, more complete numerical models based on the finite element method and the discrete element method were presented, which are able to predict the breakout in a proper way [23–32].

Setiawan and Zimmerman [33] presented a semianalytical method to investigate breakout progression and stabilization using the conformal mapping procedure, and the analysis carried out by them led to a correlation between the breakout geometry, the properties of the rock materials, and the in situ stresses.

Zhang et al. [34] investigated borehole breakout using a finite element model based on the Mogi–Coulomb failure criterion, and the results of their analysis showed that the inverse analysis using the finite element model and neural network can effectively determine the in situ stresses.

Li et al. [35] presented a hydromechanical finite element model for sand production and erosion and observed that sand production is mainly controlled by the plastic strain magnitude and flow velocity in the vicinity of the borehole.

Ma et al. [36] presented a finite element numerical simulation method based on elastic damage mechanics for progressive sand production in inhomogeneous formations, and the results obtained by them showed that the sand production area is controlled by the rock formation heterogeneity and the expansive failure process.

In addition to various studies conducted on circular boreholes, several studies have also been conducted on noncircular and elliptical boreholes. By investigating the tangential stresses around the elliptical boreholes, Aadnoy and Angell-Olsen [37] observed that the fracture initiation pressure in the elliptical borehole is different from the circular borehole, but as long as the ratio of in situ stresses is not greater than ellipticity, the position of fracture in the elliptical and circular borehole is the same, and in their study, ellipticity is the ratio of the small diameter of the ellipse to its larger diameter.

Aadnoy and Kaarstad [38] presented an elliptic geometry model to investigate sand production, and they observed that the anisotropy of the in situ stresses is a critical factor for the elliptic shape of the borehole.

Using the ellipsoidal solution of stresses acting on the borehole wall and the Mohr–Coulomb failure criterion, Aadnoy et al. [39] found that the ellipsoidal geometry presents the collapsed borehole shape better than the circular Kirsch’s analytical solution.

Using an inverse analysis, Han et al. [40] estimated the horizontal in situ stresses around the elliptical borehole using data from the leak-off test, and their studies showed that even a 2% difference in the axis of the elliptical borehole causes a 5% to 10% difference in the estimation of horizontal in situ stresses. Papamichos et al. [41] investigated hollow cylinders with different geometries experimentally and numerically, and their studies showed that holes with elliptical breakouts are more stable; in other words, the stability of the hole increases as the breakout depth increases.

What has not been addressed so far is the investigation of the impact of the eccentricity of the borehole on the stability of the borehole as well as its strength to shear failure. The purpose of this article is to present a model based on the finite element numerical method to investigate the episodic progression of breakouts around noncircular boreholes and specifically boreholes with an elliptical cross section. The presented model is two-dimensional and assumes plane strain conditions, and by it, the effect of anisotropic in situ stresses and borehole eccentricity on the final dimensions of the breakout is investigated.

#### 2. Problem Definition

There are various reasons that boreholes with noncircular cross sections may also be formed. Possible causes include the mechanical action of the drill string on the well after drilling or the horizontal cross section of the deviated boreholes becomes an ellipse.

An elliptical cavity with major axis *a* and minor axis *b* is considered in a linear elastic medium under plane strain conditions (Figure 1). is the maximum horizontal principal stress, is the minimum horizontal principal stress, and is the vertical stress along the axis of the borehole so that . It is often assumed that the major axis of the ellipse is along the minor principal stress and its minor axis is along the major principal stress . This assumption is suitable for examining holes that have previously had shear failure and are again subjected to a new stress condition.

In these conditions, local stresses including tangential stress , radial stress , vertical stress and shear stress are obtained for the points around the borehole, and it is observed that the highest concentration of compressive stress occurs in the wall of the borehole at point A because this point has the minimum radius of curvature. However, in the vicinity of point A, the concentration of compressive stress may also occur in other points and they may suffer shear failure (breakout) as a result. To find all the failure points, it is necessary to first calculate the principal stresses based on the local stresses, and then by choosing an appropriate shear failure criterion, the points that have failed can be found.

The failed points are removed from the borehole wall, and a new cross section is obtained. The stresses are redistributed around the borehole, which causes other points to reach shear failure. Again, the points that have suffered shear failure are removed and the new stresses in the environment are calculated. In this way, step by step, the breakout proceeds so that no point will fail again.

As shown in Figure 1, in breakout, the rock is separated from the borehole wall in spalls; although from a micromechanical point of view, the failure mechanism may be a combination of shear failure and tension failure.

In all breakout stages, to calculate the stresses of the environment, the environment must first be divided into fine elements, and then, using the appropriate numerical method, the stresses in the center of each element should be calculated. However, only in the wall of the borehole before the initiation of the breakout can the stresses be obtained using the analytical method that is given in the next section.

After completing the numerical analysis, the V-shaped (dog-eared) breakout is observed on both sides of the borehole symmetrically along the minimum principal stress.

The maximum breakout width is the same width that is obtained in the first step of failure, and the final depth of the breakout is obtained after the last step at the end of the analysis. It can be seen in Figure 1 that is measured from the center of the borehole.

##### 2.1. Analytical Stress Distribution around the Elliptical Borehole

In this section, the stress distribution around the elliptical holes is presented [42]:where

and are also shown in Figure 2.

##### 2.2. Numerical Model for the Episodic Breakout

A two-dimensional plane strain model is used to analyse the progressive breakout phenomenon. The presented model is solved numerically using the finite element numerical method.

The finite element method, sometimes called finite element analysis, is a computational technique used to obtain approximate solutions of boundary value problems. Boundary value problems are also sometimes referred to as field problems. The field is the domain of interest and often represents a physical structure. The field variables are the dependent variables of interest that are governed by the differential equation. Boundary conditions are the specified values of field variables at field boundaries. The finite element method is a numerical technique for solving a system of equations governing the domain of a continuous physical system, which is discretized into simple geometric forms called finite elements. Modeling a body is conducted by dividing it into an equivalent system of finite elements that are connected at a finite number of points on each element called nodes. There are fundamental unknowns in engineering problems, and if they are found, the behavior of the entire system is predictable. The basic unknowns or field variables encountered in engineering problems are displacements. In a continuum body, these unknowns are infinite. The finite element method reduces such unknowns to a finite number by dividing the solution area into small parts called elements and by expressing the unknown field variables in terms of hypothetical approximate functions (interpolation functions/shape functions) in each element. Approximate functions are defined in terms of field variables of specified points called nodal points. Thus, in the finite element method, the unknowns are the field variables of the nodes. Once these are found, the field variables at each node can be found using the interpolation functions. After selecting elements, the next step in the finite element method is to collect element properties for each element; in other words, the stiffness characteristics of each element must be found. Mathematically, this relationship is as follows [43]:where is the element stiffness matrix, is the gradient matrix, and is the element thickness. implies the elastic stiffness matrix which in plane strain conditions is given bywhere *E* and are the modulus of elasticity and Poisson’s ratio, respectively.

Derivation of the element stiffness matrix is based on equilibrium conditions. The same procedure can be applied by writing the equilibrium equation for each node for all connected elements in the model. This process is described as “assembly” because the system equations are obtained by taking the individual stiffness components and putting them together; therefore, the main relation is written as follows:where is the global stiffness matrix, is the nodal displacement vector, and is the nodal force vector.

The criterion of shear failure used in this article is the nonlinear Mogi–Coulomb criterion, in which the effect of the intermediate principal stress is also present in the failure function [44]:

In this equation,where , , and are the maximum, intermediate, and minimum principal stresses, respectively, and where is the uniaxial compressive strength of the rock material. The failure function for the Mogi–Coulomb failure criterion is defined as follows [45]:

Figure 3 shows the Mogi–Coulomb failure criterion in the principal stress space.

##### 2.3. Algorithm, Problem Domain, Boundary Conditions, and Meshing

To simulate the episodic breakout around the elliptical borehole, a computer program was coded based on the governing equations and failure criteria presented in the previous section.

Coding is conducted in MATLAB software. To increase the accuracy of the simulation, every element that fails is removed from the borehole wall by 0.1 of its length. In this way, the number of numerical analysis iterations increases until the breakout reaches stability but the breakout geometry is obtained with proper accuracy. The problem domain, boundary conditions, and meshing shown in Figure 4 are considered, in which only a quarter of the model is analyzed due to symmetry.

The dimensions of the model are and a fine mesh has been used around the borehole. The number of 4-node quadrilateral elements for the model is equal to 2560. According to Figure 4, the bottom mesh boundary is restricted for vertical displacement while it is left free for horizontal displacement. Also, the right vertical boundary is free to move vertically while its horizontal movement is restricted. The maximum horizontal principal stress is applied to the right vertical boundary, and the minimum horizontal principal stress is applied to the upper horizontal boundary. For the numerical solution of the problem shown in Figure 4 and the specifications given previously, the following simple algorithm is presented:(1)First, based on the finite element method formulation, local stresses are calculated for each element(2)Then, the principal stresses are calculated according to the local stresses(3)Using the principal stresses, the failure function is calculated for each element according to the selected failure criteria(4)If the failure function for a certain element is greater than one , that element is removed from the model geometry and the rest of the elements remain(5)For the new geometry, steps 1 to 4 are repeated until the failure function for all elements becomes smaller than one and the breakout reaches stability

#### 3. Validation

##### 3.1. Stress around the Borehole and Comparison of the Numerical Method with Analytical Solution

To validate the presented numerical method, the stresses around the borehole have been obtained using the numerical method as well as the analytical relations presented in the section “Analytical stress distribution around elliptical borehole.” The comparison of the stresses is shown in Figures 5–8. Stresses are presented for three different eccentricities, and the eccentricity of an ellipse is equal towhere according to Figure 1, a and b are the major and minor axes of the ellipse, respectively.

*m* is between zero and one , such that the circular borehole has zero eccentricity. In this article, to create boreholes with different eccentricities, the major axis of the ellipse is assumed to be fixed and its minor axis is changed (Figure 1).

Figure 5 shows the tangential stress in the borehole wall for different values of eccentricity, and as can be seen in the figure, there is a good agreement between the stresses obtained from the numerical model and the stresses obtained from the analytical relations.

It should be noted that in Figures 5–8, the values of the field stresses are as follows:

Figure 5 also shows that the maximum tangential stress or the highest compressive stress concentration occurs for and , that is, along the minor principal stress; therefore, these two points will be the starting points of the breakout. It can also be seen that with the increase of eccentricity, the concentration of compressive stress also increases. The lowest tangential stress also occurs at and , and with the increase of eccentricity, the tangential stress decreases at these points. In Figure 6, similar results have been obtained for vertical stress. Figures 6 and 7 also show the stresses for and along the minor principal stress for 0 and 0.6 eccentricity, respectively, there is a good agreement between the numerical model and the analytical relationships, and by moving away from the center of the hole, the stresses converge to the field stresses.

##### 3.2. Failure Shape and Comparison of the Numerical Method with Other Models

Ma et al. [36] presented a numerical simulation method for sand production in inhomogeneous formations, in which the effects of heterogeneity, progressive fracture process, and borehole pressure were investigated. Zhang et al. [34] also used the machine learning method to invert the relationship between in situ stresses and borehole breakout shape. Figure 9 shows the comparison between the depth and the width of the breakout failure zone obtained from the numerical model presented in this article with the model provided by Zhang et al. [34] and Ma et al. [36] for circular borehole . The analysis of this section has been performed for the mechanical and geometrical specifications presented in Table 1. According to Figure 9, the failure width is formed in the first iteration and does not increase in subsequent iterations. The half width of failure obtained from the numerical model presented in this article is and the half widths of failure obtained from Zhang et al.’s model and Ma et al.’s model are and respectively, which can be seen that there is a relatively good agreement between the models.

However, in all three models, the depth of failure increases with increasing iterations until a stable state is established. The evolution of is linear in the model presented in this article, but it is nonlinear in other models. This is because in the numerical model in this article, to increase the accuracy of the results with the failure of each element as much as one tenth of the element length, the borehole wall is removed, but in other models, the failed element is completely removed from the model. obtained from the numerical method presented in this article is and obtained from Zhang et al.’s model and Ma et al.’s model is and respectively. The difference between the values obtained from different models can be due to the different failure criteria used in each of the models. The failure criterion used in this article is the Mogi–Coulomb criterion, while the failure criteria used in Zhang et al.’s model and Ma et al.’s model are the Mohr–Coulomb criterion and the Drucker–Prager criterion, respectively.

#### 4. Results and Discussion

In this section, the results of the numerical analysis performed using the program written in MATLAB software are given. The properties selected in the analysis for the rock materials are those provided for Tablerock sandstone. The Tablerock sandstone belongs to the Cloverdale Nursery area in the United States and belongs to a group of sandstone layers in the Upper Miocene, the Lower Idaho Group; this sandstone consists of 55% quartz and 37% feldspar. Due to the high percentage of feldspar, this sandstone is classified as arkosic. For this type of rock material, the internal friction angle and cohesion are, respectively, equal to and [13].

For the borehole with different eccentricities, breakout progression steps were obtained using the presented model and some of these stages are shown in Figure 10. The cross section shown in the last column is the final breakout shape after reaching stability. The values of in situ stresses in this analysis are equal to .

In the figure, it can be seen that for the borehole with zero eccentricity, the breakout reaches stability after 383 iterations, and as the eccentricity increases, the number of breakout progress iterations decreases. For example, for the borehole with an eccentricity of 0.8, the number of iterations is reduced to 223. It can be seen that the final shape of the breakout is V-shaped, and its formation and propagation occurred along the minor principal stress, and also, with the increase in the eccentricity of the borehole, the final dimensions of the breakout became smaller.

Figure 11 shows the increase in the depth and width of the breakout failure area as the breakout progressed to stability. As can be seen in the figure, the failure width is formed in the first iteration and does not increase in subsequent iterations, but the depth of failure increases linearly until the breakout reaches stability.

Figure 12 shows the final normalized breakout depth versus eccentricity, and Figure 13 shows the breakout failure width versus eccentricity. In Figure 12, the breakout depth is normalized to the semimajor axis of the ellipse. It can be seen in the two figures that the dimensions of the breakout become smaller as the eccentricity increases. The general result is that although the stress concentration at the two vertices of the ellipse along the minimum principal stress increases with the increase of borehole eccentricity, the final dimension of the breakout becomes smaller. In other words, increasing the ovality strengthens the borehole against shear failure.

Another influencing parameter on the dimensions of the breakout is the ratio of in situ stresses. Numerous numerical analyses were performed with different in situ stress ratios and eccentricities, and their results are shown in Figures 14–16. In Figures 14 and 15, it can be seen that with the increase in the ratio of in situ stresses, the normalized breakout depth and the breakout width increase, respectively. Also, for a certain in situ stress ratio, the depth and width of failure decrease with the increase of eccentricity. Figures 14 and 15 can also be used to determine the final dimensions of the breakout after increasing the ratio of in situ stresses to another given ratio. This application is of particular importance in breakout laboratory studies.

Figure 16 is obtained from the combination of Figures 14 and 15. Figure 16 shows a significant relationship between the final breakout failure depth and the breakout failure width for a given eccentricity. The meaningful relationship between the depth and the width of breakout failure in geomechanics is important because if the horizontal in situ stresses are to be determined based on the dimensions of the breakout, only one of them can be determined. In other words, the minimum principal stress can be obtained from the hydraulic fracturing test [46, 47], and the maximum horizontal principal stress can be determined based on the dimensions of the breakout.

The stress concentration factor (SCF) is the ratio of the tangential stress at the breakout tip to the maximum in situ stress, which is as follows:where is the tangential stress at point A at the breakout tip (Figure 1). Point A also moves as the breakout progresses.

Figure 17 shows the evolution of the stress concentration factor (SCF) versus the breakout depth for three different eccentricities. The highest concentration of compressive stress occurs in the borehole wall and at the tip of the breakout, which is the cause of shear failure progression.

The remarkable point that can be concluded from Figure 17 is that, before the breakout failure and at the beginning of the analysis, the greater the eccentricity, the greater the stress concentration factor, but with the progress of breakout failure and at the end of the analysis, it can be seen that the stress concentration factor for boreholes with different eccentricities has become the same with a relatively small difference. This is due to the fact that in the studied problem, the tangential stress in the borehole wall is the only principal stress in the Mogi–Coulomb failure criterion, and the breakout failure progresses to the point where the tangential stress can no longer satisfy the failure criterion. Therefore, the stress concentration factor at the breakout tip, which is dependent on the tangential stress, at the end of the breakout, is the same for the boreholes with any eccentricity, even though the breakout depth is different.

#### 5. Conclusion

The work performed in this article was the numerical analysis of the progressive failure around noncircular boreholes. This problem, which has an unconventional shape, was investigated in the field of nonisotropic stresses and using a nonlinear failure criterion. For these reasons, solving the problem had to be performed numerically. In addition to the abovementioned, the progress of failure occurs episodically, which means that at each step, the geometry of the problem changes and meshing must be conducted again. To solve this problem, a simple algorithm based on the finite element method is presented. Using the algorithm, a computer program is coded in MATLAB software.

To investigate the effect of the borehole eccentricity parameter on the breakout, boreholes with different eccentricities are analyzed, and it was observed that with the increase of eccentricity, the depth and width of the V-shaped breakout failure area decrease, but the stress concentration factor at the breakout tip becomes the same for all models. Also, the rate of decrease in the breakout depth and width with increasing eccentricity is low until *m* = 0.3 and then increases. In other words, eccentricity strengthens the borehole against shear failure.

With the increase of the ratio of in situ stresses, the depth and width of the breakout increase and the rate of increase of the failure depth is higher for the borehole with smaller eccentricity. Also, according to the obtained results, it was observed that the depth and width of the breakout have a significant relationship with each other. Therefore, by using breakout dimensions, only one of the in situ stresses can be obtained; for example, the maximum horizontal stress can be obtained using breakout dimensions, and the minimum horizontal stress can be obtained using the hydraulic fracturing test.

#### Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.