© N.A. Batalha et al., published by IFP Energies nouvelles, 2020

1 Introduction

In order to predict the formation behavior in drilling operations, numerical models have been developed to simulate the stress state around vertical and inclined wellbores [1–10]. These models intend to define the best drilling fluid density (or mud weight) to be used. The mud weight must be designed to keep the well stable and not cause rock failure. In-situ stresses magnitude and orientation, as well as rock strength and pore pressure, are some of the parameters required to a geomechanical approach in wellbore stability, as presented by [11]. The precision and variation of these parameters over a domain are generally not considered, but average values are usually obtained from multiple laboratories or on-site tests. Because of that, most models assume that the properties of the material are constant in the entire simulated domain. Thus, wellbore stability analysis may involve many uncertainties; to list a few: definition of formation’s mechanical properties, pore pressure, mud weight, system temperature, and in-situ stresses. Therefore, analyzing the influence of these variables requires probabilistic methods and a consistent framework.

A probabilistic analysis considering the uncertainties of in-situ stresses, pore pressure, and rock strength as input data in a geomechanical model was evaluated by Udegbunam et al. [12]. Monte Carlo simulations were used, and triangular and constant distribution types were applied for the random variable distributions in order to predict critical fracturing and collapse pressures. It was observed in this study that the minimum horizontal stress has a high impact on the cumulative uncertainty in the fracture pressure prediction. Cohesive rock strength was pointed as the variable that most influences the collapse pressure.

Tran and Abousleiman [13] studied how the mechanical properties, in-situ stresses, and formation pore pressure conditions affect wellbore stability when using Mohr-Coulomb, Drucker-Prager, Modified Lade, and tensile failure criterias. Monte Carlo method was used in this study, and different far-field regimes were also evaluated. Regarding the mechanical properties, cohesion and friction angle presented the most impact on collapse mud weight. It was noted in this work that the in-situ stresses and formation pore pressure had more influence in tensile failure than rock tensile strength.

Fontoura et al. [14] analyzed the influence of the uncertainties in wellbore stability using Statistical Error Analysis Method (SEAM), First Order Second Moment (FOSM) and First Order Reliability Model (FORM) methods, which were compared with Monte Carlo method. Based on a sensitivity analysis, it was noted that formation cohesion and angle of internal friction are essential parameters for collapse failure, and, on the other hand, tensile strength is relevant for fracturing failure. Also, in-situ stresses and formation pore pressure are considered important parameters in wellbore stability analysis, as well as Biot’s constant. In the case of vertical wellbore stability, the vertical stress was not relevant.

A numerical model for wellbore stability coupled to hydromechanical process was presented by Muller et al. [15], in which stochastic and deterministic analyses were compared; the spatial variability of initial stresses, pore pressure, and hydraulic and mechanical properties were taken into account. A reliability analysis was carried out limited to a pressure range, where lower limit, compressive failure, is based on a user defined plastic region, while the upper limit is based on the rock tensile failure stress. This study was later analyzed for three-dimensional scenarios [16].

Quantitative risk assessment has also been used in wellbore stability uncertainties analysis [11, 17, 18]. The latter reference presented an analytical approach for wellbore stability uncertainties analysis, which takes into account the uncertainties of in-situ stresses, operational activities, and rock mechanical properties. According to the author, their results have shown that cohesion, internal friction, and pump pressure were not significant parameters in the total risk of the drilling operation.

Probabilistic models have led to high impact in wellbore stability analysis, especially when sensitivity analysis is taken into account to define input parameters. The influence of elastic modulus heterogeneity on in-situ stresses and damage in shale gas reservoirs was recently studied by Wang et al. [19], in which variability of von Mises stress, strain energy density, and damage factor was observed as heterogeneity increased.

Therefore, to obtain a more realistic wellbore stability analysis, a heterogeneous field shall be considered. Since the variation of rock mechanical properties along the field is seldom known due to difficulties to precisely define these properties, a random but correlated field can be assumed to carry out a probabilistic analysis. The present work introduces a geomechanical model coupled to a geostatistics algorithm to predict heterogeneous formation behavior in drilling operations.

The geomechanical model presented herein assumes linear elastic and isotropic mechanical properties and plane-strain state. A yield criterion is selected to identify the formation’s failure area, plastic region, around the borehole wall. The stochastic field is based on a covariance function, which specifies the correlation between the mechanical characteristics of neighboring elements of a domain. This stochastic analysis approach enables a more realistic evaluation of the wellbore failure; the model is also extended to analyze vertical, horizontal, and inclined wellbore.

2 Rock and drilling fluid: a mechanical interaction

A stable and well-designed wellbore guarantees operation safety and prevents several issues such as stuck pipe, cavity enlargement, fracturing formation, operation delays, and substantial expenditures. Note that the stability of a wellbore depends on equilibrium of stresses (in-situ stresses), material properties, and drilling fluid’s density, also called mud weight.

In-situ stresses are naturally developed during a rock formation, and they are a consequence of gravity interactions, tectonic processes, and others. These stresses are composed by the overburden pressure σ_ov, defined by the weight of all rock layers over a rock element, and, in response to the overburden stress (or vertical stress σ_ov = σ_V), there are the horizontal stresses, Minor Horizontal Stress (σ_h) and Major Horizontal Stress (σ_H), which avoid the lateral deformation induced by vertical loading. Estimating the ratio between in situ rock stresses, it is assumed that vertical and horizontal stresses are principal values as reported and discussed in [20]. In this research, it is assumed the same considerations with a fixed stress ratio. In-situ stresses are subjected to the static equilibrium equation:

$$\mathrm{div}\left(\mathit{\sigma }\right)+\mathit{b}=0,$$(1)

where, σ is the stress tensor and b is the vector of the body forces.

A new state of stress is created around the wellbore wall by the drilling process. This new stress state may generate significant change in formation integrity, leading to rock failure and operational problems [6, 21]. The formation around the borehole wall experiences deformation due to the relief of stresses, once the rock is removed, and then, replaced by the drilling fluid. The drilling fluid (or drilling mud) is a mixture of a fluid (water or oil based) and solid [22], used in the drilling operation to avoid the borehole closure and to fill the open space once occupied by the rock. The drilling fluid temporarily supports the borehole wall while drilling. The drilling fluid’s density must be high enough to prevent flow of formation fluids into the wellbore (a phenomenon named as “kick”) as well as wellbore closure, but low enough to prevent hydraulic fracture, formation failure and circulation loss [23]. Thus, the mud weight (or drilling fluid density) is responsible for keeping the wellbore walls stable. In order to minimize the cited operational risks, geomechanical models have been used in the design and drilling operations stages. Understanding the new stress state of vertical, inclined, and horizontal wellbore is essential to ensure a safe operation. The stress state simulation allows the definition of the optimal drilling fluid density to be used in a certain formation.

3 Methodology

A geomechanical model is introduced in order to simulate the formation behavior when mud pressure is applied to the borehole wall. The numerical model presented herein is developed using the Finite Element Method (FEM) and assumes linear elastic and isotropic material behavior under the plane-strain state. The method is implemented in NeoPZ, which is an open-source library for the development of finite element simulations [24]. A yield criterion is selected to describe the formation’s plastic region at the borehole wall. The geostatistics analysis is developed based on a stochastic field model. The methodology used to implement this model is presented in this section.

3.1 The geomechanical model

Once the wellbore length is very long (kilometers long sometimes) compared to the size of the studied domain in x and y directions, a plain strain hypothesis is relevant considering null displacement in the z-direction. The simulation is, therefore, two-dimensional, which significantly reduces computational time. As depicted in Figure 1, the finite element mesh is selected as a cross-section of the wellbore. The borehole cavity is the central circumference, and the remaining area refers to the original formation being drilled.

The FEM algorithm starts at finding the displacement u (x,y) from the equilibrium equation,

$$\{\begin{array}{c}\mathrm{div}\left(\mathit{\sigma }\right)+\mathit{b}=0\\ \mathit{u}={\mathit{u}}_{\mathbf{D}}\hspace{0.17em}\hspace{0.17em}\mathrm{on}\hspace{0.17em}\hspace{0.17em}\delta {\mathrm{\Omega }}_{\mathrm{D}}\\ \mathit{\sigma }\cdot \mathit{n}={\mathit{\sigma }}_{\mathbf{InSitu}}\cdot \mathit{n}\hspace{0.17em}\mathrm{on}\hspace{0.17em}\delta {\mathrm{\Omega }}_{\mathrm{N}}\end{array}$$

where,

• σ is the stress tensor;

• b is the body force vector;

• u is the displacement vector;

• u _D is a known function;

• σ _InSitu is the in-situ stress;

• Ω is the domain, and applying:

• $\begin{array}{ccc}\mathit{\sigma }\cdot \mathit{n}={\mathit{\sigma }}_{\mathbf{InSitu}}\cdot \mathit{n}& \mathrm{on}& \delta {\mathrm{\Omega }}_{\mathrm{N}}\end{array}$

  δΩ_D and δΩ_N are the boundary conditions, Dirichlet and Neumann, respectively.

The space of testing functions is defined by:

$$\begin{array}{ccc}V\left(\mathrm{\Omega }\right)=\{\mathit{v}\in H^1(\mathrm{\Omega }); \mathit{v}=0& \mathrm{on}& \delta {\mathrm{\Omega }}_{\mathrm{D}}\end{array} \},$$

where,

$$H^1\left(\mathrm{\Omega }\right)=\left\{\mathit{v}\in L^2\left(\mathrm{\Omega }\right);\nabla \mathit{v}\in L^2\left(\mathrm{\Omega }\right)\right\}.$$

By multiplying the equilibrium equation by the testing function v , integrating over the domain Ω, and using the divergence theorem, we have that:

$$-\underset{\mathrm{\Omega }}{\int \mathrm{}}\nabla \mathit{v}\cdot \mathit{\sigma } \mathrm{d}\mathrm{\Omega }+\underset{\delta \mathrm{\Omega }}{\int \mathrm{}}\mathit{v}\cdot \mathit{\sigma }\mathit{n} \mathrm{d}s+\underset{\mathrm{\Omega }}{\int \mathrm{}}\mathit{v}\cdot \mathit{b} \mathrm{d}\mathrm{\Omega }=0,$$

Assuming that between the elements,

$$\underset{\mathrm{\Gamma }}{\overset{}{\int }}⟦\mathit{v}\cdot \mathit{\sigma }\mathit{n}⟧\approx 0,$$

and applying,

$$\begin{array}{ccc}\mathit{\sigma }\cdot \mathit{n}=g& \mathrm{on}& \delta {\mathrm{\Omega }}_{\mathrm{N}}\\ \mathit{v}=0& \mathrm{on}& \delta {\mathrm{\Omega }}_{\mathrm{D}.}\end{array}$$

The forming function can be stated as:

$$\begin{array}{cc}a(\mathit{u},\mathit{v})=f\left(\mathit{v}\right)& \forall \mathit{v}\in V\left(\mathrm{\Omega }\right).\end{array}$$(2)

Given that:

$$a\left(\mathit{u},\mathit{v}\right)=\underset{\mathrm{\Omega }}{\sum \mathrm{}\int \mathrm{}}\nabla \mathit{v}\cdot \mathit{\sigma }\mathrm{d}\mathrm{\Omega },$$

$$\mathrm{f}\left(\mathit{v}\right)=\underset{\mathrm{\Omega }}{\sum \mathrm{}\int \mathrm{}}\mathit{v}\cdot \mathit{b} \mathrm{d}\mathrm{\Omega }+\underset{\delta {\mathrm{\Omega }}_N}{\int \mathrm{}}\mathit{v}\cdot g \mathrm{d}s,$$

and by the stress–strain relation for linear elastic condition, we find the stress tensor:

$$\mathit{\sigma }=\mathbf{C} \mathit{\epsilon }+{\mathit{\sigma }}_{\mathbf{0}},$$(3)

with C being the elastic constant matrix, σ ₀ the initial stress tensor and ε the strain tensor.

DiMaggio-Sandler yield criterion [25] is selected and implemented in the geomechanical model to describe the failure area/plastic region surrounding the wellbore. This failure criteria is chosen once it describes well the constitutive behavior of rocks [26–28]. The failure criteria yield function Φ is composed of two functions. In the present analysis, however, only the failure surface envelope function is considered. Therefore, the envelope which limits the elastic deformation under confinement is disregarded. Plastic deformation occurs when Φ > 0. In the present study, only the collapse failure is analyzed herein. For the sign convention, compressive stress is negative and tensile stress is positive.

It is important to mention that the present geomechanical model has been validated. An analytical approach is fully described in [29] to determine the stress state around wellbore for isotropic homogeneous rocks. This analytical approach, however, does not take into account rock elastic parameters such as Young’s Modulus and Poisson’s ratio, except by the axial stress that considers Poisson’s ratio in its equation. Thus, this work also aims to verify the influence of rock elastic modulus uncertainty as an input parameter. The geomechanical model solution presented herein was compared to this analytical approach; the comparison presented satisfactory agreement, as fully detailed in [30].

All analyses are based on a FE mesh of 750 elements in total, consisting of 25 radial elements and 30 circumferential elements; a quadratic approximation is used, and Dirichlet boundary condition is applied in the external part of the mesh to avoid rigid body displacement. The mud weight, or drilling fluid density pressure, is applied at the wellbore wall, and the formation’s initial stress state is defined by its in-situ stresses.

3.2 Spatially correlated field for heterogeneous formation

Heterogeneous rock mechanics properties may significantly influence the formation behavior, and, as a consequence, also affect the mud weight to be used in the operation. Probabilistic evaluation of a given field allows a more realistic analysis of the rock stress state, and assist, with better accuracy, selecting the best density for the drilling fluid. The pore pressure and fracture gradients are very narrow, and a small change can affect wellbore stability [5]. Therefore, this stochastic field modeling represents an important issue to be addressed in wellbore stability analysis.

Since the model is a plane model, the simulation is set for a specific layer or rock formation. In other words, it assumes one average mechanical property for the whole field. Therefore, the stochastic model is built with one layer at a time.

The finite element mesh is built using a geometric progression, which leads to a greater mesh refinement close to the wellbore wall, where the greatest change in stress state is found. The stochastic grid follows the same geometry of the finite element mesh. Thus, each element receives a random but correlated value of Young’s Modulus. Following the geometric progression, the variation is greater closer to the wellbore wall than over the remaining domain. Nevertheless, it is also possible to not use the geometric progression, and then the elements would be equally distributed, and the Young’s Modulus would also be uniformly defined in the mesh.

A spatially correlated field is used to evaluate the variability of the formation’s material properties. The correlation between each pair of elements is defined by a covariance function. In this model, the random variable is the formation’s elastic modulus (Young’s Modulus). The Squared Exponential covariance function is chosen for this model, which is often used to describe soil properties autocovariance [31]. The covariance function, also known as “kernel,” is as follows, from [32]:

$$k=\exp \left(-\frac{r^2}{{\gamma }^2}\right).$$(4)

The selected kernel allows the user to control the distance for which the correlation between the elastic modulus of two elements is not relevant. The control is set by the constant γ, which is associated with the spatial frequency in the field [33]. Therefore, the scale γ defines the smoothness of the stochastic field.

The covariation function also depends on the distance between the centroids of the elements, defined by r. A matrix R_m,m holds all distances between each pair of elements, calculated by their elements centroids coordinates (X):

$$R_{\mathrm{ij}}=\left(\begin{array}{llll}\sqrt{{\left(X_1-X_1\right)}^2}& \sqrt{{\left(X_1-X_2\right)}^2}& \cdots & \sqrt{{\left(X_1-X_{\mathrm{m}}\right)}^2}\\ \sqrt{{\left(X_2-X_1\right)}^2}& \sqrt{{\left(X_2-X_2\right)}^2}& \cdots & \sqrt{{\left(X_2-X_{\mathrm{m}}\right)}^2}\\ ⋮& ⋮& \ddots & ⋮\\ \sqrt{{\left(X_{\mathrm{m}}-X_1\right)}^2}& \sqrt{{\left(X_{\mathrm{m}}-X_2\right)}^2}& \cdots & \sqrt{{\left(X_{\mathrm{m}}-X_m\right)}^2}\end{array}\right).$$

The correlation matrix K_m,m, holds the correlation between all pair of elements, and it is a function of the matrix R and the scale γ:

$$K{(\gamma ,R)}_{\mathrm{ij}}=\left(\begin{array}{llll}{k}_{\mathrm{1,1}}& k_{\mathrm{1,2}}& \cdots & k_{1,m}\\ k_{\mathrm{2,1}}& k_{\mathrm{2,2}}& \cdots & k_{2,m}\\ \cdots & \cdots & \ddots & ⋮\\ k_{m,1}& k_{m,1}& \cdots & k_{m,m}\end{array}\right).$$

The matrix K is then decomposed by the singular value decomposition method, such that $K=U S V^{\mathrm{T}}$. The product of the left singular vector U and the square root of the diagonal matrix S, multiplies a Gaussian random and uncorrelated vector d _rand , generating a correlated vector d _corr .

The random but correlated vector d _corr is then scaled multiplying it by a standard deviation (SD_E) and summed to an average Young’s Modulus (E_avg).

$$E_{\text{stochastic}\hspace{0.17em}\text{field}}=d_{\mathbf{corr}}\ast S{D}_{\text{E}}+E_{\text{avg}},$$(5)

E_{stochastic field} refers to a stochastic but correlated field of Young’s Modulus. The variable SD_E is the standard deviation of the Young’s Modulus, and E_avg is the Young’s Modulus, both given as input. Figures 2 and 3 depict for the Young’s Modulus variable, the difference between a random correlated field, which will be considered in the analysis, and a simple random and not correlated field.

3.3 A 3D spatially correlated field for inclined wellbore

A new approach must be considered for inclined wellbores. Inclined wellbores are obtained from drilling in a given direction and inclination angle. Inclination angle is the angle between the wellbore axis tangent and the local gravitational vector. Therefore, a vertical wellbore has an inclination angle of 0°, and a horizontal wellbore has an inclination angle of 90°. Direction (or azimuth) angle consists of the angle between the horizontal wellbore projection and the true geographic north. This angle is commonly referred to as “wellbore azimuth.” It varies from 0° to 360° and is measured in the clockwise direction starting at the geographic north [34].

Fjaer et al. [29] present a transformation of the in-situ stresses to the local stress tensor of an inclined wellbore that can be obtained from a rotation between the z and y axes. This operation can be represented in matrix notation as follows:

$$\left[{\mathit{\sigma }}^o\right]=\left[\mathit{Q}\right]\left[{\mathit{\sigma }}_{\mathbf{InSitu}}\right][\mathit{Q}{]}^T,$$(6)

where,

• σ _InSitu is the in-situ stresses;

• α is the wellbore direction/azimuth;

• β is the wellbore inclination;

• Q is the rotation matrix, composed by the rotation between the z and y axis.

In which Q and σ _InSitu are composed, in matrix notation, by:

$$\mathit{Q}=\left(\begin{array}{lll}\cos \left[\alpha \right]\cos \left[\beta \right]& \sin \left[\alpha \right]\cos \left[\beta \right]& -\sin \left[\beta \right]\\ -\sin \left[\alpha \right]& \cos \left[\alpha \right]& 0\\ \cos \left[\alpha \right]\sin \left[\beta \right]& \sin \left[\alpha \right]\sin \left[\beta \right]& \cos \left[\beta \right]\end{array}\right),$$

$${\mathit{\sigma }}_{\mathbf{InSitu}}=\left(\begin{array}{lll}{\sigma }_{\mathrm{H}}& 0& 0\\ 0& {\sigma }_{\mathrm{h}}& 0\\ 0& 0& {\sigma }_{\mathrm{V}}\end{array}\right).$$

A two-dimensional correlated field is considered for vertical wellbore since the mesh is perpendicular to the well. Inclined wellbore’ mesh is also perpendicular to the wellbore orientation, but embedded in the layer that the geomechanical model is being analyzed. Figures 4 and 5 illustrate the difference between the two cases.

The geomechanical model also assumes a two-dimensional approach for inclined wellbores. The simplification is possible due to the local stress tensor transformation, previously described. For the stochastic analysis, however, the mesh is inclined compared to the horizontal formation layer, according to the wellbore inclination and direction. In order to evaluate the heterogeneous field of a drilled layer and determine the correlation between elements, a three-dimensional correlated field is developed.

The three-dimensional approach follows five steps: (1) definition of the height and radius of a cylindrical layer, (2) selection of a two-dimensional mesh as if it had been in the vertical model, (3) selection of the number of cross sections inside the cylindrical layer using the same two-dimensional mesh, (4) rotation of an extra two-dimensional mesh inside the cylindrical layer following the wellbore direction and inclination, and (5) definition of the coordinates of all elements centroids. All the steps can be seen in Figure 6. All meshes are discretized in the same way, as well as the inclined mesh.

The correlation matrix K_m,m is, therefore, calculated for all pairs of elements centroids. The stochastic field over the inclined mesh (in red) is correlated with inclined neighboring elements, which allows the use of the two-dimensional geomechanical model. The methodology described herein to create three-dimensional correlated layers, allows the user to generate drilling fields expeditiously.

4 Results and discussions

This section presents the results of a case study carried out using the numerical model and analysis method proposed herein. The wellbore in question is a carbonatic rock, and it is already known that the selected drilling fluid pressure (P_w) would induce a plastic zone at the wellbore wall. The wellbore parameters and rock features for the yield criterion are reported in [26]; the model proposed herein is validated against the deterministic result in [30].

The input values used in the simulations are summarized in Table 1, where r_w and r_e are the wellbore radius and mesh’s external radius, respectively; the same values are used in the vertical and inclined wellbore. The carbonatic material parameters for DiMaggio-Sandler criterion are A: 152.52; B: 0.0015489, and C: 146.29. The plastic zone area is calculated using the determinant of the Jacobian matrix from each integration point.

4.1 Plastic zone for vertical and inclined wellbores

The plastic zone (obtained based on Φ – yield criterion envelope) at the borehole wall are compared between vertical and inclined wellbores considering homogeneous media. The analyses input follow Table 1, and the inclination and azimuth angles analyzed are reported in Table 2.

Figures 7–10 depict the yield function for V1 and I1 wellbores according to Table 2 nomenclature. Note that for the evaluated mud pressure of P_w = 19.5 MPa, the plastic area for the inclined wellbore is smaller than for the vertical wellbore.

Figure 11 presents the yield functions and Figure 12 the plastic zones for wellbores with same direction (α = 30°), but different inclination, as reported in Table 2. The inclined meshes are plotted with the vertical mesh, and their yield functions and plastic zone are compared. It can be seen that the plastic area is gradually decreasing as the wellbore increase its inclination. In all cases, the same wellbore mud pressure is assumed.

The plastic zone of a heterogeneous media is compared to a homogeneous one for a vertical wellbore. Both models assume the same input parameters, except for Young’s Modulus. Note that a homogeneous field simulation – with a single value of Young’s Modulus for the entire domain – is a common practice by petroleum engineers. This average Young’s Modulus value is used as the mean value necessary to define the stochastic field. The standard deviation considered in this study is based on a Coefficient of Variation (CoV) of 10%. For a homogeneous and a heterogeneous field, Figure 13 depicts Young’s Modulus distribution in the field, the plastic zone at the borehole wall, and the plastic region areas. The total area of the plastic zone for the homogeneous field is 13.901 × 10⁻⁴ m² and 12.631 × 10⁻⁴ m² for this heterogeneous field example case. It can be noted that the plastic region around the wellbore for a homogeneous field is symmetric. For the heterogeneous case, however, the plastic zone surrounding the wellbore is not symmetric, presenting 56% of the total area on the top side and 44% of the total plastic zone area on the downside of the third column of Figure 13.

The scale factor used in this work for all heterogeneous case analysis is γ = 0.5, which is the ratio of the first element at the wellbore wall. The asymmetry commented above may lead to operational issues during a drilling operation. Once mud properties are designed to stabilize the most vulnerable side of the well, a drilling fluid density considered safe if drilling a homogeneous formation may present a higher and asymmetric plastic zone when a heterogeneous analysis is performed.

4.2 Stochastic analysis: inclined wellbores

The analyzed inclined wellbores present a decrease in the plastic zone area as the inclination increases when considering homogeneous media. Thus, an analysis is made in order to evaluate the heterogeneity effects in wellbores’ inclination and direction. The analysis assumes the same input parameters as Table 1, but for this inclined study, the three-dimensional stochastic field generation method is used. Therefore, a cylinder is built for the stochastic analysis. For the present case, the cylinder has eight equal cross-sections, and each section includes a two-dimensional and horizontal mesh, while the inclined mesh is embedded within the cylinder.

In order to reduce computational processing time, the mesh is changed to 2 m of the external radius. The number of elements and the wellbore radius are kept the same as in all analyses. The space between each mesh in the cylinder is 0.5 m, and the cylinder’s total height is 4 m. Ten thousand cases for each wellbore inclination is presented in Figure 14. The vertical wellbore considered in the analyses have the same direction as considered in the inclined wellbores, but no inclination angle (wellbore named V2). The direction was assumed in the stochastic analysis, in order to evaluate the influence of the direction in the yield function. One can note that as well as in the homogeneous domain, the heterogeneous field cases also reduce its total plastic zone area (A_T) as inclination increases. Figure 16 compares the total plastic zone area histogram for the same vertical well with and without direction angle. The histograms depict that the azimuth has little effect in the vertical distribution dispersion. The Complementary Cumulative Distribution Function (CCDF) of each histogram are depict in Figures 15 and 17, which allows the user to verify the failure probability of A_T > x, 0 < x < 20 × 10⁻⁴ m².

Table 3 presents the total plastic area mean (${\mu }_{A_{\mathrm{T}}}$) of each stochastic case, as well as their CoV and the deterministic total area for the homogeneous field ($A_{\mathrm{T}}^{\mathrm{H}}$). The highest CoV is found in the most inclined wellbores. One, however, may associate the dispersion of the variable under analysis to the three-dimensional approach for stochastic inclined fields; thus, the same study presented in Figure 14 is evaluated disregarding the 3D field model. The result gives no significant change in ${\mu }_{A_{\mathrm{T}}}$ and CoV values, with differences lower than 2%. Therefore, the small influence on the plastic zone area under the analyzed mud pressure is likely due to the very small plastic zone area under analysis. Also, one of the covariance function main feature is the distance between elements, and at the borehole wall, the length of the elements is significantly small due to the geometric progression. Thus, around the borehole wall, Young’s Modulus variability is not as high as it is in the farfield, even if the mesh has an inclination, once the plastic zone area is still very close to the borehole wall with very small elements.

4.3 Probabilistic analysis: maximum plastic deformation area

As reported by Muller et al. [15], probabilistic analyses allow the user to verify the failure probability of wellbores based on a maximum value of the plastic zone area, according to the simulated drilling fluid’s density. If the maximum allowed area (A_y) of the plastic region at the borehole wall is a percentage of the wellbore cavity area (A_w), the failure probability of reaching the maximum failure area could be verified by:

$$P_{\mathrm{f}}=P\left(A_{\mathrm{T}}>A_{\mathrm{y}}\right),$$(7)

where, A_y = p% of A_w, and A_T is the total area of plastic region.

In order to analyze the uncertainty of the drilling fluid pressure, a probabilistic analysis is carried out for a pressure range from 10 MPa to 35 MPa for the vertical wellbore studied in this work, and considering the same failure criteria. Ten thousand simulations are carried out for each mud weight. Additionally, random normal distributions are assigned for the in-situ stresses. The distribution mean values are the respective stresses from Table 1, with a CoV of 10%. The random distributions follow the criteria:

$${\sigma }_{\mathrm{H}}>{\sigma }_{\mathrm{V}}>{\sigma }_{\mathrm{h}}.$$

Figure 18 depicts the plastic area for each simulation of ten thousand simulations for each mud pressure, while Figure 19 represents the failure probability for a given mud pressure, failure herein is given in a performance base manner based on account of plastic area being higher than 10%, 5%, and 3% of A_w. Considering 3% of A_w as an expected performance for the considering wellbore, no collapse would occur if using mud drilling’s pressure over 31 MPa since failure probability is close to zero. On the other hand, if 10% of A_w is the expected performance of the wellbore, the drilling mud weight shall be greater than 24 MPa to avoid shear failure. Note that the mud weight’s upper limit is defined by rock tensile strength, which is not evaluated in this work.

5 Conclusion

A stochastic geomechanical model is presented to predict formation behavior in drilling operations considering heterogeneous fields. The simulation method proposed herein predicts the total area of the plastic region surrounding the wellbore once it is submitted to a mud/drilling fluid’s pressure. The results show that, when considering a heterogeneous field, a non-symmetric plastic region may occur. This lack of symmetry may lead to operational difficulties.

A new procedure is proposed to evaluate stochastic fields of the inclined wellbore; this procedure allows the modeling of thick formation layers based on two-dimensional geomechanical models with a low computational cost if compared to full 3D random field analyses.

A wellbore case is evaluated, and its results considering homogeneous and heterogeneous field are compared. The studied wellbore, when subject to a specific mud weight, exhibits different total plastic zone area according to its inclination. It is observed that the entire plastic area becomes smaller as the wellbore inclination increases. Therefore, considering the wellbore features analyzed herein, the mud pressure evaluated is safer for a wellbore with greater inclination, because the total plastic zone area is more significant for the vertical well compared to the inclined wellbore. This behavior can be seen in the homogeneous and heterogeneous field. Wellbore azimuth depicted a small influence in the distribution dispersion for vertical wells.

The total plastic zone area distributions determined in the stochastic analysis lead to an evaluation of the probabilistic failure. Moreover, it leads to information regarding the safety of a drilling operation, once an acceptable plastic region area at the borehole wall is defined by the group responsible for the operation. Note that a mud pressure considered safe in a homogeneous field analysis, may not be indicated based on the probabilistic study framework detailed herein.

The failure probability for a mud weight range can help in decision-making during drilling operations in heterogeneous formations. Besides, one is able to define a target wellbore performance based on plastic deformation area and verify the failure probability when using a specific drilling fluid’s density.

Acknowledgments

The authors are indebted to National Agency of Petroleum, Natural Gas and Biofuels (ANP) and to Brazilian Petroleum Corporation – Petrobras, Grant: 2014/00090-2, for funding this research project. Luiz C.M. Vieira Jr. is indebted to national council for scientific and technological development (CNPq) for the research productivity fellowship grant 304005/2017-7. Philippe R. B. Devloo is indebted to CNPq for the research productivity fellowship grant 305823/2017-5, and to São Paulo Research Foundation (FAPESP) for the research grant 2017/15736-3. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors only and do not necessarily reflect the views of the sponsors or affiliates. The authors would like to thank Dr. Thiago Dias dos Santos and Mr. Patrick Kherlakian for helpful discussions and contributions.

References

All Tables

All Figures