Thermoporoelastoplastic Wellbore Breakout Modeling by Finite Element Method

Drilling a hole into rock results in stress concentration and redistribution close to the hole. When induced stresses exceed the rock strength, wellbore breakouts will happen. Research on wellbore breakout is the fundamental of wellbore stability. A wellbore breakout is a sequence of stress concentrations, rock falling, and stress redistributions, which involve initiation, propagation, and stabilization sequences. Therefore, simulating the process of a breakout is very challenging. Thermoporoelastoplastic models for wellbore breakout analysis are rare due to the high complexity of the problem. In this paper, a fully coupled thermoporoelastoplastic finite element model is built to study the mechanism of wellbore breakouts. The process of wellbore breakouts, the influence of temperature and the comparison between thermoporoelastic and thermoporoelastoplastic models are studied in the paper. For the finite element modeling, the D-P criterion is used to determine whether rock starts to yield or not, and the maximum tensile strain criterion is used to determine whether breakouts have happened.


Introduction
When a hole is drilled into the crust, the rock mass removed from the hole will not support the surrounding rocks, which leads to the stress concentration and redistribution. If the induced stresses exceed the rock strength, parts of rock fall from the wellbore wall, which is called a wellbore breakout [1,2]. Wellbore breakouts result in problematic wellbore stability and affects drilling efficiency, so the research of wellbore breakouts is important in petroleum engineering [3,4].
Numerous studies have been made to explain the mechanism of wellbore breakouts. Zoback found the analytical solution of initial breakout zone by incorporating the Mohr-Coulomb criterion to the Kirsch equation [2]. The studies of wellbore breakouts show that the wellbore breakouts occur by a series of successive spalls that result from shear failure subparallel to the direction of the local minimum principal stress [3,[5][6][7][8][9]. In addition, some studies also show that a wellbore breakout can also influence the wellbore breakdown pressure [10]. To consider the effect of pore pressure and obtain effective stress distribution in a porous medium, Biot [11] proposed a linear poroelastic model based on Terzaghi's theory and introduced Biot's coefficient. Later, Fluid seepage and pressure diffusion between the wellbore and the formation were considered for both vertical and inclined wells [12][13][14]. Recently, numerical methods were also used to analyze the breakout mechanisms [15][16][17][18][19]. Most of these models do not consider the influence of plastic damage of the rock. For the plastic model, different yield criteria are used to analyze the character of rock and soil. The Mohr-Coulomb criterion is the most popular one, which relates the shearing resistance to the contact forces and friction [20]. The Drucker-Prager criterion is an extended version of the Von Mises criterion [21]. Whittle and Kavvadas [22] presented the MIT-E3 soil model to describing the behavior of overconsolidated clays that obey normalized behavior and are rate-independent. Akl and Whittle [23] analyzed horizontal wellbore stability in clay shale based on MIT-E3 soil model. Zhang and Yin [24] made a wellbore study using the tangent stiffness matrix method based on the Drucker-Prager criterion.
However, in reality, a wellbore breakout is a complex time-dependent process including initiation, propagation, and stabilization sequences [25]. The breakout initiation and its propagation are the result of stress concentration on the wellbore wall, further stress concentration at the tip of breakout, and the formation of a plastic zone around the tip [26]. Poroelastoplastic analysis of wellbore breakouts is still not well studied due to the complexity of the problem.
With the development of a deep well, consideration of the thermal effects is essential. Coupled thermal-hydraulic-mechanical processes play an important role in the stability of wells in thermal reservoir formations [27]. Lewis [28] used FE simulation to study thermal recovery processes and heat losses problems to surrounding strata. Aboustit et al. [29] used a general variational principle to investigate thermo-elastic consolidation problems, and Vaziri [30] also presented a fully coupled thermo-hydro-mechanical FE model.
In this paper, a fully coupled thermoporoelastoplastic finite element model is built to study the mechanism of wellbore breakouts. The Drucker-Prager criterion is used to determine whether rock starts to yield or not, and the maximum tensile strain criterion is used to determine whether breakouts happened.

Model Structure and Methodology
To simulate wellbore breakouts in plastic rock, the Drucker-Prager yield criterion and maximum tensile strain criterion are used to analyze the yield and failure of the rock. The finite element implementation of thermoporoelasticity is based on Biot's theory with the compressible fluid flowing through the saturated non-isothermal deformable porous medium.

Rock Failure Criteria
In this paper, the Drucker-Prager yield criterion and maximum tensile strain criterion are used to analyze the yield and failure of the rock. Figure 1 shows a simplified sketch for the relationship between stress and strain in the condition of yield and failure of the rock, where σ yield and ε yield are the yield stress and strain of rock, ε breakout is the maximum allowable tensile strain of rock. The Drucker-Prager yield criterion is an extended version of the Von Mises criterion. To obtain a smooth yield surface approximate to Mohr-Coulomb surface, Drucker and Prager [21] put forward the following yield criterion:  The Drucker-Prager yield criterion is an extended version of the Von Mises criterion. To obtain a smooth yield surface approximate to Mohr-Coulomb surface, Drucker and Prager [21] put forward the following yield criterion: where is F is the yield criterion, I 1 is the first stress invariant, J 2 is the second deviatoric stress invariant, and α, k 0 are the material constants. Two of the most common approximations used are obtained by making the yield surfaces of the Drucker-Prager and Mohr-Coulomb criteria coincident either at the outer or inner edges of the Mohr-Coulomb surface. Coincidence at the outer edges is obtained when whereas, coincidence at the inner edges is given by the choice where c is the cohesion of the material; ϕ is the internal friction angle of the material. The yield condition of rock can be determined by the following equation.
The failure (breakout) condition of rock can be determined by the following equation.

Finite Element Implementation
The general theory of thermoporoelasticity is based on Biot's theory. With the compressible fluid flowing through the saturated non-isothermal deformable porous medium, in the form of displacements, pressure and temperature as unknowns, the governing equations can be described as follows [31].
where K, K f and K m are the bulk moduli of the skeleton, fluid and matrix, respectively. k is the permeability of the porous medium, µ is the viscosity of the fluid, G and λ are Lamé constants. u, p and p t denote the displacement of the porous medium, the pore pressure and its time derivative, φ is the porosity of the porous medium, λ T is the thermal conductivity matrix of the porous media, T and T t are the temperature and time derivative, ρ s c s is the heat capacity of the solid phase, ρ f c f is the heat capacity of the fluid phase, β s is the thermal expansion coefficient of the matrix, and β f is the thermal expansion coefficient of the fluid. Furthermore, i T = [1, 1, 1, 0, 0, 0], and D is the elastic stiffness matrix, the subscript t denotes time derivative.
Using Galerkin finite element method to approximate the governing equations, the final form of the FEM solution is as follows.
T t ] T are the vectors of unknown variables and corresponding time derivatives, and f u , f p , f T T is the vector for the nodal loads, flow source and heat sources. The explicit expressions of above matrixes are as follows.
To integrate the above equations with respect to time (θ method), the equation becomes: where ∆u ∆p is the increment of displacement, pore pressure, ∆ f u ∆ f p is the increment vector for the nodal loads and flow source, ∆t is the increment of time, θ is the time integration coefficient, f p 0 is the initial vector for the nodal loads in a time step, and p 0 is the initial vector for the pore pressure in a time step.
If the stresses are in the plastic state at a special time t a based on the yield criterion, Equation (19) becomes Equation (20) at the special time t a [32].
where Σ e is the set of element nodal forces round the node, and Ψ n is the unbalanced force.
For an ideal elastoplastic body, A = 0.
In the process of iteration, the convergence condition can be described as: where norm (Ψ n ) is the norm of vector Ψ n , and V is the convergence precision. For different iterative methods, the global stiffness matrixes are different, which are shown in Figure 2. For the constant stiffness method, the values in Equation (20) can be calculated based on Equations (10)- (18). For the tangent stiffness method, the explicit expressions of the values in Equation (20) are as follows.
Mining 2022, 2, FOR PEER REVIEW 6 where is elastic-plastic matrix [32,33]. For an ideal elastoplastic body, = 0. The tangent stiffness method is adopted in this paper.

Numerical Experiments
Thermoporoelastoplastic wellbore breakouts are simulated for a vertical well, and there are four parts in this section: the relationship between breakout shape and in situ stresses, verification of wellbore breakout process, influence of drilling fluid temperature, The tangent stiffness method is adopted in this paper.

Numerical Experiments
Thermoporoelastoplastic wellbore breakouts are simulated for a vertical well, and there are four parts in this section: the relationship between breakout shape and in situ stresses, verification of wellbore breakout process, influence of drilling fluid temperature, and a comparison between thermoporoelastic and thermoporoelastoplastic modeling.

Thermoporoelastoplastic Finite Element Modeling of Wellbore Breakouts
For a vertical wellbore shown in Figure 3 that is subjected to horizontal in-situ stresses σ H and σ h , the shape of breakouts φ b and r b can be acquired by finite element modeling, where r b is the depth of breakouts, and φ b is the width of breakouts (Figure 4). The tangent stiffness method is adopted in this paper.

Numerical Experiments
Thermoporoelastoplastic wellbore breakouts are simulated for a vertical well, and there are four parts in this section: the relationship between breakout shape and in situ stresses, verification of wellbore breakout process, influence of drilling fluid temperature, and a comparison between thermoporoelastic and thermoporoelastoplastic modeling.

Thermoporoelastoplastic Finite Element Modeling of Wellbore Breakouts
For a vertical wellbore shown in Figure 3 that is subjected to horizontal in-situ stresses and ℎ , the shape of breakouts and can be acquired by finite element modeling, where is the depth of breakouts, and is the width of breakouts ( Figure  4).  The yield condition of rock is determined by Equation (4), and the failure (breakout) condition of rock is determined by Equation (5).
The geometric and mechanical parameters are shown in Table 1.  The yield condition of rock is determined by Equation (4), and the failure (breakout) condition of rock is determined by Equation (5).
The geometric and mechanical parameters are shown in Table 1.  The results of the breakout shape are shown in Table 2. From Table 2, Figures 5 and 6, it can be seen that the relationship between breakout shape and in situ stresses is nonlinear. With relatively large σ H and difference between σ H and σ h , as σ h increases, the breakout width and depth increase. With relatively small σ H and difference between σ H and σ h , as σ h increases, the difference between σ H and σ h is smaller, the in situ stress is closer to symmetric in situ stresses state, and the breakout width and depth decrease. Table 2. Result of the shape of wellbore breakouts r b and φ b by FEM.

Verification of Wellbore Breakout Process
For a vertical well shown in Figure 5, associated with horizontal in situ stresses σ H = 87.5 MPa, σ h = 72.5 MPa and other parameters are listed in Table 1, the process of well breakout is shown in Table 3 and Figures 7-9. Table 3. Data of effective principal stresses and breakout region in the process of breakouts.

Iteration1
Iteration2 Iteration3         From Table 3 and Figures 7-9, the principal stresses of elements close to the tip of wellbore breakouts are increasing in the process of breakouts, but the maximum tensile strain changes to less than . As a breakout develops, the breakouts depth increases and the breakout width decreases, which means the depth of breakouts increases till a stable state, but the width of breakouts remains unchanged [3].

Influence of Drilling Fluid Temperature
For a vertical well shown in Figure 3, associated with horizontal in situ stresses = 85 MPa, ℎ = 72.5 MPa, and the temperature difference between drilling fluid and surrounding rock ∆ = −50~50K, and other parameters are listed in Table 1, the breakout shape for different drilling fluid temperature is shown in Figure 10. From Table 3 and Figures 7-9, the principal stresses of elements close to the tip of wellbore breakouts are increasing in the process of breakouts, but the maximum tensile strain changes to less than ε breakout . As a breakout develops, the breakouts depth increases and the breakout width decreases, which means the depth of breakouts increases till a stable state, but the width of breakouts remains unchanged [3].

Influence of Drilling Fluid Temperature
For a vertical well shown in Figure 3, associated with horizontal in situ stresses σ H = 85 MPa, σ h = 72.5 MPa, and the temperature difference between drilling fluid and surrounding rock ∆T = −50 ∼ 50K, and other parameters are listed in Table 1, the breakout shape for different drilling fluid temperature is shown in Figure 10. Based on Figure 10, it can be seen that as temperature of drilling fluid increases, the breakout width and depth will increase.

Comparison between Thermoporoelastic and Thermoporoelastoplastic Modeling
For a vertical well, shown in Figure 3, which is subjected to horizontal in-situ stresses Based on Figure 10, it can be seen that as temperature of drilling fluid increases, the breakout width and depth will increase.

Comparison between Thermoporoelastic and Thermoporoelastoplastic Modeling
For a vertical well, shown in Figure 3, which is subjected to horizontal in-situ stresses σ H and σ h , thermoporoelastic and thermoporoelastoplastic modeling are compared in this section, and parameters are listed in Table 1.
For the thermoporoelastoplastic model, the yield condition of rock is determined by Equation (4), and the failure (breakout) condition of rock is determined by Equation (5). For the thermoporoelastic model, the failure (breakout) condition of rock can be determined by Equation (4). Figures 11 and 12 show different breakout depths and widths corresponding to different in situ stresses by thermoporoelastic and thermoporoelastoplastic modeling. Based on Figure 10, it can be seen that as temperature of drilling fluid increases, the breakout width and depth will increase.

Comparison between Thermoporoelastic and Thermoporoelastoplastic Modeling
For a vertical well, shown in Figure 3, which is subjected to horizontal in-situ stresses and ℎ , thermoporoelastic and thermoporoelastoplastic modeling are compared in this section, and parameters are listed in Table 1.
For the thermoporoelastoplastic model, the yield condition of rock is determined by Equation (4), and the failure (breakout) condition of rock is determined by Equation (5). For the thermoporoelastic model, the failure (breakout) condition of rock can be determined by Equation (4). Figures 11 and 12 show different breakout depths and widths corresponding to different in situ stresses by thermoporoelastic and thermoporoelastoplastic modeling.  From Figures 11 and 12, it can be found that for different I -situ stresses, the width and depth of wellbore breakouts for thermoporoelastoplastic model are smaller than those for thermoporoelastic model. For the plastic model, the yield condition of rock is determined by Equation (4), and the failure (breakout) condition of rock is determined by From Figures 11 and 12, it can be found that for different I -situ stresses, the width and depth of wellbore breakouts for thermoporoelastoplastic model are smaller than those for thermoporoelastic model. For the plastic model, the yield condition of rock is determined by Equation (4), and the failure (breakout) condition of rock is determined by Equation (5). For the elastic model, the failure (breakout) condition of rock can be determined by Equation (4). Therefore, in elastic model, when F > 0 in Equation (4), a breakout happens, but in plastic model, rocks just enter plastic state and no breakouts happened. Therefore, the breakout widths and depths are greater for the elastic model compared to the plastic model.

Discussion
In this paper, wellbore breakout process is simulated by the finite element method. Plasticity and temperature influence are considered, and results are compared with elastic model. There are some limits for this study. First, mud cake is assumed to be perfect, and there is no liquid transfer across the mud cake. Second, ideal elastictoplastic stress-strain relationship is used. Third, results are not verified by experimental data and field data. Future work is expected to address this in the next step.

Conclusions
In this paper, the finite element method is employed to simulate wellbore breakouts based on the thermoporoelastoplastic model. Numerical experiments on finite element modeling of wellbore breakouts show the contrasting tendency of breakouts shape with different in situ stress, and some conclusions can be made as follows: 1.
The relationship between breakout shape and in situ stresses is nonlinear. If the σ h is constant, as σ H increases, the breakout width and depth become greater. If the σ H is kept constant while the difference between σ h and σ H is relatively large, as σ h increases, the breakout width and depth increase greater. However, if the σ H is kept constant while the difference between σ h and σ H is relatively small, as σ h increases, the breakout width and depth decrease instead.

2.
In the process of wellbore breakouts, the breakout depth increases till a stable state, but the breakout width remains unchanged.

3.
As temperature of drilling fluid increases, the breakout width and depth will increase.

4.
For different in situ stresses, the width and depth of wellbore breakouts for thermoporoelastoplastic model are smaller than those for thermoporoelastic model.