Numerical Investigation of the Inﬂuence of the Drill String Vibration Cyclic Loads on the Development of the Wellbore Natural Fracture

: 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 ﬁnite element software (ABAQUS) as the numerical simulator, a poro-elasto-plastic model has been established which aimed at analyzing the inﬂuence 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 proﬁles followed a sinusoidal behavior similar to the drill string vibration cyclic load proﬁles. 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. these results are obtained with a relatively small simulation time ( t = 100 s). Further stud ‐ ies need to be conducted to extrapolate these results with larger simulation time.


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 Energies 2021, 14,2015 3 of 30 software is utilized as the numerical simulator and the results obtained were compared with small asymptotic analytical solutions to validate the numerical model.

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.

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 investigation of the influence of the dri velopment of the wellbore natural fracture is a c lowing phenomena: (1) Flow of the drilling fluid inside the wellbore (2) Fracture propagation and fracturing fluid fl (3) Porous medium deformation and pore fluid

Flow of the Drilling Fluid inside the Wellbore
The flow of the drilling fluid inside the we 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 loss = c l ρV 2 2 where, ∆P is the pressure difference between P 1 and P 2 , Pa; ∆P hydro is the hydrostatic drilling fluid pressure between P 1 and P 2 , Pa; ∆P loss 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.

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.

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: where, t n , t s1 and t s2 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 , ε s1 and ε s2 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]: where, t 0 n , t 0 s1 and t 0 s2 are respectively the peak values of the normal nominal stress (t n ), the first shear (t s1 ) and the second shear (t s2 ) 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 C s = G C t ) the governing equation of the BK mixed mode behavior model is given as: where, G C n , G C s and G C t 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 st = G s + G t , G T = G st + 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.

Fracturing Fluid Flow
The Reynold's lubrication theory governs the flow of the fluid within the fracture and is given by [22]: ∂w ∂t 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]: 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.

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.

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]: 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.

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]: where, σ ij and σ 0 ij are respectively the effective stress and the initial effective stress, Pa; G is the shear modulus, Pa; ε ij and ε kk are respectively the strain tensor and the volumetric strain tensor; δ ij 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.

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: • anisms 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: • The dimensionless parameter 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 ≪ 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 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.

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 lc . C L is the fluid leakoff 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 lc ); • 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]: where, E is the plane strain modulus, Pa; µ , K and C are respectively the alternate viscosity, Pa.s; fracture toughness, Pa. √ m and leak-off coefficient, m/s/Pa.

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:

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]: where, V inject (t), V crack (t) and V leak (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: 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.

Elasticity Equation
The relationship between the fracture width w and the net pressure p is described with the elasticity equation expressed by [25]:

Lubrication Equation
The flow of incompressible Newtonian fluid inside the fracture is governed by the lubrication equation given as [27]: 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.

Linear Elastic Fracture Mechanics
The fracture propagation is expressed in the framework of linear elastic fracture mechanics theory with the following equation [27]: 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 lc ), are used to complete the previous equations:

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]: The small parameter ε, the timescale t * , the length scale L and the dimensionless parameter M k for radial crack geometry are expressed as [27]: Under the above governing equations and general scaling, the problem of hydraulically driven fracture is resumed to the following nonlinear equations: when M k 1, and assuming the fracture is always propagating (K = K lc ), the opening Ω and the pressure Π in Equation (22) can be related to the fracture length γ using fracture elastic mechanics: Finally, the global fluid balance equation governing the scaled fracture length as a function of the dimensionless scaled time is given as:

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]: For radial fracture geometry, the initial boundary conditions are given as: The other coefficients (for) and the scaled fracture length at small time are determined using the regular perturbation theory: By combining Equations (14)

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.

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.

PEER REVIEW
12 of 30 Figure 4. Schematic of the numerical model redrawn from [17]. ℎ and ℎ 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.
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. Twodimensional 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].

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: where, is the stress exerted by the drill string vibration cyclic loads on the wellbore surface at the time t, Pa; is the maximal magnitude of , Pa; is the initial phase of the vibration, rad; is the angular frequency, rad/s. is given by: 2 (47) Figure 4. Schematic of the numerical model redrawn from [17]. Sh max and Sh min 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.
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. Twodimensional 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].

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: 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: 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 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 max of the stress vibration at the impact point of the bottom hole is obtained by dividing the force due to the WOB by the crosssectional area A of the drill collar. The maximal magnitude S max 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 max .
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.

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.

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.

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  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.

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.

R REVIEW 14 of 30
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.

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 are obtained for successive decreasing values of cyclic load magnitudes . 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 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 ro is reached for a normalized ratio of the magnitude of the cyclic load to the mon tonic strength equals to 0.7 [15]. In this paper, it can be observed that the enduran limit of the rock formation is reached for a normalized ratio equals to 0.33 (the monoto strength of the rock equals 34.62 MPa). In the literature, most of the S-N curves were tablished using laboratory experiments. During the laboratory experiments, the rock sp imen is only subjected to the confining pressure and the axial deformation pressure. the numerical simulation conducted in this paper, apart from the confining pressure a the axial deformation pressure, loss circulation which occurs at the fracture wall is a taken into consideration using the pipe elements theory. This loss circulation globa weakens the rock formation and diminishes the strength of the rock formation. After duction of the strength, the magnitude of the cyclic load necessary to fail the rock is a reduced and finally the normalized ratio is therefore smaller than the ratio present in literature.

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 ti for different cyclic load magnitude (S) and constant number of cyclic loads (N) is inve gated. It is observed in Figure 7, that the fracture width on each profile increases with cyclic load magnitude. The maximal value of the fracture width without vibration cyc loads increased by 64.77% after application of cyclic loads with parameters S = 11.3 M and N = 10. Equally, the profile of the fracture width without vibration is almost smoo then, for different increasing cyclic load magnitude, the fracture width as a function time profiles follows a sinusoidal behavior similar to the drill string vibration cyclic loa profile. The drill string vibration cyclic loads which are hitting the wellbore surface ge 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 σ max to the monotonic strength σ min 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.

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 fluid pressure inside the fracture is highest near the wellbore, therefore, a la fracture width in that region is also expected. Nevertheless, stress concentration aro the wellbore is generally observed during drilling operations due to the excavation o rock from the formation. Although, the drilling mud is utilized to redistribute stre around the wellbore and ensure wellbore stability, this redistribution is not instantane Since the time simulation of this modeling is relatively small (t = 100 s), therefore s concentration is always observed around the wellbore at the final simulation stage. stress concentration in the near wellbore region acts against the opening of the frac resulting in smaller fracture width. The effect of stress concentration on the fracture w is alleviated after integration of different increasing cyclic loads because the cyclic lo buffeting the wellbore surface generate a decohesion of the near wellbore fracture lea to a significant increase of the fracture width with the cyclic loads. 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.

Development of the Wellbore natural Fracture for Constant Cyclic Load Magnitude an Different Number of Cyclic Loads
The development of the wellbore natural fracture for constant cyclic load magn and different number of cyclic loads is presented in this section. The results of the fra width evolution in function of the time are shown in Figure 9. It can be observed th fracture width as a function of time profile without vibration is almost smooth, th different increasing values of cyclic loads number, the fracture width as a function o profiles follows a sinusoidal behavior similar to the drill string vibration cyclic load file. The fracture width without vibration cyclic loads also increased after applicat the drill string vibration cyclic loads (an increase of 63.12% of the fracture width wi vibration was observed after application of cyclic loads with parameters S = 11.3 MP N = 15). However, contrary to the case of the evolution of the fracture width in funct the time for different cyclic load magnitude and constant number of cyclic loads w the maximal values of the fracture width on each profile were increasing with the load magnitude, it is observed in this case that the maximal values of the fracture on each profile are almost the same when the cyclic load magnitude is constant bu the oscillating period of the fracture width which increases due to the increase of the ber of cyclic loads.
The fracture width evolution as a function of the fracture length for different nu of cyclic loads and constant cyclic load magnitude is presented in Figure 10. The r show that the fracture width profiles are almost the same for different increasing nu of cyclic loads and constant cyclic load magnitude. The fracture width decreased wi increase of the fracture length due to the seepage of the fracturing fluid from the fra into the rock. Equally, the fracture width increased after implementation of drill vibration with constant cyclic load magnitude and different number of cyclic loads.
At the wellbore surface, the fracture width without vibration cyclic loads incr by 70.96% after application of cyclic loads with parameters S = 11.3 MPa and N = 15

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.  The fracture length development as a function of time is a major aspect which been investigated in the research. Vibration cyclic loads were integrated in the sy then, fracture length was determined at different simulation time. The fracture length lution with the time obtained after implementation of cyclic load in the model is prese in Figure 11. The results show that the fracture length with and without vibration c  The fracture length development as a function of time is a major aspect whic been investigated in the research. Vibration cyclic loads were integrated in the sy then, fracture length was determined at different simulation time. The fracture length lution with the time obtained after implementation of cyclic load in the model is prese in Figure 11. The results show that the fracture length with and without vibration c loads are slightly the same because the effect of cyclic load is localized only on the bore surface. This result demonstrates that the cyclic loads which are hitting the wel 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.
ies 2021, 14, x FOR PEER REVIEW 19 of these results are obtained with a relatively small simulation time (t = 100 s). Further stud ies need to be conducted to extrapolate these results with larger simulation time.

Discussion and Parametric Studies
In this section, parametric studies are presented to investigate the effect of other dril ing operational parameters on the development of the wellbore natural fracture in con sidering 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 r sults are presented at the final stage of simulation at t = 100 s. Discussion and parametr studies are conducted by comparing the results obtained in this paper (poro elasto plast model coupled with drill string vibration cyclic load) with the results of the poroelast model without drill string vibration cyclic loads [17].

Effect of Mud Density on Fracture Development Integrating Drill String Vibration
The mud density is an important drilling operational parameter which affects th wellbore stability. The results of the effect of mud density on fracture length and fractur width with and without vibration are shown in Figure 12. The results show that for di ferent increasing drilling mud densities, the fracture width obtained with drill string v bration are greater than those obtained without vibration due to the enlargement of th fracture width caused by the drill string vibration cyclic loads. The fracture width for e ample increased by 64.77% after integration of cyclic loads with parameters S = 11.3 MP N = 10 and mud density equal to 1.3 g/cm . Apart from mud density being 1.2 g/cm , th fracture length with and without vibration cyclic loads are slightly the same because th effect of the cyclic loads is more localized at the wellbore surface. Therefore, it is the widt of the fracture mouth which is greatly affected by the cyclic loads, the fracture length itse

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].

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. accurately predict the fracture growth development. In fact, the model without vibra tends to underestimate the fracture width and fracture length development. This un estimation can lead to several problems such as the inaccurate selection of the particle dimension (PSD) necessary to plug the natural fracture and maintain the wellbore stab during drilling operations. The effect of mud density on the loss circulation rate and the BHP with drill st vibration cyclic loads is also investigated and the results are presented in Figure 13. straightforward that an increase of mud density will lead to an increase of the hydrost mud pressure and consequently an increase of the BHP. For different increasing dril mud density, the BHP remains slightly the same with or without vibration cyclic load fact, the BHP is mainly affected by the hydrostatic mud pressure and the viscous press losses. Therefore, the cyclic loads which are applied on the wellbore natural fracture not greatly influence the value of the BHP. For mud density between 1.1 g/cm and g/cm , the loss circulation rate without vibration cyclic loads is tiny because the mo without vibration could not predict fracture initiation and propagation for this rang density.
For mud densities between 1.2 g/cm and 1.4 g/cm , the fracture initiation press has already been reached and the loss circulation rate without vibration starts increa with the mud density due to the augmentation of the BHP in the system. The model i grating vibration cyclic loads predicted fracture initiation and propagation for the wh range of mud densities between 1.1 g/cm and 1.4 g/cm , therefore the loss circula rates increase with mud density for this range of density due to the increase of the BH the system. The loss circulation rate obtained with drill string vibration are greater t the results obtained without vibration cyclic loads. An increase of 14.3% of the loss ci lation rate was observed after integration of cyclic loads with parameters S = 11.3 M and N = 10 and mud density equal to 1.3 g/cm . This result was expected because e cyclic load applied on the natural fracture progressively enlarges the fracture width. natural fracture becoming larger and larger will lead to an increased loss circulation in 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.

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 fluid viscosity lead to wider and shorter fracture [17,35].
Concerning the influence of the drill string vibration cyclic loads on the natural frac ture 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 deco hesion of the formation matrix in the neighborhood of the fracture mouth. For mud vis cosity 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 cycli loads in the system does not significantly affect the fracture length development because its impact is more localized at the wellbore surface. 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.

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]. It is also interesting to study the effect of mud viscosity on the loss circulation and the BHP with and without vibration cyclic loads. From Figure 15, it is observed the loss circulation rate and the BHP increase with the mud viscosity because the incr of the mud viscosity will result in the elevation of the viscous pressure losses which in turn lead to the augmentation of the BHP. Importantly, the increase of the BHP enhances the risk of loss circulation which is materialized in the system by an increa the loss circulation rate.  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. the loss circulation rate and the BHP increase with the mud viscosity because the increa of the mud viscosity will result in the elevation of the viscous pressure losses which w in turn lead to the augmentation of the BHP. Importantly, the increase of the BHP al enhances the risk of loss circulation which is materialized in the system by an increase the loss circulation rate. The results also showed that the BHP values do not change too much when dr string vibration cyclic loads are considered in the simulation. The maximal percentage the BHP increase after integration of drill string vibration did not exceed 0.5%. Howeve Figure 15. Effect of mud viscosity on loss circulation rate and BHP with and without vibration cyclic loads.

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. Figures 18-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 Figures 21 and 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. that the loss circulation rate and the BHP increase with the pump rate. The loss circula rate obtained with drill string vibration cyclic loads are greater than those obtained w out vibration cyclic loads. These results confirm that the tripping operations shoul made very slowly in order to avoid high values of pump rates. This will enable to pre uncontrolled fracture growth development and loss circulation.    The asymptotic solutions were obtained by following the steps of the flowchart in Fi 3. Good agreement between asymptotic solutions and numerical ABAQUS solutio observed which allows to validate the numerical model. The numerical results and clusions obtained in this study can therefore be trusted due to the establishment o verification analytical model.

Conclusions
In this research, a poro-elasto-plastic model based on a finite element method ha been established to analyze the influence of the drill string vibration cyclic loads on th development of the wellbore natural fracture. The flow of the fluid inside the wellbore i modeled with the pipe element theory in ABAQUS. Drill string vibration is modeled wit 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 zon model is used to simulate the fracture propagation and the fracturing fluid flow, while th

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:  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.