#### Abstract

Stick-slip is very harmful to the service life of drillstring. The extended Hamilton principle is applied in the paper. Then, finite element method (FEM) is employed to describe the model. The drillstring-borehole impact and friction, fluid-structure interaction, bit-rock interaction, and gravity are considered in this model. The influence of axial and torsional excitation on stick-slip is analyzed. The nonlinear motion predicted by the model is consistent with the observation results in the experiments. The research shows that the fluctuation amplitude of the bit angular velocity also increases along with the increase of driving angular velocity (torsional excitation). However, both the ratio of the maximum angular velocity of the stick-slip vibration and the fluctuation of the angular velocity are continuously reduced. Meanwhile, the strength of the stick-slip vibration has a tendency to slow down. As the axial load (axial excitation) increases, the fluctuation of the maximum angular speed of the stick-slip vibration does not change significantly, but the smaller load causes a smaller stick duration.

#### 1. Introduction

Rotary drilling system for exploration and development of oil and gas reservoirs in petroleum engineering is widely applied. The drillstring is a series of drilling tools, which may include drill pipes, drill collars, stabilizers, and auxiliary tools. The axial load, torque, bit-rock interaction, slurry-drillstring interaction, and slurry damping have negative effects on the drillstring. Due to the complex environment, a statistic data will produce some unwanted vibrations, in which stick-slip is one of the most serious ones [1–5]. The drillstring often experiences severe torsional vibration and stick-slip vibration, in which the angular velocity near the bit can be increased to two or more times of the applied angular velocity at the top of the drillstring [6–8].

Belokobyl’skii and Prokopov research the friction that induces self-excited stick-slip in the early days [9]. They present a single degree of freedom system with dry friction and think that the nonlinear friction can cause stick-slip. Rappold also analyzes a self-excited oscillating system under dry friction and studies the stick-slip vibration with chaos [10]. A nonregularized dry friction model of bit-rock interaction is applied to analyze the stick-slip motion by Hugo. The mean deviation of angular velocity is obtained. Then the stability regions are identified [11]. Navarro-López and Licéaga-Castro develop a kind of control method on dynamical sliding mode to suppress stick-slip vibration in vertical well. And the mathematical model with four DOFs is presented. Furthermore, they present a new procedure to analyze stick-slip vibration based on stability characteristics and the relationship of various sliding motions [12]. Tang et al. investigate the influence of static and kinetic friction on stick-slip vibration. And a lumped-parameter model is developed [13]. Sarker et al. develop a bond graph dynamic model in horizontal well. The friction in this model is divided into two parts, frictions in deflecting section and horizontal section. Coupled behavior among bit-rock interaction, angular velocity of bit, and weight on bit (WOT) can be predicted in this model [14].

In some researches, slurry is taken into account in the models [15–17]. Ertas et al. develop a general dynamical model of drillstring taking into account the viscosity of slurry [15]. Yigit and Christoforou develop a mathematical model to obtain the dynamic behavior with the influence of drilling slurry. The results show that an additional axial excitation can effectively suppress stick-slip vibration [16]. Ytrehus et al. propose that the friction between drillstring and borehole wall is dependent on the physical property of drilling slurry [17].

Bit-rock interaction is also an important part in the horizontal well [4, 18–20]. Kovalyshen applies drag-bit model to study coupled axial-torsional vibration. The interaction between rock and bit includes cutting and friction. Axial and torsional vibrations are coupled through bit-rock interaction [21]. Kamel and Yigit propose a more realistic nonlinear model to describe cutting and contact components with axial and torsional motions [1]. A drag-bit model with two degrees of freedom is presented by Zamanian et al. in order to study the stick-slip vibration of drillstring. The research shows that bit-rock interaction has an effect on stick-slip vibration, which can be avoided by optimizing bit parameters [22]. Kovalyshen thinks that the stick-slip vibration is caused by cutting motion and wear between rock and bit.

In the present researches, most of them are aimed at double coupled vibration or vertical well. However, there are few studies on fully coupled vibration of horizontal well with fluid, friction, and bit-rock interaction. In this paper, the nonlinear model on fully coupled vibration in horizontal well under various complex boundary conditions is applied. And the effects of axial excitation and torsional excitation on stick-slip vibration are researched. Finally, the mathematical model is verified by experiments.

#### 2. Mathematical Model

Because the drillstring in horizontal well is subjected to shear deformation, Timoshenko beam theory is employed in this paper. Moreover, the drillstring is also subjected to the friction from borehole wall, drilling slurry, loading at the end, and the driving torque at the top, which make the mathematical model of drillstring more complicated. Therefore, the motion equation is deduced by an energy method-extended Hamilton principle, which is of great advantages in dynamic modeling with complex boundary conditions [23].

##### 2.1. Basic Model

The model presented in this paper includes three kinds of vibrations including axial vibration, torsional vibration, and lateral vibration. They are shown in Figure 1.

According to the characteristics and the purpose of model, it is a crucial step in modeling to make necessary and reasonable simplifications and assumptions. If all the factors are taken into account, the whole model will be too complicated. Thus, the consumption of time and resources are huge. And the improvement of accuracy is not so much because of the complexity of the model. The following assumptions are made in this paper with a certain calculation accuracy:(1)The drillstring is horizontal and the cross section is a regular circle(2)The material of drillstring is isotropic and in linear elasticity(3)The drilling slurry is ideal

##### 2.2. Motion Equation

The extended Hamilton principle is used to derive the motion equation of drillstring. The extended Hamilton principle expression iswhere is the work done by the external forces, *T* is the kinetic energy, and *U* is the strain energy.

The model is discretized by finite element method. The drillstring has a total length of *L* and is equally divided into *n* sections that has a length of *l*. There are two nodes and 12 degrees of freedom in each unit that is shown in Figure 2.

In Figure 2, , and denote the displacement of axial, torsional, lateral, and angular displacement, respectively.

Make its local coordinates be . Because the model is discretized, the element displacement can be given aswhere *t* denotes time. , , , , , and denote the shape functions of each displacement, respectively.

###### 2.2.1. Work Done by Gravity

In the process of modeling the drillstring, gravity is a conservative force that largely affects its dynamic behavior. Work done by gravity is not related to the path of the object, but only to its initial and final positions. The work done by the conservative force gravity can be written aswhere *A* is the cross-sectional area of drillstring, denotes acceleration of gravity, and denotes the drillstring density.

Finite element method and variation are used to discrete each element. The force vector can be written aswhere is the transposed matrix of the shape function .

###### 2.2.2. Kinetic Energy of Drillstring

The kinetic energy expression of drillstring can be given aswhere *I* is cross-sectional moment matrix of inertia, **r** is angular velocity vector, and **s** is velocity vector. The expression can be written aswhere and are polar moment of inertia and cross-sectional moment of inertia, respectively.

Finally, the kinetic energy of drillstring in variational form can be written aswhere , , , , , and denote the variation of , , and , respectively.

The second derivative term is the mass matrix of the model and the rest is the force vector:

###### 2.2.3. Strain Energy of Drillstring

The drillstring is a system of mass points whose positions have changed after deformation. Thus, the strain energy of the whole system has changed. According to the principle of elastic mechanics, the strain energy can be written aswhere denotes stress tensor, denotes strain tensor, and denotes drillstring volume.

###### 2.2.4. Friction and Impact

When the drillstring rotates in a horizontal well, it will bend in some cases. Thus, friction and impact may occur between the drillstring and the borehole. In this paper, elastic force is employed to describe the friction between drillstring and borehole wall [24]. The elastic forces in the *y* and *z* directions can be given aswhere is impact factor, is coefficient of determination, denotes external radius of drillstring, denotes borehole radius, and *r* is the radial displacement of drillstring. As shown in Figure 3.

The friction between the borehole and the drillstring during rotation is modeled by friction torque in *x* direction:where is the coefficient of kinetic friction, sign denotes symbolic function, denotes normal contact force, and .

The variation form of the work done by the impact can be written aswhere denotes coordinate at node *n* and denotes displacement vector.

###### 2.2.5. Bit-Rock Interaction

Bit-rock interaction is composed of the torque on bit and the weight on bit. They can be decomposed into the forces transmitted by the cutting interface (the subscript *c*) and the wear flats (the subscript *f*) [25]. Thus, they can be written as

and can be given aswhere *d* denotes cutting depth, , denotes the number of blade of bit, denotes maximum contact force, denotes the cutting force orientation description, denotes intrinsic energy, and *a* denotes the bit radius.

The interaction can be composed of four parts:

This section dictates that the minimum forces (*W* and *T*) are applied:

Thus, the variational form of the work done by cutting and wear is

###### 2.2.6. Drillstring-Mud Interaction

Figures 4 and 5 are the inner and outer slurry of drillstring d*x*. , , and are the inner and external diameter of drillstring and the diameter of borehole wall, respectively.

In Figure 3, and denote viscous hydrodynamic force and lateral nonviscous hydrodynamic force, respectively. and denote the tangential and normal component:where is the mass per unit, , and is the pressure of internal flow.

The force caused by flow outside the drillstring along *x* direction can be written aswhere , denotes pressure coefficient, and denotes viscous force caused by the external flow.

The resultant force along *x* direction is

The resultant force from internal flow in *z* direction can be written as

The resultant force caused by external flow along *z* direction iswhere , denotes drilling slurry density, denotes frictional viscous force, and and denote the flow velocity of the external and internal slurry.

FEM is employed to discretize the physical model. Then, the second derivative term of displacement is considered as mass matrix:

The first derivative term is considered as the damping matrix:

The first term is considered as the element stiffness matrix:

The rest of the term is the force vector:

###### 2.2.7. Discretization of Motion Equation

The motion equation can be written aswhere *M*, *C*, *K*, and *f* denote the mass matrix, damping matrix, stiffness matrix, and force vector, respectively. Each matrix can be written as follows:

A damping model is presented in this paper to estimate the dissipation of system. A comprehensive estimation of the energy dissipative matrix can be written as [26]where and are the additional damping coefficient.

#### 3. Numerical Results and Discussion

There may be many factors that affect stick-slip vibration, such as the angular velocity, axial force, friction, and the load of bit. The numerical simulation about the influence of angular velocity and axial force on stick-slip vibration of the drillstring is conducted in this paper.

##### 3.1. Relationship between Torsional Excitation and Stick-Slip Vibration

Because FEM is used for discretization, the six degrees of freedom at the last node of drillstring can be obtained. Then, the angular velocity is extracted to obtain the angular velocity curve in time domain, just as shown in Figure 6.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

As shown in Figure 6, the simulation speeds are 17.9 r/min, 35.7 r/min, 53.6 r/min, 74.4 r/min, 89.3 r/min, and 107.1 r/min, respectively. When the driving speed is at the low-speed range (17.9 r/min; 35.7 r/min), the regularity of stick-slip vibration is poor. When the speed is 17.9 r/min, the angular velocity of the bit fluctuates greatly in the initial stage, which causes a large dynamic load and the premature failure of drillstring and bit.

When the velocity is 35.7 r/min, the regularity of the bit is the worst and the angular velocity fluctuates violently. When the driving velocity is in the middle section (53.6 r/min; 74.4 r/min), the angular velocity of the bit fluctuates softly, which is more regular than the driving velocity is at low speed. When the driving speed is at the high-speed range (89.3 r/min, 107.1 r/min), the amplitude of the angular velocity drops off. In particular, when the driving velocity is 89.3 r/min, the fluctuation of the bit angular velocity is small. Thus, it is an ideal driving velocity.

In order to observe the relationship between the maximum angular velocity fluctuation of the bit and the angular velocity caused by stick-slip vibration more clearly, the maximum amplitude of numerical simulation under different driving velocities is obtained and compared with the different driving velocity, just as shown in Figure 7.

Driving velocity and the maximal angular velocity of the numerical simulation are shown in Figure 6. With the drive speed increases, the maximum amplitude of the angular velocity also increases. However, the proportion between the maximum angular velocity and the driving velocity is decreasing, whose ratios are 4.8, 4.0, 3.7, 2.6, 2.43, and 2.4, respectively. It can be seen that the angular velocity fluctuation of the bit caused by stick-slip vibration constantly decreases, and the intensity of stick-slip vibration tends to be weakened. The speed fluctuation caused by stick-slip vibration tends to be modest at the high-speed range. Thus, a higher driving angular velocity can suppress the strength of stick-slip vibration.

##### 3.2. Relationship between Axial Excitation and Stick-Slip Vibration

Axial excitation mainly refers to the axial force exerted on the top of the drillstring. In this section, the axial excitation will be exerted by changing the axial load on the top of the drillstring to study the influence of different axial excitation on the dynamic behavior of the drillstring.

Since there is no independent load term in the motion equation, an initial equivalent axial force is added when calculating the force vector of element. Finally, the angular velocity of the bit is obtained.

As can be seen from Figure 8, the fluctuation of the maximum angular velocity of the stick-slip vibration does not change significantly with the increase of the axial load. Therefore, it can be seen that the axial force does not affect the maximum of angular velocity fluctuation and does not have too much influence on the drillstring. However, in the global stick-slip vibration, the frequency of stick-slip vibration caused by smaller axial load is lower than that of higher axial load. Moreover, the stick state caused by the bigger axial load lasts longer than the smaller load.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

#### 4. Experiment and Results Discussion

##### 4.1. Experimental Setup

Similarity theory is applied to design an experimental apparatus in order to verify the mathematical model [27]. Schematic diagram of experiment apparatus of drillstring constrained in a horizontal pipe is shown in Figure 8. Based on the similarity theory and the functions of the test equipment, the drillstring dynamic behavior test system mainly includes five parts: driving system, pipe system, testing system, cycling system, and loading system, just as shown in Figures 9 and 10.

The similarity criteria are deduced by the dimension analysis method. The dimensions of the parameters applied in this paper are shown in Table 1. And the dimension matrix of the system is shown in Table 2.

The dimension matrix of the system can also be written as some equations just as shown in

The following equations of are obtained on the basis of matrix:

According to the similarity theory, there are two kinds of similarity ratios. One is the geometric similarity ratio. The other one is the physical similarity ratio. All physical quantities with unit of length in nonambiguity take the same geometric ratio, such as the inside and outside diameter and length of the drillstring. Thus, the length similarity ratio , where *l* and denote the length in simulation and experiment, respectively. The parameter with an angular unit in nonambiguity remains the same. In this paper, the similarity ratio of angle , where and denote the angle used in simulation and experiment, respectively.

Based on the geometric similarity ratio, the parameters used in simulation are given as follows: Outside diameter of drillstring, Inside diameter of drillstring,

The parameters used in experiment are given as follows: Outside diameter of drillstring, Inside diameter of drillstring,

In order to ensure that the test device is consistent with the real working condition of the drillstring, it is necessary to determine the main body of the experiment apparatus; the material of the drillstring in experiment and real material in simulation (real working condition) have the same or similar stress-strain curves in order to make sure of the accuracy of the results. Thus, the models of simulation and experiment are similar in physics in this paper. The density similarity ratio , where and denote the density of drillstring in simulation and experiment, respectively. The similarity ration of elasticity modulus , where *E* and are the elasticity modulus in simulation and experiment, respectively.

Based on the physical similarity ratio, the parameters used in simulation are given as follows: The density of drillstring, The elasticity modulus of drillstring, *E* = 210 GPa The parameters used in experiment are given as follows: The density of drillstring, The elasticity modulus of drillstring,

The other similarity ratios in this paper are and . and are the axial load and angular speed in the experiment, respectively. *P* and *ω* are the axial load and angular speed in the numerical simulation, respectively. In the paper, *c*_*ω* = 2.8 and *c*_*P* = 1/9130. The experimental and the simulation parameters, applied in this paper, are shown in Tables 3 and 4, which satisfied the similarity ratios.

The function of the loading system (Figure 11(a)) is to simulate the contact between the bit and the bottom hole. And a pressure sensor is installed at the bottom to measure the axial force acting on the bit in real time. The driving system is shown in Figure 11(b). The servo motor increases the torque and reduces the rotating speed through the planetary reducer and then transmits the rotating power through the torque sensor to drive the drillstring, which can simulate the rotation exerted on the drillstring in the drilling process of a horizontal well. In the experiment, the rotating speed of the servo motor can be controlled by adjusting the servo driver; different rotating speeds are applied in the horizontal drillstring test device. The influence of the rotating speed on the dynamic state of the drillstring can be detected. In addition, the equipment for applying WOB is also integrated into the system to achieve the purpose of convenient operation.

**(a)**

**(b)**

**(c)**

**(d)**

In order to detect the dynamic state of drillstring, 5 detection points are set in the experiment apparatus. They make up the testing system (Figure 11(c)). Due to the presence of water in the annular space, it is necessary to select eddy current displacement sensors that can be waterproof and realize noncontact measurement. The test points include a total of 5 (15 eddy current displacement sensors): 3 sensors (in *x*, *y*, and *z* directions, respectively) are installed orthogonally on each test point to, respectively, detect the displacement of different cross sections of the drillstring in different directions.

In the process of drilling, the existence of circulating slurry will also have a certain impact on the dynamic behavior of the drillstring. The circulatory system (Figure 11(d)) can simulate friction slurry at different flow speed and the effect of slurry on the dynamic state of drillstring. Water is used to simulate slurry in this paper. The circulatory system includes a tank, a pump, and a flowmeter. The drillstring and the plexiglass tube form the channel of slurry. The increase of flow and flow rate of slurry in the drillstring can be observed in the flowmeter through adjusting the pump. The workflow of the system is as follows: pump the slurry into the drillstring through the flow meter, slurry flows into the annular space between the drillstring and the plexiglass tube, slurry returns to the tank through ball valve, and complete the slurry circulation.

##### 4.2. Relationship between Torsional Excitation and Stick-Slip Vibration

In order to verify the model, an experiment is carried out about the influence of the driving velocity on the stick-slip vibration according to the simulation results of the mathematical model. The angular velocity of the bit is measured by a sensor at the bottom of the drillstring. The curve is shown in Figure 12.

**(a)**

**(b)**

**(c)**

Experimental angular velocity is 100 r/min, 200 r/min, and 300 r/min, respectively. It can be observed that when the driving velocity is at the low-speed range, the angular velocity of the bit fluctuates violently. However, when the driving velocity is 200 r/min, although the amplitude fluctuation is large, the ratio of the maximum fluctuation to the driving speed is still small. When the driving speed is 300 r/min, the test speed rises and fluctuates regularly. The driving speed, the fluctuation of the test speed, and the relationship between them are shown in Figure 13 in order to identify the influence of the driving velocity on stick-slip vibration.

The comparison between the experiment and the simulation is shown in Figure 13. The numerical results in the brackets are the driving velocity of the numerical simulation under the similarity theory. It can be found that, with the driving velocity increases, the maximum angular velocity in experiment and simulation will also increase. A larger driving velocity at the top of the drillstring will inevitably lead to a larger rotating speed of the bit. However, the ratio between them is decreasing. It can be considered that the angular velocity fluctuation caused by stick-slip is decreasing and closer to the driving angular velocity. And the detriment to the drillstring and bit is reduced. The numerical simulation is consistent with the experimental results. Thus, the result proves the correctness of the mathematical model.

##### 4.3. Relationship between Axial Excitation and Stick-Slip Vibration

According to the structure of the experiment apparatus, the dynamic behavior of drillstring under different axial excitations is studied through changing the axial force, just as shown in Figure 10.

Figure 14 shows the angular velocity of the bit with different loads. The angular velocity curve is almost matched when the axial load is 0.5 kgf∼2.5 kgf. And there are no significant changes on the maximum of the angular velocity (only short-term fluctuation occurs at the initial stage). The stick-slip strength of the bit cannot be clearly distinguished under the three conditions. Therefore, the total stick-slip time is applied to measure the intensity of the stick-slip. Figure 15 shows the stick-slip duration under simulation and experiment.

**(a)**

**(b)**

**(c)**

According to statistic data, the total stick durations of the bit are 11.5 s, 13.6 s, and 14.3 s under different axial load, respectively. It can be found that, with the initial axial load increase, the total stick duration shows an upward trend. It also indicates that the strength of stick-slip is increasing. The numerical simulation results are performed at different driving velocities according to the similarity theory. The trends of simulation and experiment are exactly the same. Thus, the experimental results are consistent with the numerical simulation. Experiments have proved the correctness of the numerical simulation.

#### 5. Conclusions

The dynamic behavior of the drillstring under axial and torsional excitation is researched in this paper. The mathematical model is discretized by FEM. A numerical simulation has been conducted in MATLAB and verified by the experiment. The following conclusions can be obtained:(1)With the drive velocity increases, the maximum of the numerical simulation angular velocity also increases, but the ratio of the maximum stick-slip angular velocity to the driving angular velocity decreases. The intensity of stick-slip vibration tends to be weak. And the velocity fluctuation is modest at the high speed.(2)The axial load does not affect the maximum of angular velocity fluctuation. However, in the global stick-slip vibration, the stick-slip frequency of drillstring caused by smaller axial load is lower than that of higher axial load. Moreover, the stick state caused by the larger axial load lasts longer than that at the low-speed range.(3)The design of the experimental apparatus in this paper is feasible. And the experimental results are in good agreement with the simulation results. Finally, the mathematical model is proved to be correct.

#### Data Availability

Data are not made available as the NSF supporting the project is confidential.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 11902072), Natural Science Foundation of Northeast Petroleum University (Grant no. 2019QNL-13), and China Postdoctoral Science Foundation (Grant no. 2020M670880) to Baojin Wang.