Next Article in Journal
Component-in-the-Loop Testing of Automotive Powertrains Featuring All-Wheel-Drive
Previous Article in Journal
Design and Implement of Three-Phase Permanent-Magnet Synchronous Wave Generator using Taguchi Approach
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Investigation of the Influence of the Drill String Vibration Cyclic Loads on the Development of the Wellbore Natural Fracture

1
College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China
2
State Key Laboratory of Petroleum Resources, China University of Petroleum, Beijing 102249, China
3
State Key Laboratory of Heavy Oil Processing, China University of Petroleum, Beijing 102249, China
4
School of Materials and Metallurgy, Wuhan University of Science and Technology, Wuhan 430081, China
5
School of Engineering Science, University of Science and Technology of China, Hefei 230027, China
6
College of Petroleum Engineering, Liaoning Petrochemical University, Fushun 113001, China
*
Author to whom correspondence should be addressed.
Energies 2021, 14(7), 2015; https://doi.org/10.3390/en14072015
Submission received: 4 January 2021 / Revised: 17 March 2021 / Accepted: 19 March 2021 / Published: 6 April 2021

Abstract

:
Wellbore instability is one of the most serious issues faced in the drilling process. During drilling operations, the cyclic loads applied on the fractured formation progressively modify the initial parameters (i.e., length and width) of the fractured formation, thus resulting into undesirable wellbore instability. In this paper, using a nonlinear finite element software (ABAQUS) as the numerical simulator, a poro-elasto-plastic model has been established which aimed at analyzing the influence of drill string vibration cyclic loads on the development of the wellbore natural fracture. The numerical results showed that the fracture width as a function of time profiles followed a sinusoidal behavior similar to the drill string vibration cyclic load profiles. For different cyclic load magnitudes with constant number of cyclic loads, the highest percentage increase of the fracture width after integration of cyclic loads was 64.77%. Interestingly, the fracture width increased with the fracture length in the near wellbore region while it globally decreased in the region far away from the wellbore. But for constant cyclic load magnitude with different number of cyclic loads, the biggest percentage increase of the fracture width after integration of cyclic loads was slightly lower representing 63.12% while the oscillating period of the fracture width increased with the number of cyclic loads. The parametric study revealed that the drill string vibration cyclic loads were relatively independent of the fracture length and the bottom hole pressure. However, the fracture width and the loss circulation rates were considerably impacted by the drill string vibration and the highest percentage increase of the loss circulation rate after integration of cyclic loads was 14.3%. This study provides an insight into the coupling of the fracture rock development and the continuous cyclic loads generated by drill string vibrations which is an aspect that has been rarely discussed in the literature.

1. Introduction

The conventional methods used to assess the wellbore stability analysis generally assume that the loads acting on the wellbore are in steady-state [1,2,3]. Nonetheless, recent researches have revealed that the magnitudes of the loads affecting the wellbore stability are more likely to change over time; therefore, the wellbore pressure is prone to fluctuations [4,5]. The drill string vibration cyclic load is one of the most common sources of wellbore pressure fluctuations and numerous researches have demonstrated its impact on wellbore stability: ref. [6] were among the first researchers to investigate the rock failure under dynamic conditions. They found the critical rotary speed at which the fluctuating axial loads applied on the drill string will induce excessive wellbore enlargement. Ref. [7] concluded in their research that when drilling in hard rocks such as basalt, the drill string acceleration generates lateral load that will impact the wellbore wall and later results into wellbore instability. Ref. [8] after conducting fatigue tests experiments concluded that: (i) if the stress exerted by the cyclic loads on the rock formation was above 87% of the rock formation collapse strength, the rock collapsed within 100 cyclic loads, (ii) when the stress generated by the cyclic loads acting on the rock formation ranged between 75% to 87% of the rock formation collapse strength, the rock formation was weakened but the rock itself remained uncollapsed, and (iii) if the stress applied by the cyclic loads on the rock formation was below 75% of the rock formation collapse strength, the rock formation neither weakened nor collapsed. Ref. [9] reported that thermal induced stresses caused by daily fluctuations in temperature can lead to crack propagation in the rock mass.
More so, ref. [10] studied the impact of drill string vibration on wellbore stability. They made a correlation between the drill string critical speed and the rock fragmentation volume. They also established a rock fragmentation law related to the drill string impact on the wellbore. Ref. [11] investigated the effect of formations compressive strength on drill string vibrations and demonstrated that uniaxial compressive strength of harder formation increases with the magnitude of vibrations and reversely. Ref. [12] investigated the effects of the drill string vibration cyclic loads on the wellbore stability and they obtained slightly similar results with those of [8]. They found that wellbore failure occurred due to drill string vibration cyclic loads when: (i) the drill string vibration applied stresses on the formation exceeded the strength of the rock formation, (ii) the drill string vibration applied continuous cyclic loads on the wellbore resulted in the rock fatigue, (iii) the drill string generates cyclic loads which did not result in the rock fatigue, however, the strength of the rock reduced. Ref. [13] investigated the effect of vibrations in drilling systems and concluded that revolutions per minute (RPM) and weight on bit (WOB) are the major parameters that affect drill string vibrations. Ref. [14] studied the mechanical behavior of shale rock under uniaxial cyclic loading and unloading condition and found that the residual strains was related to the number of cyclic loads to failure by an exponential function. In the same year, ref. [15] presented a general review of cyclic and fatigue behavior of rock materials from literature published over the preceding 50 years. They reported that cyclic loads generated a decohesion of rock grains and matrix loosening of the rock. A wider fracture zone was also observed on the samples after cyclic loadings.
Although some studies in the past have already investigated the effect of the drill string vibration on the wellbore stability, some of those studies [6,7,10] only demonstrated that drill string vibration indeed had affected the wellbore stability. Other researchers [8,12] focused on the effect of the cyclic loads on the wellbore collapse pressure. Few studies investigated how the continuous cyclic loads can change the parameters of a naturally fractured formation such as its length and its width. Meanwhile, it is well known that one of the primary functions of the drilling mud is to maintain the stability of the wellbore by outward hydrostatic pressure and mud filter cake [16]. In a naturally fractured formation, plugging and sealing materials are generally added to the drilling mud in order to seal the fractures and to ensure the wellbore stability [17,18,19]. To successfully seal the fractures and prevent loss circulation, it is vital to accurately predict the growth of the natural fractures especially when the wellbore pressure is subjected to fluctuations caused by some factors such as the drill string vibrations cyclic loads.
Due to the aforementioned reasons, a plane strain poro-elasto-plastic finite element analysis is carried out in this paper to investigate the influence of the drill string vibrations cyclic loads on the development of the wellbore natural fracture rock. The ABAQUS software is utilized as the numerical simulator and the results obtained were compared with small asymptotic analytical solutions to validate the numerical model.

2. Statement of Theory and Definitions

The investigation of the influence of the drill string vibration cyclic loads on the development of the wellbore natural fracture is a complex problem which involves the following phenomena:
(1)
Flow of the drilling fluid inside the wellbore;
(2)
Fracture propagation and fracturing fluid flow;
(3)
Porous medium deformation and pore fluid flow.

2.1. Flow of the Drilling Fluid inside the Wellbore

The flow of the drilling fluid inside the wellbore follows a U tube configuration as shown in Figure 1.
The fluid flows inside the drill pipe, reaches the bottom hole and finally flows back to the surface passing through the wellbore annulus. The pipe element model in Abaqus is used to simulate the flow of the fluid inside the wellbore. If P 1 and P 2 represent two nodes of the pipe element, then the governing equation for the pipe element model is given by [20] as:
Δ P = Δ P h y d r o   +   Δ P l o s s
Δ P h y d r o = ρ g Δ Z
Δ P l o s s = c l ρ V 2 2
where, Δ P is the pressure difference between P 1 and P 2 , Pa; Δ P h y d r o is the hydrostatic drilling fluid pressure between P 1 and P 2 , Pa; Δ P l o s s is the viscous pressure loss, Pa; ρ is the drilling fluid density, kg / m 3 ; g is the gravitational acceleration, m / s 2 ; Δ Z is the elevation difference between P 1 and P 2 , m; V is the fluid velocity in the pipe, m   / s ; c l is the fluid loss coefficient.

2.2. Fracture Propagation and Fracturing Fluid Flow

The Abaqus coupled pressure/deformation Cohesive Zone Model is used to model the fracture propagation and the fracturing fluid flow. The Cohesive Zone Model is defined with a traction separation law and a fracturing fluid flow law.

2.2.1. Traction Separation Law

The traction separation model is subdivided into three phases:
(1)
Initial loading before damage. In this phase, the material is supposed to have a linear elastic behavior. The linear elastic behavior is described with an elastic constitutive matrix which relates the nominal stresses to the nominal strains. The nominal stresses are the force components divided by the original area at each integration point, while the nominal strains are the separations divided by the initial thickness. The elastic constitutive matrix is given as:
t = { t n t s 1 t s 2 } = [ E n n E n s E n t E n s E s s E s t E n t E s t E t t ] { ε n ε s 1 ε s 2 } = E ε
where, t n , t s 1 and t s 2 respectively represent the normal nominal stress, the first shear and second shear nominal stresses, Pa; E is the tensile stiffness matrix of the cohesive element, ε n , ε s 1 and ε s 2 are respectively the normal nominal strain, the first shear and second shear nominal strains.
(2)
Damage initiation phase: The damage initiation refers to the beginning of degradation of the material point. The damage is assumed to be initiated when the stresses and/ or strains fulfill certain criteria. In this paper, the damage is assumed to be initiated when the following quadratic stress criterion is satisfied [21]:
( t n t n 0 ) 2   +   ( t s 1 t s 1 0 ) 2   +   ( t s 2 t s 2 0 ) 2 = 1
where, t n 0 , t s 1 0 and t s 2 0 are respectively the peak values of the normal nominal stress ( t n ) , the first shear ( t s 1 ) and the second shear ( t s 2 ) nominal stresses, Pa.
(3)
Damage evolution phase: The damage evolution phase corresponds to the progressive degradation of the interface stiffness after damage initiation. Several models can be used to describe damage evolution. In this research, a damage evolution law based on energy combined with the Benzeggagh-Kenane (BK) mixed mode behavior is defined to model the material stiffness degradation after damage initiation. By assuming equal critical fracture energies in the first and second shear directions ( G s C = G t C ) the governing equation of the BK mixed mode behavior model is given as:
G n C   +   ( G s C         G n C ) { G s t G T } η = G C
where, G n C , G s C and G t C respectively refer to the critical fracture energies required to cause failure in the normal, first and second shear directions,   J / m 2 ; G C is the total energy dissipated due to failure, J / m 2 ; G s t = G s   +   G t , G T = G s t   +   G n where G n , G s and G t are respectively the work done by the tractions in the normal, first and second shear directions, J; and η is the BK power law parameter.

2.2.2. Fracturing Fluid Flow

The Reynold’s lubrication theory governs the flow of the fluid within the fracture and is given by [22]:
w t   +   q f s   +   v t   +   v b = 0
where w is the fracture aperture, m; q f is the longitudinal flow rate per unit of width, m 2 / s ; v t and v b are respectively the normal velocities at the top and the bottom of the fracture, m / s . The fracturing fluid flow and the normal fluid velocities are given by [22]:
q f =     w 3 12 μ f p
v T = c T ( p f     p T )
Δ P l o s s = c l ρ V 2 2
where, μ f is the viscosity of the fracturing fluid, Pa. s; p   is the gradient of the fracturing fluid pressure along the fracture surface, Pa/m; c T and c B are respectively the leak off coefficients at the top and the bottom of the fracture, m/s/Pa; p f   is the fracturing fluid pressure along the surface of the fracture, Pa; p T and p B are respectively the pore fluid pressures at the top and the bottom of the fracture, Pa.

3. Porous Medium Deformation and Pore Fluid Flow

A coupled pore-pressure/deformation continuum finite elements model is used to model the pore fluid flow inside the rock formation and the porous medium deformation.

3.1. The Pore Fluid Flow

The continuity equation combined with the Darcy’s law is used to describe the flow of the fluid in the porous medium. Assuming small volumetric strains, the pore fluid diffusion is given by the continuity equation [23]:
1 M p ˙   +   α . u ˙   +   . v d = 0
v d =     k μ p p
where, M and α are respectively the Biot’s modulus and the Biot’s coefficient, u is the displacement of the solid phase, m; p is the pore fluid pressure, Pa; v d is the fluid flow velocity of the pore fluid, m / s ; k is the permeability tensor, m 2 ; μ p is the viscosity of the pore fluid, Pa. s; p is the gradient of the fracturing fluid pressure along the fracture surface, Pa/m.

3.2. The Porous Medium Deformation

The formation is expected to be isotropic and homogeneous. The porous medium deformation is defined with the Biot’s theory using the following equations [22]:
σ i j     σ i j , 0 = 2 G ε i j   +   ( K     2 3 G ) ε k k δ i j     ( α     1 ) ( p     p o ) δ i j
where, σ i j and σ i j , 0 are respectively the effective stress and the initial effective stress, Pa; G is the shear modulus, Pa; ε i j and ε k k are respectively the strain tensor and the volumetric strain tensor; δ i j is the Kronecker delta, K is the bulk modulus, Pa; α is the Biot’s coefficient, p and p o are respectively the pore pressure and the initial pore pressure, Pa.

4. Model Verification

In the previous section, theories and definitions for the numerical modeling of hydraulically driven fracture problem has been presented. Fluid driven fracture problems are governed by two competing energy dissipation mechanisms associated respectively with viscous fluid flow and fracture propagation; and two competing fluid balance mechanisms associated with fluid storage within the fracture and fluid leakage from the fracture into the surrounding material. Consequently, four limiting propagation regimes as illustrated in Figure 2 (each regime is characterized by the dominance of one dissipation mechanism and one fluid balance mechanism) define the problem of hydraulically driven fracture:
  • Storage-viscosity dominated regime;
  • Leak off-viscosity dominated regime;
  • Storage-toughness dominated regime;
  • Leak off-toughness dominated regime.
The dimensionless parameter M k which is the ratio of energy dissipated during the viscous fluid flow to energy expanded in fracturing the rock is introduced to classify the two energy dissipation mechanisms. For example, in the leak- off toughness dominated regime, the fluid storage inside the fracture is negligible compared to the fluid leak off and the parameter M k 1 .
Despite the fact strong simplifying assumptions were made during the numerical modeling, there are no available closed-form analytical solutions for the problem of fluid driven fracture coupling drill string vibration cyclic loads, fluid flow in porous medium and fracturing fluid leaking into the surrounding rock materials. However, several researchers [24,25,26,27,28,29] have developed small- and large-time asymptotic solutions to predict the fracture development and the net fluid pressure for a hydraulic fracture propagating in an elastic rock driven by a Newtonian fluid. They demonstrated that the asymptotic solutions of the hydraulic driven fracture problem are function of the four limiting propagating regimes.
This research investigates the problem in the storage-toughness dominated regime (near-K). The boundary conditions, material parameters and loads utilized during the numerical modeling are always used to determine the near-K asymptotic solutions. Then, comparisons between Abaqus numerical solutions and analytical asymptotic solutions are conducted to validate the numerical model. The input parameters in Table 1 were specifically selected in order to render the Abaqus numerical results comparable with the asymptotic solutions. For example, fracture geometries are much smaller than the domain dimensions, then, cohesive properties are selected to maintain fracture size larger than cohesive zone, finally, the permeability is defined to reduce the impact of poroelastic ahead of the fracture tip. In the following sections, the different equations and theories governing the small-and large-time asymptotic analytical models are presented.

4.1. Problem Formulation

An incompressible Newtonian fluid with viscosity μ is injected at the center of the fracture at constant flow rate Q o in a brittle elastic fractured formation characterized by Young’s modulus E, Poisson’s ratio ν and fracture toughness K l c . C L is the fluid leak-off coefficient. The crack is loaded by the fluid pressure p f ( x , t ) and the formation is subjected to far field stress σ o perpendicular to the fracture plane. The assumptions utilized during the numerical modeling of the problem are reconducted for the establishment of the analytical model:
  • The rock formation is homogeneous (uniform values of Young’s modulus E, Poison’s ratio ν and fracture toughness K l c );
  • The fracture is designed with radially symmetric geometries;
  • The lag between the fracture front and the fracturing fluid is negligible;
  • Laminar and unidirectional regime govern the flow of incompressible fluid inside the fracture.
Under the above assumptions, the problem is solved by determining the net fluid pressure   p ( x , t ) (difference between fluid pressure inside the crack p f ( x , t ) and far field stress σ o ), the fracture aperture w ( x , t   ) and the fracture half length l ( t ) ; where x is the spatial coordinate and t is the time. The following parameters are first defined to simplify the upcoming governing equations [28]:
E = E 1     ν 2
μ = 12 μ
K = 4 ( 2 π ) 1 / 2 K l c
C = 2 C L
where, E   is the plane strain modulus, Pa; μ ,   K   and   C are respectively the alternate viscosity, Pa.s; fracture toughness, P a . m and leak-off coefficient, m/s/Pa.

4.2. Governing Equations

The solutions w ( x , t ) , p ( x , t ) and l ( t ) of the problem are obtained by solving a set of governing equations summarized as follows:

4.2.1. Global Fluid Conservation

The global fluid conservation assumes that the volume of fluid injected inside the fracture remains equal to the sum of the fracturing fluid volume inside the crack and the fracturing fluid volume leaking into the surrounding rock [25]:
V i n j e c t ( t ) = V c r a c k ( t )   +   V l e a k ( t )
Q o t = 2 0 l ( t ) w ( x , t ) ( π x ) d x + 2 0 t 0 l ( τ ) g ( x , τ ) d x d τ
where,   V i n j e c t ( t ) ,   V c r a c k ( t ) and V l e a k ( t ) are respectively the volume injected inside the fracture, the fracturing fluid volume inside the crack and the fluid leak-off volume, m 3 ; Q o is the fluid flow rate m 3 / s ; w ( x , t ) is the fracture aperture at the time t and the position x, m;   l ( t ) is the fracture length at the time t, m. The function g ( x , t ) denotes the rate of fluid leak-off into the surrounding rock expressed using the Carter’s theory equation:
g ( x , t ) = C t     t o ( x ) , t     t o ( x ) > 0
where, C is the filter cake leak-off coefficient, m . s 1 2 ; t is the injection time, s; t o ( x )     is the time it takes for the fluid front to reach the point x ,   s .

4.2.2. Elasticity Equation

The relationship between the fracture width w and the net pressure p is described with the elasticity equation expressed by [25]:
p ( x , t ) = E 4 π l l w ( s , t ) s d s s x

4.2.3. Lubrication Equation

The flow of incompressible Newtonian fluid inside the fracture is governed by the lubrication equation given as [27]:
w t + g = 1 μ x x ( x w 3 p x ) , t > 0 , x < l ( t )
where, w ,   g   and   p   are respectively the fracture aperture, m; the rate of fluid leak-off, m/s and the net fluid pressure, Pa.

4.2.4. Linear Elastic Fracture Mechanics

The fracture propagation is expressed in the framework of linear elastic fracture mechanics theory with the following equation [27]:
w = 8 π E l 0 l G ( x l , x l ) p ( x , t ) x d x
G ( x l , x l ) = { 1 ξ F ( s i n 1 1 ξ 2 1 ξ 2 , ξ 2 ξ 2 ) , ξ > ξ 1 ξ F ( s i n 1 1 ξ 2 1 ξ 2 , ξ 2 ξ 2 ) , ξ < ξ
where, G is the elastic kernel [30,31] for radial crack geometry. The following boundary conditions which include the assumption that the fracture is always propagating ( K = K l c ), are used to complete the previous equations:
w K E l x , 1 x l 1
w 3 p x = 0 , x = l

4.3. General Scaling

The fracture length l ( t ) , the fracture aperture w ( x , t ) , and the net fluid pressure p ( x , t ) are obtained using the following scaling formulas [26]:
l ( t ) = γ ( τ ) L
w ( x , t ) = ε L Ω ( ξ , τ )
p ( x , t ) = ε E Π ( ξ , τ )
t = t * τ
x = l ( t ) ξ
The small parameter ε , the timescale t * , the length scale L and the dimensionless parameter M k for radial crack geometry are expressed as [27]:
t * = ( K 4 Q o E 4 C 5 ) 2 / 3
L = ( K Q o E C 2 ) 2 / 3
ε = ( C 2 K 2 Q o E 2 ) 1 / 3
M k = μ ( C 4 E 11 Q o K 14 ) 1 / 3
Under the above governing equations and general scaling, the problem of hydraulically driven fracture is resumed to the following nonlinear equations:
1 2 τ = γ m ( τ ) 0 1 Ω ( ξ , τ ) ( π ξ ) m 1 d ξ + ( π 2 ) m 1 0 τ γ m τ τ d τ , τ > 0
Ω τ ξ γ d γ d τ Ω ξ + 1 τ τ o ( x ) = 1 M k γ 2 ξ ξ ( ξ · Ω 3 Π ξ ) , τ > 0 , | ξ | < 1
Ω = 8 γ π 0 1 G ( ξ , ξ ) Π ( ξ , τ ) ξ d ξ , τ > 0 , | ξ | < 1
Ω γ 1 / 2 ( 1 ξ ) 1 / 2 , 1 ξ 1
when M k 1 , and assuming the fracture is always propagating ( K = K l c ), the opening Ω and the pressure Π in Equation (22) can be related to the fracture length γ using fracture elastic mechanics:
Ω = 2 1 / 2 γ 1 / 2 ( 1 ξ 2 ) 1 / 2
Π = ( π 2 ) 2 5 / 2 γ 1 / 2
Finally, the global fluid balance equation governing the scaled fracture length as a function of the dimensionless scaled time is given as:
1 2 τ = π 3 2 γ 7 / 2 + π 2 τ 0 1 γ 2 ( η τ ) 1 η d η

4.4. Asymptotic Analytical Solutions and Comparisons with Abaqus Numerical Solutions

The solution of the global fluid balance equation is supposed to exhibit an asymptotic behavior for small values of expressed as follows [27]:
τ 1 :   γ ( τ ) = τ β i = 0 n γ k i τ i α + O ( τ ( n + 1 ) α + β )
For radial fracture geometry, the initial boundary conditions are given as:
α = 3 10   ; β = 2 5   a n d   γ k 0 = ( 3 2 2 π ) 2 / 5
The other coefficients (for) and the scaled fracture length at small time are determined using the regular perturbation theory:
γ k = γ τ β ,   C k = τ α ,   τ 1
By combining Equations (14)–(17), (27)–(35), (40)–(41) and (43)–(45); the small-time asymptotic solutions w ( x , t ) , p ( x , t ) and l ( t ) are determined and can be compared to the Abaqus numerical solutions. The flowchart of algorithm for the determination of the asymptotic solutions of the hydraulic driven fracture problem is presented in Figure 3.

5. Description and Application of the Numerical Model

5.1. Geometry and Material Properties

The schematic of the numerical model is illustrated in Figure 4. The model consists of the wellbore, the rock formation and the natural fracture. The wellbore radius is 0.1 m. The well starts from the surface and extends up to the depth of 1000 m. The dimensions of the rock formation are 20   ×   120   m . Due to the symmetry of the model; only half of the formation is represented with a 2D geometry. The rock formation is modeled using plane strain coupling pore fluid stress CPE4P elements. The maximum and minimum horizontal stresses applied on the rock formation are respectively in the x-and y-directions.
The natural fracture is initially opened as it can be seen in Figure 1. The natural fracture is designed with the Perkins-Kern-Nordgren (PKN) model with an initial length of 0.5 m and an initial width of 2   ×   10 4   m . Due to the fluctuating wellbore pressure caused by drill string vibration, the natural fracture will re-open and the fracture tip propagates following cohesive zone model in the direction of the maximum horizontal stress. Two-dimensional pore pressure cohesive elements COH2D4P are used to model the cohesive zone elements. The initial pore pressure within the natural fracture is 10 Mpa.
The drill pipe is modeled with 2-node linear fluid pipe connector elements FPC2D2 while the wellbore annulus is modeled with 2-node linear fluid pipe elements FP2D2. The standard connection type model includes constant bidirectional loss terms that allow the accurate prediction of the total pressure loss inside the fluid pipe connector. Therefore, this model is used in this study to compute the fluid loss inside the fluid pipe connector elements. The Churchill model is utilized to compute the fluid loss inside the fluid pipe elements since this model takes into account the pipe roughness and the fluid flow regimes [20,32].

5.2. Drill String Vibration Modeling

The drill string vibration is modeled using a sinusoidal function of time. The mathematical equation of the stress generated by the drill string vibration cyclic loads on the wellbore surface is given by:
σ ( t ) = S   c o s ( ω t + φ )
where, σ ( t ) is the stress exerted by the drill string vibration cyclic loads on the wellbore surface at the time t, Pa; S is the maximal magnitude of σ ( t ) , Pa; φ is the initial phase of the vibration, rad; ω is the angular frequency, rad/s. ω is given by:
ω =   2 π N T
where, N is the number of cyclic loads and T is the total simulation time, s. The parameters φ ,   N ,   S   totally define the stress vibration σ ( t ) . In this research, the initial phase φ of the drill string vibration stress is assumed to be zero. The number of cyclic loads N is related to the drill string RPM by the equation:
N =   ( R P M ) × T 60
The maximal magnitude S of the stress vibration σ ( t ) is related to the WOB according to the following procedure: The drill string applies WOB through the cross-sectional area A of the drill collar. After each shock of the drill string at the bottom hole, vibration cyclic loads are generated and the maximal magnitude S m a x of the stress vibration at the impact point of the bottom hole is obtained by dividing the force due to the WOB by the cross-sectional area A of the drill collar. The maximal magnitude S m a x of the stress vibration at the impact point of the bottom hole progressively decreases with the time and space due to friction losses inside the wellbore. The vibrations cyclic loads finally hit the wellbore surface with a stress magnitude S equal to 60% of the initial maximal magnitude S m a x .
S =   0 . 6 × ( W O B ) × g A
where, g is the gravitational acceleration, m / s 2 .The RPM used in this paper ranges from 1 to 10 and the WOB varies between 6500 lbs to 70,000 lbs. For different values of the drilling operational parameters WOB and RPM, the number of cyclic loads N and the magnitude of the stress vibration S hitting the wellbore surface can be determined respectively using equations (48) and (49), therefore the stress exerted by the drill string vibration cyclic loads on the wellbore surface is totally defined since the initial phase of the vibration φ is assumed to be zero.

5.3. Simulation Steps

The analysis is run in two steps: an initial step where the initial pore pressure, the initial in situ stresses and the initial void ratio are applied on the model. The next step is the geostatic step where the equilibrium of the different loads applied during the initial step is achieved. The next step is a transitional soils consolidation analysis step where the dynamic mud circulation, the drill string vibration cyclic loads and the loss circulation are simulated. An unsymmetric matrix storage is also used in this step to improve the convergence rate of the non-linear solution. The incremental size of this step varies between 1 to 10 s and the total step time is t = 100 s.

5.4. Loading and Boundary Conditions

The drilling fluid is injected in the wellbore from the node of the drill pipe lying at the surface. A periodic distributed surface load is applied on the wellbore surface to simulate the drill string vibration cyclic loads. The wellhead pressure is defined on the node of the wellbore annulus lying at the surface with a zero-boundary pore fluid pressure. The ‘TIE’ constraint is used to attach the node of the drill pipe lying at the bottom hole with the nodes of the formation located at the natural fracture tip and the circumference of the wellbore formation. The ‘PORMECH’ constraint is equally used to apply the bottom hole pressure at the wellbore annulus onto the wellbore surface as a mechanical surface pressure. The left edge is modeled with the X symmetry boundary condition and a zero-pore fluid pressure. The three other edges are modeled with constrained normal displacements and pore pressure equals to 10 MPa.

6. Presentation of Data and Results

In this section, numerical results are presented to demonstrate the importance of the drill string vibration cyclic loads on the development of the wellbore natural fracture. The input parameters used during this simulation are the same parameters used by Feng and Gray (2018) in their research. However, the cyclic load magnitude (S) and the number of cyclic loads to failure (N) have been added to model the drill string vibration cyclic loads. The cohesion and the friction angle have also been considered to define the plastic behavior of the rock formation. The input parameters corresponding to a sandstone type rock are presented in Table 1.

6.1. Establishment of the S-N Curves

The Mohr Coulomb plasticity model is used to investigate the post elastic behavior of the rock. The shear failure of the rock is assumed to occur when the Plastic Element Equivalent strain (PEEQ) is greater than zero. The detailed algorithm of the establishment of the S-N curve is given in Figure 5.
By applying the different steps of the flow chart given in Figure 5, using the parameters provided in Table 1, different number of cyclic loads to failure ( N o ) are obtained for successive decreasing values of cyclic load magnitudes ( S o ) . The results obtained are plotted in the S-N curve presented in Figure 6. It is observed that the cyclic load magnitude (S) decreased with the number of cyclic loads to failure (N). Each cyclic load generates a decohesion of the matrix formation. Therefore, the greater the magnitude of the cyclic load, the smaller the number of cyclic loads required to fail the rock. These results are in agreement with the results of other researchers [8,33,34].
It is equally noticed that the S-N curve encompasses three regions: the plastic region, the elastic region and the infinite life region. In the infinite life region, the natural rock formation did not fail even if the rock formation is subjected to an infinite number of cyclic loads. The elastic region is the region where the material recovers its initial position after the cyclic load is released. In the plastic region, the rock material undergoes permanent deformation. In this region, few cyclic loads with high stress amplitude are sufficient to fail the rock formation. The different regions of the S-N curve are separated with three threshold values. The ultimate strength of the rock (28.85 MPa) is the magnitude of the cyclic load for which the rock fails after one cyclic load. The yield strength (23.5 MPa) is at the limit between the plastic region and the elastic region. Below a specific magnitude of the cyclic load called the endurance limit of the rock (11.3 MPa), the rock did not fail even if the number of cyclic loads is increased.
Numerous authors in the literature have assumed that the endurance limit of the rock is reached for a normalized ratio of the magnitude of the cyclic load σ m a x to the monotonic strength σ m i n equals to 0.7 [15]. In this paper, it can be observed that the endurance limit of the rock formation is reached for a normalized ratio equals to 0.33 (the monotonic strength of the rock equals 34.62 MPa). In the literature, most of the S-N curves were established using laboratory experiments. During the laboratory experiments, the rock specimen is only subjected to the confining pressure and the axial deformation pressure. In the numerical simulation conducted in this paper, apart from the confining pressure and the axial deformation pressure, loss circulation which occurs at the fracture wall is also taken into consideration using the pipe elements theory. This loss circulation globally weakens the rock formation and diminishes the strength of the rock formation. After reduction of the strength, the magnitude of the cyclic load necessary to fail the rock is also reduced and finally the normalized ratio is therefore smaller than the ratio present in the literature.

6.2. Development of the Wellbore Natural Fracture for Different Cyclic Load Magnitude and Constant Number of Cylic Loads

In this section, the fracture width and fracture length evolution in function of the time for different cyclic load magnitude (S) and constant number of cyclic loads (N) is investigated. It is observed in Figure 7, that the fracture width on each profile increases with the cyclic load magnitude. The maximal value of the fracture width without vibration cyclic loads increased by 64.77% after application of cyclic loads with parameters S = 11.3 MPa and N = 10. Equally, the profile of the fracture width without vibration is almost smooth, then, for different increasing cyclic load magnitude, the fracture width as a function of time profiles follows a sinusoidal behavior similar to the drill string vibration cyclic loads profile. The drill string vibration cyclic loads which are hitting the wellbore surface generate a decohesion in the preexisting wellbore natural fracture, then, the trend of the opening and closure of the wellbore natural fracture in function of the time also follows the trend of the drill string vibration oscillating profiles. Other studies in the literature have also highlighted the decohesion of the rock structure caused by the drill string vibration cyclic load [8,12,15].
The fracture width evolution as a function of the fracture length for different cyclic load magnitude and constant number of cyclic loads is presented in Figure 8. The results show that the values of the fracture width decreased with the increase of the fracture length due to the leakage of the fracturing fluid from the fracture into the surrounding rock. Equally, for different increasing cyclic load magnitude and constant number of cyclic loads: (i) The fracture width increased with fracture length in the near wellbore region (fracture length less than 5.16 m; (ii) For fracture length between 5.16 m and 22.83 m, the fracture width decreased with the fracture length and (iii) for fracture length greater than 22.83 m, the profiles of the fracture width in function of the fracture length remain the same.
The fluid pressure inside the fracture is highest near the wellbore, therefore, a larger fracture width in that region is also expected. Nevertheless, stress concentration around the wellbore is generally observed during drilling operations due to the excavation of the rock from the formation. Although, the drilling mud is utilized to redistribute stresses around the wellbore and ensure wellbore stability, this redistribution is not instantaneous. Since the time simulation of this modeling is relatively small (t = 100 s), therefore stress concentration is always observed around the wellbore at the final simulation stage. The stress concentration in the near wellbore region acts against the opening of the fracture resulting in smaller fracture width. The effect of stress concentration on the fracture width is alleviated after integration of different increasing cyclic loads because the cyclic loads buffeting the wellbore surface generate a decohesion of the near wellbore fracture leading to a significant increase of the fracture width with the cyclic loads.

6.3. Development of the Wellbore Natural Fracture for Constant Cyclic Load Magnitude and Different Number of Cyclic Loads

The development of the wellbore natural fracture for constant cyclic load magnitude and different number of cyclic loads is presented in this section. The results of the fracture width evolution in function of the time are shown in Figure 9. It can be observed that the fracture width as a function of time profile without vibration is almost smooth, then for different increasing values of cyclic loads number, the fracture width as a function of time profiles follows a sinusoidal behavior similar to the drill string vibration cyclic loads profile. The fracture width without vibration cyclic loads also increased after application of the drill string vibration cyclic loads (an increase of 63.12% of the fracture width without vibration was observed after application of cyclic loads with parameters S = 11.3 MPa and N = 15). However, contrary to the case of the evolution of the fracture width in function of the time for different cyclic load magnitude and constant number of cyclic loads where the maximal values of the fracture width on each profile were increasing with the cyclic load magnitude, it is observed in this case that the maximal values of the fracture width on each profile are almost the same when the cyclic load magnitude is constant but it is the oscillating period of the fracture width which increases due to the increase of the number of cyclic loads.
The fracture width evolution as a function of the fracture length for different number of cyclic loads and constant cyclic load magnitude is presented in Figure 10. The results show that the fracture width profiles are almost the same for different increasing number of cyclic loads and constant cyclic load magnitude. The fracture width decreased with the increase of the fracture length due to the seepage of the fracturing fluid from the fracture into the rock. Equally, the fracture width increased after implementation of drill string vibration with constant cyclic load magnitude and different number of cyclic loads.
At the wellbore surface, the fracture width without vibration cyclic loads increased by 70.96% after application of cyclic loads with parameters S = 11.3 MPa and N = 15.
The fracture length development as a function of time is a major aspect which has been investigated in the research. Vibration cyclic loads were integrated in the system, then, fracture length was determined at different simulation time. The fracture length evolution with the time obtained after implementation of cyclic load in the model is presented in Figure 11. The results show that the fracture length with and without vibration cyclic loads are slightly the same because the effect of cyclic load is localized only on the wellbore surface. This result demonstrates that the cyclic loads which are hitting the wellbore surface weakly affect the fracture length development. However, it should be noted that these results are obtained with a relatively small simulation time (t = 100 s). Further studies need to be conducted to extrapolate these results with larger simulation time.

7. Discussion and Parametric Studies

In this section, parametric studies are presented to investigate the effect of other drilling operational parameters on the development of the wellbore natural fracture in considering the effect of the drill string vibration cyclic loads. The fracture width, the bottom hole pressure and the loss circulation are all measured at the fracture mouth. All the results are presented at the final stage of simulation at t = 100 s. Discussion and parametric studies are conducted by comparing the results obtained in this paper (poro elasto plastic model coupled with drill string vibration cyclic load) with the results of the poroelastic model without drill string vibration cyclic loads [17].

7.1. Effect of Mud Density on Fracture Development Integrating Drill String Vibration

The mud density is an important drilling operational parameter which affects the wellbore stability. The results of the effect of mud density on fracture length and fracture width with and without vibration are shown in Figure 12. The results show that for different increasing drilling mud densities, the fracture width obtained with drill string vibration are greater than those obtained without vibration due to the enlargement of the fracture width caused by the drill string vibration cyclic loads. The fracture width for example increased by 64.77% after integration of cyclic loads with parameters S = 11.3 MPa, N = 10 and mud density equal to 1.3 g / cm 3 . Apart from mud density being 1.2 g / cm 3 , the fracture length with and without vibration cyclic loads are slightly the same because the effect of the cyclic loads is more localized at the wellbore surface. Therefore, it is the width of the fracture mouth which is greatly affected by the cyclic loads, the fracture length itself does not increase too much after integrating the drill string vibration in the system.
For mud densities between 1.1 g / cm 3 and 1.2 g / cm 3 , the fracture length and fracture width without vibration remain equals to zero which implies that for this range of mud density, the bottom hole pressure (BHP) is not sufficient to initiate the fracture propagation. However, for this same range of mud density, the model with drill string vibration already predicts the fracture propagation because of the cyclic loads which are hitting the wellbore surface. The fracture length for example increases from zero to 12.54 m after integration of cyclic loads with parameters S = 11.3 MPa and N = 10 and mud density equal to 1.2 g / cm 3 . These results point out the importance of the newly built model that can accurately predict the fracture growth development. In fact, the model without vibration tends to underestimate the fracture width and fracture length development. This underestimation can lead to several problems such as the inaccurate selection of the particle size dimension (PSD) necessary to plug the natural fracture and maintain the wellbore stability during drilling operations.
The effect of mud density on the loss circulation rate and the BHP with drill string vibration cyclic loads is also investigated and the results are presented in Figure 13. It is straightforward that an increase of mud density will lead to an increase of the hydrostatic mud pressure and consequently an increase of the BHP. For different increasing drilling mud density, the BHP remains slightly the same with or without vibration cyclic loads. In fact, the BHP is mainly affected by the hydrostatic mud pressure and the viscous pressure losses. Therefore, the cyclic loads which are applied on the wellbore natural fracture do not greatly influence the value of the BHP. For mud density between 1.1 g / cm 3 and 1.2 g / cm 3 , the loss circulation rate without vibration cyclic loads is tiny because the model without vibration could not predict fracture initiation and propagation for this range of density.
For mud densities between 1.2 g / cm 3 and 1.4 g / cm 3 , the fracture initiation pressure has already been reached and the loss circulation rate without vibration starts increasing with the mud density due to the augmentation of the BHP in the system. The model integrating vibration cyclic loads predicted fracture initiation and propagation for the whole range of mud densities between 1.1 g / cm 3 and 1.4 g / cm 3 , therefore the loss circulation rates increase with mud density for this range of density due to the increase of the BHP in the system. The loss circulation rate obtained with drill string vibration are greater than the results obtained without vibration cyclic loads. An increase of 14.3% of the loss circulation rate was observed after integration of cyclic loads with parameters S = 11.3 MPa and N = 10 and mud density equal to 1.3 g / cm 3 . This result was expected because each cyclic load applied on the natural fracture progressively enlarges the fracture width. The natural fracture becoming larger and larger will lead to an increased loss circulation in the system.

7.2. Effect of Mud Viscosity on the Fracture Development Integrating Drill String Vibration

The mud viscosity is another major operational drilling parameter which affects the wellbore stability. The impact of the mud viscosity on fracture length and fracture width with and without vibration cyclic loads is investigated and presented in Figure 14. The results show that the fracture width with and without vibration cyclic loads increase with the mud viscosity because the higher the fluid viscosity, the larger the resistance of the fluid flow into the fracture and finally the larger the fracture width. However, the fracture length is not greatly affected by the mud viscosity increase. For mud viscosity less than 10 cp, the fracture length linearly increases with the mud viscosity up to reach 30.3 m, and finally stabilizes around this value for mud viscosity between 10 cp and 100 cp. Those results are consistent with hydraulic practices theories which stipulate that higher fluids viscosity lead to wider and shorter fracture [17,35].
Concerning the influence of the drill string vibration cyclic loads on the natural fracture development, it is observed that the fracture width is more affected by the drill string vibration cyclic loads because the cyclic loads hitting the wellbore surface cause the decohesion of the formation matrix in the neighborhood of the fracture mouth. For mud viscosity equal to 1 cp for example, the fracture width without vibration increased by 64.78% after integration of drill string vibration cyclic loads with parameters S = 11.3 MPa and N = 10. However, like the investigations of the effect of mud density on the fracture length development with and without vibration, the integration of the drill string vibration cyclic loads in the system does not significantly affect the fracture length development because its impact is more localized at the wellbore surface.
It is also interesting to study the effect of mud viscosity on the loss circulation rate and the BHP with and without vibration cyclic loads. From Figure 15, it is observed that the loss circulation rate and the BHP increase with the mud viscosity because the increase of the mud viscosity will result in the elevation of the viscous pressure losses which will in turn lead to the augmentation of the BHP. Importantly, the increase of the BHP also enhances the risk of loss circulation which is materialized in the system by an increase of the loss circulation rate.
The results also showed that the BHP values do not change too much when drill string vibration cyclic loads are considered in the simulation. The maximal percentage of the BHP increase after integration of drill string vibration did not exceed 0.5%. However, the loss circulation rates increase when drill string vibration cyclic loads are considered in the simulation. The loss circulation rate increased by 14.3% after drill string vibration cyclic loads with parameters S = 11.3 MPa and N = 10 was added in the system. Some researchers such as [35,36] arrived at contrasting conclusions when they demonstrated that the higher the viscosity of the drilling mud, the lower the loss circulation rate inside the fracture. In this paper, the dimensions of the wellbore (depth, length and width) are far greater than the dimensions of the wellbore natural fracture, therefore it is the BHP which is the key driver of the loss circulation inside the fracture, so when the mud viscosity will increase, the viscous pressure loss will also increase and it will finally lead to an increase of the BHP and the loss circulation rate.

7.3. Effect of Pump Rate on the Fracture Development Integrating Drill String Vibration

The pump rate is one of the dominating factors that control the fluid flow velocity inside the wellbore. The effect of pump rate on fracture length and fracture width with and without vibration cyclic loads is presented in Figure 16. The results showed that the fracture width and fracture length increase with pump rate due to the augmentation of the wellbore pressure. The fracture width obtained when the drill string vibration cyclic loads are implemented in the system are greater than those obtained without vibration cyclic loads. An average increase of 64.94% of the fracture width is observed when the drill string vibration are considered during the simulation. It is equally observed that the fracture length did not evolve too much after integration of vibration cyclic load in the system because the impact of cyclic loads is localized on the wellbore surface.
The effect of pump rate on loss circulation rate and the BHP when the drill string vibration are integrated in the system have equally been investigated. Figure 17 shows that the loss circulation rate and the BHP increase with the pump rate. The loss circulation rate obtained with drill string vibration cyclic loads are greater than those obtained without vibration cyclic loads. These results confirm that the tripping operations should be made very slowly in order to avoid high values of pump rates. This will enable to prevent uncontrolled fracture growth development and loss circulation.
Figure 18, Figure 19 and Figure 20 present the asymptotic solutions of time evolution of fracture aperture, fracture length and fluid pressure inside the crack. The asymptotic solutions of the pressure and fracture distribution along the fracture are also presented in Figure 21 and Figure 22. The asymptotic solutions were obtained by following the steps of the flowchart in Figure 3. Good agreement between asymptotic solutions and numerical ABAQUS solutions is observed which allows to validate the numerical model. The numerical results and conclusions obtained in this study can therefore be trusted due to the establishment of the verification analytical model.

8. Conclusions

In this research, a poro-elasto-plastic model based on a finite element method has been established to analyze the influence of the drill string vibration cyclic loads on the development of the wellbore natural fracture. The flow of the fluid inside the wellbore is modeled with the pipe element theory in ABAQUS. Drill string vibration is modeled with a sinusoidal function of time characterized by the magnitude of the stress vibration S and the number of cyclic loads N. The ABAQUS coupled pressure/deformation cohesive zone model is used to simulate the fracture propagation and the fracturing fluid flow, while the pore fluid flow inside the rock formation and the porous medium deformation is described with coupled pore-pressure/deformation continuum finite elements model. The main results of the research showed that:
(1)
The S-N curve for the formation is subdivided into three main regions: a plastic region with an ultimate strength of 28.85 MPa, the elastic region with a yield strength of 23.5 MPa and an infinite life region with an endurance limit 11.3 MPa;
(2)
The fracture width as a function of the time profiles follows a sinusoidal behavior similar to the drill string vibration cyclic loads profile; For different increasing cyclic load magnitude S with constant number of cyclic loads N, the maximal values of the fracture width on each fracture profile increased with the cyclic load magnitude; the fracture width increased by 64.77% after application of cyclic loads with parameters S = 11.3 MPa and N = 10;
(3)
For constant cyclic load magnitude with different number of cyclic loads, an increase of 63.12% of the fracture width is observed after application of cyclic loads with parameters S = 11.3 MPa and N = 15. The maximal values of the fracture width on each fracture profile are almost the same and the oscillating period of the fracture width increases with the number of cyclic loads;
(4)
For different increasing cyclic load magnitude and constant number of cyclic loads: (i) The fracture width increased with fracture length in the near wellbore region (fracture length less than 5.16 m; (ii) For fracture lengths between 5.16 m and 22.83 m, the fracture width decreased with the fracture length and (iii) for fracture lengths greater than 22.83 m, the profiles of the fracture width in function of the fracture length remain the same;
(5)
The fracture width profiles in function of the fracture length are almost the same for different increasing number of cyclic loads and constant cyclic load magnitude;
(6)
The drill string vibration cyclic loads which are hitting the wellbore surface do not greatly affect the fracture length and the bottom hole pressure; However, the fracture width and the loss circulation rates are considerably impacted by the drill string vibration cyclic loads; The fracture width and loss circulation rate increased by 64.77% and 14.3%, respectively, after integration of cyclic loads with parameters S = 11.3 MPa, N = 10 and mud density equal to 1.3 g / cm 3 .
The newly built model can accurately predict the natural fracture growth when the wellbore pressure is subjected to fluctuations caused by drill string vibration cyclic loads. The results obtained in this research can be used for wellbore strengthening studies to determine the optimal PSD necessary to plug the natural fractures and mitigate the lost circulation during drilling operations.

Author Contributions

Conceptualization: A.R.K.L. and J.D.; methodology: A.R.K.L., Y.F. and H.L.; software and Investigation: A.R.K.L. and Y.F.; resources and data curation: A.R.K.L. and H.L.; writing—original draft preparation: A.R.K.L.; N.B.S.S. and Y.S.; writing—review and editing: A.O. and M.M.; supervision: J.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Science Foundation of China University of Petroleum, Beijing, grant number 2462019YJRC008 and National Natural Science Foundation of China, grant number 52074312. The APC was funded by 2462019YJRC008.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The study did not report any data.

Acknowledgments

This research was supported by China University of Petroleum (Beijing) through China Scholarship Council (CSC).

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

A Drill Collar Cross-Sectional Area, m 2
c T , c B Leak off coefficients respectively at the top and the bottom of the fracture, m/s/Pa
c l Fluid loss coefficient
B H P Bottom Hole Pressure
g Gravitational acceleration, m / s 2
k Permeability tensor, m 2
G Shear modulus, Pa
K Bulk modulus, Pa
M Biot’s modulus
NNumber of cyclic loads to failure
p , p o Pore pressure and initial pore pressure respectively, Pa
P K N Perkins-Kern-Nordgren
P S D Particle Size Dimension
p f Fracturing fluid pressure along the surface of the fracture, Pa
P 1 , P 2 Pipe element nodes
p T , p B Pore fluid pressures respectively at the top and the bottom of the fracture, Pa
Δ P l o s s Viscous pressure loss, Pa
Δ P Pressure difference between P 1 and P 2 , Pa
Δ P h y d r o Hydrostatic drilling fluid pressure between P 1 and P 2 , Pa
p     Gradient of the fracturing fluid pressure along the fracture surface, Pa/m
q f Longitudinal flow rate per unit of width, m 2 / s
RPMRevolution Per Minute
SMaximal magnitude of the stress exerted by the drill string vibration cyclic loads on the wellbore surface at the time t, Pa
S m a x Maximal magnitude of the stress exerted by the drill string vibration cyclic loads at the impact point of the bottom hole, Pa
TTotal simulation time, s
t n ,   t s 1 , t s 2 Normal nominal stress, First and second shear nominal stress respectively, Pa
t n 0 , t s 1 0 , t s 2 0 Normal nominal stress, First and second shear nominal stress peak values respectively, Pa
uDisplacement of the solid phase, m
v t , v b Normal velocities respectively at the top and the bottom of the fracture, m/s.
v d Fluid flow velocity of the pore fluid, m/s
V Fluid velocity, m/s
w Fracture aperture, m
WOBWeight on bit
Δ Z Elevation difference, m
α Biot’s coefficient
φ Initial phase of drill string vibration, rad
ε i j , ε k k Strain tensor and volumetric strain tensor respectively
ω Angular frequency of drill string vibration, rad/s
σ i j , , σ i j , 0 Effective stress and initial effective stress respectively, Pa
σ ( t ) Stress exerted by the drill string vibration cyclic loads on the wellbore surface at the time t, Pa
δ i j Kronecker delta
ρ Drilling fluid density, kg / m 3
μ f Viscosity of the fracturing fluid, Pa·s.
μ p Viscosity of the pore fluid, Pa·s.

SI Metric Conversion Factors

cp×1.0 *10−3=Pa·s
in.×2.54 *10−2=m
lbm×0.45359237100=kg
md×9.86923310−16=m2
* conversion factor is exact.

References

  1. Meng, M.; Zamanipour, Z.; Miska, S.; Yu, M.; Ozbayoglu, E. Dynamic stress distribution around the wellbore influenced by surge/swab pressure. J. Pet. Sci. Eng. 2019, 172, 1077–1091. [Google Scholar] [CrossRef]
  2. Jaeger, J.C.; Cook, N.G.; Zimmerman, R. Fundamentals of Rock Mechanics; John Wiley&Sons: Hoboken, NJ, USA, 2009. [Google Scholar]
  3. Fjar, E.; Holt, R.M.; Raaen, A.M.; Risnes, R.; Horsrud, P. Petroleum Related Rock Mechanics, 2nd ed.; Elsevier: Amsterdam, The Netherlands, 2008. [Google Scholar]
  4. Zhang, F.; Kang, Y.; Wang, Z.; Miska, S.; Yu, M.; Zamanipour, Z. Real-Time Wellbore Stability Evaluation for Deepwater Drilling During Tripping. In Proceedings of the SPE Deepwater Drilling and Completions Conference, Galveston, TX, USA, 14–15 September 2016. [Google Scholar]
  5. Zamanipour, Z.; Miska, S.Z.; Hariharan, P.R. Effect of Transient Surge Pressure on Stress Distribution around Directional Wellbores. In Proceedings of the IADC/SPE Drilling Conference and Exhibition, Fort Worth, TX, USA, 1–3 March 2016; p. 178831. [Google Scholar]
  6. Dunayevsky, V.; Abbassian, F.; Judzis, A. Dynamic Stability of Drillstrings Under Fluctuating Weight on Bit. SPE Drill. Complet. 1993, 8, 84–92. [Google Scholar] [CrossRef]
  7. Santos, H.; Placido, J.; Wolter, C. Consequences and Relevance of Drillstring Vibration on Wellbore Stability. In Proceedings of the SPE/IADC Drilling Conference, Amsterdam, The Netherlands, 9–11 March 1999; Society of Petroleum Engineers (SPE): Danbury, CT, USA, 1999; p. 52820. [Google Scholar]
  8. Ewy, R.T.; Bovberg, C.A.; Chen, G.; Jansson, R.; Pepin, G. Fatigue Testing of Hollow Cylinders and Application to Injection Well Cycling. In Proceedings of the Gulf Rocks 2004, the 6th North America Rock Mechanics Symposium, Houston, TX, USA, 5–9 June 2004. [Google Scholar]
  9. Vargas, E.A., Jr.; Chavez, E.; Gusmão, L.; Amaral, C. Is thermal fatigue a possible mechanism for failures of some rocks slopes in Rio de Janeiro, Brazil? In Proceedings of the 43th US Rock Mechanics Symposium, Asheville, NC, USA, 28 June–1 July 2009. [Google Scholar]
  10. Zhu, X.; Liu, W. The effects of drill string impacts on wellbore stability. J. Pet. Sci. Eng. 2013, 109, 217–229. [Google Scholar] [CrossRef]
  11. Esmaeili, A.; Elahifar, B.; Fruhwirth, R.K.; Thonhauser, G. Effect of Formations Compressive Strength on Drill String Vibrations. In Proceedings of the IPTC 2013: International Petroleum Technology Conference, Beijing, China, 26–28 March 2013. [Google Scholar]
  12. Khaled, M.S.; Shokir, E.M. Effect of Drillstring Vibration Cyclic Loads on Wellbore Stability. In Proceedings of the SPE Middle East Oil & Gas Show and Conference, Manama, Bahrain, 6–9 March 2017; Society of Petroleum Engineers (SPE): Danbury, CT, USA, 2017. [Google Scholar]
  13. Osogba, O.J.; Otamere, B.; Abdullahi, G.S.B. Analysis of Vibrations in Drilling Systems. In Proceedings of the SPE Nigeria Annual International Conference and Exhibition, Lagos, Nigeria, 6–8 August 2018; Society of Petroleum Engineers (SPE): Danbury, CT, USA, 2018. [Google Scholar]
  14. Zhao, B.; Liu, D.; Li, Z.; Huang, W.; Dong, Q. Mechanical Behavior of Shale Rock under Uniaxial Cyclic Loading and Unloading Condition. Adv. Civ. Eng. 2018, 2018, 1–8. [Google Scholar] [CrossRef] [Green Version]
  15. Cerfontaine, B.; Collin, F. Cyclic and Fatigue Behaviour of Rock Materials: Review, Interpretation and Research Perspectives. Rock Mech. Rock Eng. 2018, 51, 391–414. [Google Scholar] [CrossRef]
  16. Well Control for the Rig-Site Drilling Team; Aberdeen Drilling School &Well Control Training Centre: Aberdeenm, UK, 2002.
  17. Feng, Y.; Gray, K.E. Modeling Lost Circulation Through Drilling-Induced Fractures. SPE J. 2018, 23, 205–223. [Google Scholar] [CrossRef]
  18. Liu, W.; Zhou, B.; Li, Y.; Yu, B.; Deng, J. Numerical Analysis of Wellbore Instability in Weakly Consolidated Rock Formations. In Proceedings of the ISRM International Symposium- 10th Asian Rock Mechanics Symposium, Singapore, 29 October–3 November 2018. [Google Scholar]
  19. Fekete, P.; Bruno, L.A.; Dosunmu, A.; Odagme, S.; Sanusi, A.; Bowe, E. The effect of wellbore stability in naturally fractured reservoirs. In Proceedings of the Effect of Wellbore Stability in Naturally Fractured Reservoirs, Lagos, Nigeria, 4–6 August 2015; Society of Petroleum Engineers (SPE): Danbury, CT, USA, 2015; p. 178267. [Google Scholar]
  20. SIMULIA. Abaqus Version 2016 Analysis User’s Guide; Dassault Systems: Providence, RI, USA, 2016. [Google Scholar]
  21. Wu, R.; Deng, J.; Liu, W.; Mao, S.; Sun, J.; Yan, M.; Li, M.; Li, Y. Numerical investigation of fluid injection into poorly consolidated geomaterial considering shear dilation and compaction. J. Pet. Sci. Eng. 2018, 168, 119–132. [Google Scholar] [CrossRef]
  22. Zielonka, M.G.; Searles, K.H.; Ning, J.; Buechler, S.R. Development and Validation of Fully-Coupled Hydraulic Fracturing Simulation Capabilities. In Proceedings of the 2014 SIMULIA Community Conference, Providence, RI, USA, 22 May 2014. [Google Scholar]
  23. Feng, Y.; Gray, K.E. Effects of Porous Properties of Rock on Near-Wellbore Hydraulic Fracture Complexity. In Proceedings of the 4th Unconventional Resources Technology Conference, Houston, TX, USA, 23–25 July 2018; American Association of Petroleum Geologists: Tulsa, OK, USA, 2018. [Google Scholar]
  24. Savitski, A.; Detournay, E. Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: Asymptotic solutions. Int. J. Solids Struct. 2002, 39, 6311–6337. [Google Scholar] [CrossRef]
  25. Detournay, E.; Adachi, J.I.; Garagash, D.I.; Savitski, A.A. Interpretation and Design of Hydraulic Fracturing Treatments. U.S. Patent Application No. 7111681, 26 September 2006. [Google Scholar]
  26. Garagash, D.I. Plane-strain propagation of a fluid-driven fracture during injection and shut-in: Asymptotics of large toughness. Eng. Fract. Mech. 2006, 73, 456–481. [Google Scholar] [CrossRef]
  27. Bunger, A.P.; Detournay, E.; Garagash, D.I. Toughness-dominated Hydraulic Fracture with Leak-off. Int. J. Fract. 2005, 134, 175–190. [Google Scholar] [CrossRef]
  28. Peirce, A.; Detournay, E. An implicit level set method for modeling hydraulically driven fractures. Comput. Methods Appl. Mech. Eng. 2008, 197, 2858–2885. [Google Scholar] [CrossRef]
  29. Hu, J.; Garagash, D.I. Plane-Strain Propagation of a Fluid-Driven Crack in a Permeable Rock with Fracture Toughness. J. Eng. Mech. 2010, 136, 1152–1166. [Google Scholar] [CrossRef]
  30. Sneddon, I.N.; Lowengrub, M. Crack Problems in the Classical Theory of Elasticity; John Wiley & Sons: New York, NY, USA, 1969. [Google Scholar]
  31. Barr, D.T. Leading-Edge Analysis for Correct Simulation of Interface Separation and Hydraulic Fracturing. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, USA, 1991. [Google Scholar]
  32. Churchill, S.W. Friction Factor Equation Spans All Fluid-Flow Regimes. Chem. Eng. 1977, 84, 91–92. [Google Scholar]
  33. Erarslan, N.; Williams, D. Mechanism of rock fatigue damage in terms of fracturing modes. Int. J. Fatigue 2012, 43, 76–89. [Google Scholar] [CrossRef]
  34. Momeni, A.; Karakus, M.; Khanlari, G.; Heidari, M. Effects of cyclic loading on the mechanical properties of a granite. Int. J. Rock Mech. Min. Sci. 2015, 77, 89–96. [Google Scholar] [CrossRef]
  35. Yew, C.H.; Weng, X. Mechanics of Hydraulic Fracturing; Gulf Professional Publishing: Houston, TX, USA, 2014. [Google Scholar]
  36. Majidi, R.; Miska, S.Z.; Yu, M.; Thompson, L.G.; Zhang, J. Quantitative Analysis of Mud Losses in Naturally Fractured Reservoirs: The Effect of Rheology. SPE Drill. Complet. 2010, 25, 509–517. [Google Scholar] [CrossRef]
Figure 1. Illustration of the fluid flow inside the wellbore.
Figure 1. Illustration of the fluid flow inside the wellbore.
Energies 14 02015 g001
Figure 2. Parametric diagram of the four limiting propagating regimes of the hydraulically driven fracture problem.
Figure 2. Parametric diagram of the four limiting propagating regimes of the hydraulically driven fracture problem.
Energies 14 02015 g002
Figure 3. Flowchart algorithm for the determination of asymptotic solutions.
Figure 3. Flowchart algorithm for the determination of asymptotic solutions.
Energies 14 02015 g003
Figure 4. Schematic of the numerical model redrawn from [17]. S h m a x and S h m i n respectively represent the maximal and minimal horizontal stresses. The drill string and the wellbore annulus are in the vertical Z direction, but they are represented in the Y direction for better visualization. The natural fracture opens in the horizontal X-Y plane. The formation is in plane strain condition in the horizontal X-Y plane.
Figure 4. Schematic of the numerical model redrawn from [17]. S h m a x and S h m i n respectively represent the maximal and minimal horizontal stresses. The drill string and the wellbore annulus are in the vertical Z direction, but they are represented in the Y direction for better visualization. The natural fracture opens in the horizontal X-Y plane. The formation is in plane strain condition in the horizontal X-Y plane.
Energies 14 02015 g004
Figure 5. Detailed algorithm of the establishment of the S-N curve.
Figure 5. Detailed algorithm of the establishment of the S-N curve.
Energies 14 02015 g005
Figure 6. S-N curve of the fractured formation with a poro elasto plastic model coupling drill string vibration and loss circulation using the base case simulation parameters.
Figure 6. S-N curve of the fractured formation with a poro elasto plastic model coupling drill string vibration and loss circulation using the base case simulation parameters.
Energies 14 02015 g006
Figure 7. Fracture width evolution in function of the time for different magnitude cyclic load and constant number of cyclic loads before failure.
Figure 7. Fracture width evolution in function of the time for different magnitude cyclic load and constant number of cyclic loads before failure.
Energies 14 02015 g007
Figure 8. Fracture width versus fracture length for different cyclic load magnitude and constant number of cyclic loads at the final stage of simulation (t = 100 s).
Figure 8. Fracture width versus fracture length for different cyclic load magnitude and constant number of cyclic loads at the final stage of simulation (t = 100 s).
Energies 14 02015 g008
Figure 9. Fracture width evolution in function of the time for constant cyclic load magnitude and different number of cyclic loads.
Figure 9. Fracture width evolution in function of the time for constant cyclic load magnitude and different number of cyclic loads.
Energies 14 02015 g009
Figure 10. Fracture width versus fracture length for constant cyclic load magnitude and different number of cyclic loads aft the final stage of simulation (t = 100 s).
Figure 10. Fracture width versus fracture length for constant cyclic load magnitude and different number of cyclic loads aft the final stage of simulation (t = 100 s).
Energies 14 02015 g010
Figure 11. Fracture length as a function of time profiles for different vibration cyclic loads.
Figure 11. Fracture length as a function of time profiles for different vibration cyclic loads.
Energies 14 02015 g011
Figure 12. Effect of mud density on fracture length and fracture width with and without vibration cyclic loads.
Figure 12. Effect of mud density on fracture length and fracture width with and without vibration cyclic loads.
Energies 14 02015 g012
Figure 13. Effect of mud density on loss circulation rate and BHP with and without vibration cyclic loads.
Figure 13. Effect of mud density on loss circulation rate and BHP with and without vibration cyclic loads.
Energies 14 02015 g013
Figure 14. Effect of mud viscosity on fracture width and fracture length with and without vibration cyclic loads.
Figure 14. Effect of mud viscosity on fracture width and fracture length with and without vibration cyclic loads.
Energies 14 02015 g014
Figure 15. Effect of mud viscosity on loss circulation rate and BHP with and without vibration cyclic loads.
Figure 15. Effect of mud viscosity on loss circulation rate and BHP with and without vibration cyclic loads.
Energies 14 02015 g015
Figure 16. Effect of pump rate on fracture length and fracture width with and without vibration cyclic loads.
Figure 16. Effect of pump rate on fracture length and fracture width with and without vibration cyclic loads.
Energies 14 02015 g016
Figure 17. Effect of pump rate on loss circulation rate and BHP with and without vibration cyclic loads.
Figure 17. Effect of pump rate on loss circulation rate and BHP with and without vibration cyclic loads.
Energies 14 02015 g017
Figure 18. Time evolution of fracture aperture near the injection point (K vertex).
Figure 18. Time evolution of fracture aperture near the injection point (K vertex).
Energies 14 02015 g018
Figure 19. Time evolution of fracture length near the injection point (K vertex).
Figure 19. Time evolution of fracture length near the injection point (K vertex).
Energies 14 02015 g019
Figure 20. Time evolution of pressure near the injection point (K vertex).
Figure 20. Time evolution of pressure near the injection point (K vertex).
Energies 14 02015 g020
Figure 21. Fracture width distribution along fracture at the final time near the injection point (K vertex).
Figure 21. Fracture width distribution along fracture at the final time near the injection point (K vertex).
Energies 14 02015 g021
Figure 22. Pressure distribution along fracture at the final time near injection point (K vertex).
Figure 22. Pressure distribution along fracture at the final time near injection point (K vertex).
Energies 14 02015 g022
Table 1. Input parameters for base case simulation.
Table 1. Input parameters for base case simulation.
ParametersValueParametersValue
Formation size20 m × 120 mGravitational acceleration10 m/s2
Formation depth1000 mRock Young’s modulus7 000 MPa
Wellbore radius10 cmRock Poisson’s ratio0.2
Drill-pipe radius5 cmRock permeability5 md
Drill collar OD8 inRock porosity0.25
Drill collar ID3 inPeak value of the nominal stress in the damage initiation criterion0.4 MPa
Initial pore pressure10 MPaFracture energy28 J/m2
Minimum horizontal stress13 MPaBK power law parameter2.284
Maximum horizontal stress15 MPaLeakoff coefficient 5 × 10 9 m / s / P a
Cohesion 10 MPaInterface stiffness80 000 MPa
Friction angle30 DegreesPumping rate 0.36   m 3 / m i n
Pore-fluid density1.0 g / c m 3 Pore fluid viscosity1 cp
Drilling-mud density1.3 g / c m 3 Fracture pore pressure10 MPa
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Kamgue Lenwoue, A.R.; Deng, J.; Feng, Y.; Li, H.; Oloruntoba, A.; Songwe Selabi, N.B.; Marembo, M.; Sun, Y. Numerical Investigation of the Influence of the Drill String Vibration Cyclic Loads on the Development of the Wellbore Natural Fracture. Energies 2021, 14, 2015. https://doi.org/10.3390/en14072015

AMA Style

Kamgue Lenwoue AR, Deng J, Feng Y, Li H, Oloruntoba A, Songwe Selabi NB, Marembo M, Sun Y. Numerical Investigation of the Influence of the Drill String Vibration Cyclic Loads on the Development of the Wellbore Natural Fracture. Energies. 2021; 14(7):2015. https://doi.org/10.3390/en14072015

Chicago/Turabian Style

Kamgue Lenwoue, Arnaud Regis, Jingen Deng, Yongcun Feng, Haitao Li, Adefarati Oloruntoba, Naomie Beolle Songwe Selabi, Micheal Marembo, and Yuanxiu Sun. 2021. "Numerical Investigation of the Influence of the Drill String Vibration Cyclic Loads on the Development of the Wellbore Natural Fracture" Energies 14, no. 7: 2015. https://doi.org/10.3390/en14072015

APA Style

Kamgue Lenwoue, A. R., Deng, J., Feng, Y., Li, H., Oloruntoba, A., Songwe Selabi, N. B., Marembo, M., & Sun, Y. (2021). Numerical Investigation of the Influence of the Drill String Vibration Cyclic Loads on the Development of the Wellbore Natural Fracture. Energies, 14(7), 2015. https://doi.org/10.3390/en14072015

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop