# Stick-Slip Analysis of a Drill String Subjected to Deterministic Excitation and Stochastic Excitation

Authors: Hongyuan Qiu; Jianming Yang; Stephen Butt
Journal: Shock and Vibration
 (2016)

---

## Abstract

Using a finite element model, this paper investigates the torsional vibration of a drill string under combined deterministic excitation and random excitation. The random excitation is caused by the random friction coefficients between the drill bit and the bottom of the hole and assumed as white noise. Simulation shows that the responses under random excitation become random too, and the probabilistic distribution of the responses at each discretized time instant is obtained. The two points, entering and leaving the stick stage, are examined with special attention. The results indicate that the two points become random under random excitation, and the distributions are not normal even when the excitation is assumed as Gaussian white noise.


---

## Body

## 1. Introduction

In oil and gas industry, wells are drilled for either exploration or production purposes. The drill string, a very long and slender structure, experiences various vibrations during drilling operations, which can have significant detrimental effects on the drilling system [1]. Those detrimental effects include decrease of rate of penetration (ROP), interference with measurement while drilling (MWD) tools, and causing fatigue of the drilling components. In general, three vibration modes exist in the drill string, namely, axial, torsional, and lateral. Among the three vibrations, torsional vibration has attracted significant research effort partially due to its severe negative effects on drilling efficiency and life of drill string components. In the open literature, the torsional vibration/dynamics of a drill string has been modeled in several ways [2–4]. A torsional pendulum with one or two degrees of freedom is often used [5–7] for modeling stick-slip, one special type of torsional vibration. These models are good to look at the dynamics qualitatively; however, they may not be adequate to provide more quantitative insight into the problems under investigation. Thus, some other researchers modeled the system with the Finite Element Method (FEM) [8–10]. In modeling the torsional vibration, an important factor is the excitation to the system. Some researchers accounted for the friction between the drill string and the wellbore; others focused on the resistant torque on the bit [1, 2, 5, 6, 11–13]. Existing research has revealed that the excitation from the bit-rock interaction is especially complicated [1, 6, 13], involving both friction and cutting mechanisms. Under certain conditions, the friction mechanism may cause stick-slip motion of the bit which has particularly negative effect on the drilling system. Some researchers examined the stick-slip motion and the complicated dynamics in torsion [1, 2, 5, 11–13]. Undoubtedly, these research works are helpful to understand the complex dynamics of the drill string in rotation; however, limitations exist. First, both field test and theoretical analysis have indicated that the friction mechanism between two surfaces is very complicated; the friction coefficient, in reality, is related to various different factors, such as the profile of the surface, the materials, and the lubrication conditions. The value of friction coefficient is always highly scattered. Second, the downhole condition is highly unpredictable due to the many uncertain factors in the wellbore. These factors determine that the drill string vibration and dynamics cannot be well understood with deterministic theory; rather methods from random vibration and/or stochastic dynamics would be much more powerful tools. As a matter of fact, [14] realized this point already in as early as 1950s and proposed a probabilistic model. However, research work along this direction has progressed very little since then, probably due to the conceptual complexity of random vibration. Among the very few researchers working on random vibrations of drill string, [15] investigated lateral vibration with a nonlinear random model. In the work of Chevallier, the excitation for a tricone bit and the excitation for a PDC bit were modeled as a Kanai-Tajimi process and as band limited white noise, respectively. The nonlinearity was handled with a stochastic linearization technique. In recent years, [16–19] also investigated drill string dynamics with probabilistic models. The authors’ focus was placed on the bit-rock interaction and the drilling fluid effect on dynamics. In some other papers, [20, 21] also studied the uncertainties in the weight-on-hook.In general, there is a lack of work on torsional vibration of drill strings from the random perspective. In view of this fact, this paper, which is based on the work of [22], focuses on the random torsional vibration of a drill string. FEM is used to build the dynamic model. The central difference method is then used to find the solution, and Monte Carlo (MC) simulation is carried out to obtain the statistics of the responses. The paper is organized as follows. In Section 2, a dynamic model is developed. Following that, the solution strategies used in both deterministic and random cases are presented in Section 3. Simulation results are presented and analyzed in Section 4. Finally, conclusions are drawn in Section 5.


## 2. Formulation

### 2.1. Dynamic Model

The drill string investigated in this paper is schematically shown in Figure1. For convenience of mathematical derivation, it is assumed that the top is clamped, while the ground rotates with a constant speed. This does not change the nature of problem in terms of the relative motion. The drill string, including drill pipes and drill collars, is discretized into finite elements using linear Euler-Bernoulli beam theory. If only the rotation is considered, the element stiffness matrix Ke and mass matrix Me are given by(1)
K
e
=
G
J
l
e
-
G
J
l
e
-
G
J
l
e
G
J
l
e
,
M
e
=
ρ
J
l
e
3
ρ
J
l
e
6
ρ
J
l
e
6
ρ
J
l
e
3
,where G is the shear modulus of the drill string material, ρ is the drill string density, J is the polar moment of inertia of drill string cross section, and le is the element length.Figure 1
The dynamic model of the system.By assembling the local stiffness and mass matrices, the global mass and stiffness matrices of the system can be obtained. After some mathematical manipulation, the governing equation of motion of the whole system can be represented by the following form:(2)Mq¨t+Cq˙t+Kqt=T,where q(t) denotes the global rotation displacement vector; T represents the torque excitation to the system; M, C, and K are the system global mass, damping, and stiffness matrices, respectively; C is assumed to be a linear combination of K and M as follows:(3)C=αM+βK,where α and β are constants to be selected.


### 2.2. Excitation

In (2), the excitation to the system is generalized as torque T which consists of several components explained below. The first important component of excitation is the periodic torque created by the cutting mechanism, denoted by T0. It is straightforward to assume this excitation as a periodic function of time with a frequency of Nω. For simplicity, it can be represented as a sinusoidal function as follows:(4)T0=a+bsin⁡Nωt, where a and b are constant parameters, ω is the rotary speed of the ground, and N is the number of cutters on the bit.Another important component of the excitation is torqueT1 created by the friction between the cutter and the bottom [2]. Depending on the two possible motion conditions, torque T1 can be represented as(5)T1≤Wrbμ1ω-ωb=0,=Wrbμ2ω-ωb>0, where W is the weight on bit, ωb denotes the instantaneous rotary speed of the bit, rb is a number related to the radius of the cutter, and μ1 and μ2 represent the static and kinetic friction coefficients, respectively.Given the fact that it is almost impossible to capture the exactω-ωb=0 in simulation, (5) can be recast into an alternative form as follows:(6)T1≤Wrbμ1ω-ωb≤e,=Wrbμ2ω-ωb>e, where e is a preset small number, the proper value of which can be determined by trial and error in simulation.There are several different expressions forμ1 and μ2 available in the literature [2, 5, 12]. The kinetic friction coefficient μ2 used in this paper is taken from Puebla and Alvarez-Ramirez [5] and is given as follows: (7)μ2=μcb+μsb-μcbexp⁡-γbω-ωb,where μsb and μcb are the static and Coulomb friction coefficients, respectively, γb is a positive constant, and μ2 is substantially smaller than μ1. In general. The relationship between μ1 and μ2 in (7) is graphically shown in Figure 2.Figure 2
Torque created by the friction between the cutter and the bottom.Switching between static friction and dynamic friction is the primary cause for stick-slip. Numerous papers have investigated stick-slip of the drilling bit from the deterministic point of view [2, 5, 11–13]. However, the deterministic representation of the excitation above is highly idealized compared to the reality. The friction, in particular, is hard to model with deterministic theory. Experiments [23, 24] in both drilling industry and other areas have indicated that friction is stochastic or random in nature. According to Kilburn’s work [23], friction force can be represented by combining a deterministic component and a Gaussian form random component. As such, to account for the effect of random component in the bit-rock friction, a Gaussian white noise ξ(t) with a constant spectral intensity S0 is added to the deterministic torque T1 in this paper. Excitation T, therefore, is expressed as(8)T=T0+T1+ξt.


## 2.1. Dynamic Model

The drill string investigated in this paper is schematically shown in Figure1. For convenience of mathematical derivation, it is assumed that the top is clamped, while the ground rotates with a constant speed. This does not change the nature of problem in terms of the relative motion. The drill string, including drill pipes and drill collars, is discretized into finite elements using linear Euler-Bernoulli beam theory. If only the rotation is considered, the element stiffness matrix Ke and mass matrix Me are given by(1)
K
e
=
G
J
l
e
-
G
J
l
e
-
G
J
l
e
G
J
l
e
,
M
e
=
ρ
J
l
e
3
ρ
J
l
e
6
ρ
J
l
e
6
ρ
J
l
e
3
,where G is the shear modulus of the drill string material, ρ is the drill string density, J is the polar moment of inertia of drill string cross section, and le is the element length.Figure 1
The dynamic model of the system.By assembling the local stiffness and mass matrices, the global mass and stiffness matrices of the system can be obtained. After some mathematical manipulation, the governing equation of motion of the whole system can be represented by the following form:(2)Mq¨t+Cq˙t+Kqt=T,where q(t) denotes the global rotation displacement vector; T represents the torque excitation to the system; M, C, and K are the system global mass, damping, and stiffness matrices, respectively; C is assumed to be a linear combination of K and M as follows:(3)C=αM+βK,where α and β are constants to be selected.


## 2.2. Excitation

In (2), the excitation to the system is generalized as torque T which consists of several components explained below. The first important component of excitation is the periodic torque created by the cutting mechanism, denoted by T0. It is straightforward to assume this excitation as a periodic function of time with a frequency of Nω. For simplicity, it can be represented as a sinusoidal function as follows:(4)T0=a+bsin⁡Nωt, where a and b are constant parameters, ω is the rotary speed of the ground, and N is the number of cutters on the bit.Another important component of the excitation is torqueT1 created by the friction between the cutter and the bottom [2]. Depending on the two possible motion conditions, torque T1 can be represented as(5)T1≤Wrbμ1ω-ωb=0,=Wrbμ2ω-ωb>0, where W is the weight on bit, ωb denotes the instantaneous rotary speed of the bit, rb is a number related to the radius of the cutter, and μ1 and μ2 represent the static and kinetic friction coefficients, respectively.Given the fact that it is almost impossible to capture the exactω-ωb=0 in simulation, (5) can be recast into an alternative form as follows:(6)T1≤Wrbμ1ω-ωb≤e,=Wrbμ2ω-ωb>e, where e is a preset small number, the proper value of which can be determined by trial and error in simulation.There are several different expressions forμ1 and μ2 available in the literature [2, 5, 12]. The kinetic friction coefficient μ2 used in this paper is taken from Puebla and Alvarez-Ramirez [5] and is given as follows: (7)μ2=μcb+μsb-μcbexp⁡-γbω-ωb,where μsb and μcb are the static and Coulomb friction coefficients, respectively, γb is a positive constant, and μ2 is substantially smaller than μ1. In general. The relationship between μ1 and μ2 in (7) is graphically shown in Figure 2.Figure 2
Torque created by the friction between the cutter and the bottom.Switching between static friction and dynamic friction is the primary cause for stick-slip. Numerous papers have investigated stick-slip of the drilling bit from the deterministic point of view [2, 5, 11–13]. However, the deterministic representation of the excitation above is highly idealized compared to the reality. The friction, in particular, is hard to model with deterministic theory. Experiments [23, 24] in both drilling industry and other areas have indicated that friction is stochastic or random in nature. According to Kilburn’s work [23], friction force can be represented by combining a deterministic component and a Gaussian form random component. As such, to account for the effect of random component in the bit-rock friction, a Gaussian white noise ξ(t) with a constant spectral intensity S0 is added to the deterministic torque T1 in this paper. Excitation T, therefore, is expressed as(8)T=T0+T1+ξt.


## 3. Solution Strategy

### 3.1. Solution of Dynamic Model

The response of (2) to the excitation is numerically solved by using the central difference method. In this method, the time is first discretized into short time steps Δt; then the acceleration and velocity vectors at time ti can be written as [25](9)q˙i=12Δtqi+1-qi-1,q¨i=1Δt2qi+1-2qi+qi-1,where subscript i indicates the time instant ti. For the algorithm to be stable and convergent, the time step Δt has to meet the following condition [26]:(10)Δt<Δtcrit=2ωn,where Δt and ωn represent the time step and the maximum natural frequency of the system, respectively.By substituting (9) into (2) and rearranging the terms, one has [25](11)qi+1=Δt2N1Ti+N2qi+N3qi-1,with(12)N1=M+12ΔtC-1,N2=N12M-Δt2K,N3=N112ΔtC-M.Equation (2) is highly nonlinear due to the presence of friction. For comparison purpose, it is solved for both deterministic and random cases. In the random case, the dynamic equation is solved to a certain time with randomly generated excitation samples, and the response statistics are obtained by Monte Carlo simulation.The random excitationξ(ti) in simulation is generated as follows [26]:(13)ξti=2πS0ΔtUi,where Ui are random number series generated by the computer which is normally distributed with zero mean and unit standard deviation.


### 3.2.  Strategy for Probability Density Estimate

In order to obtain the statistics of response, the ranges of the responses, bothqi and q˙i, are divided into certain numbers of small intervals. Assuming the sample size of Monte Carlo simulation is N1, the probabilistic density of the response falling into a specific small area in slip stage is computed by(14)p1=n1N1δqδq˙.In the stick case, the probabilistic density is calculated by(15)p2=n2N1δq,where δq and δq˙ are the interval length along q and q˙ directions, respectively, n1 and n2 are the numbers of the solution falling into the slip area and the stick interval, respectively, and p1 and p2 are the PDF at the specific slip area and stick interval, respectively.


## 3.1. Solution of Dynamic Model

The response of (2) to the excitation is numerically solved by using the central difference method. In this method, the time is first discretized into short time steps Δt; then the acceleration and velocity vectors at time ti can be written as [25](9)q˙i=12Δtqi+1-qi-1,q¨i=1Δt2qi+1-2qi+qi-1,where subscript i indicates the time instant ti. For the algorithm to be stable and convergent, the time step Δt has to meet the following condition [26]:(10)Δt<Δtcrit=2ωn,where Δt and ωn represent the time step and the maximum natural frequency of the system, respectively.By substituting (9) into (2) and rearranging the terms, one has [25](11)qi+1=Δt2N1Ti+N2qi+N3qi-1,with(12)N1=M+12ΔtC-1,N2=N12M-Δt2K,N3=N112ΔtC-M.Equation (2) is highly nonlinear due to the presence of friction. For comparison purpose, it is solved for both deterministic and random cases. In the random case, the dynamic equation is solved to a certain time with randomly generated excitation samples, and the response statistics are obtained by Monte Carlo simulation.The random excitationξ(ti) in simulation is generated as follows [26]:(13)ξti=2πS0ΔtUi,where Ui are random number series generated by the computer which is normally distributed with zero mean and unit standard deviation.


## 3.2.  Strategy for Probability Density Estimate

In order to obtain the statistics of response, the ranges of the responses, bothqi and q˙i, are divided into certain numbers of small intervals. Assuming the sample size of Monte Carlo simulation is N1, the probabilistic density of the response falling into a specific small area in slip stage is computed by(14)p1=n1N1δqδq˙.In the stick case, the probabilistic density is calculated by(15)p2=n2N1δq,where δq and δq˙ are the interval length along q and q˙ directions, respectively, n1 and n2 are the numbers of the solution falling into the slip area and the stick interval, respectively, and p1 and p2 are the PDF at the specific slip area and stick interval, respectively.


## 4. Simulation Results

The parameters of the drill string in the simulation are listed in Table1. In the simulation, the drill string system is divided into 30 elements, 20 elements for the drill pipes and 10 elements for the drill collars. Some other parameters used in the simulation are listed in Table 2.Table 1
Drill string parameters.
Parameters
Drill pipe
Drill collar
Length
1,000 m
200 m
Outer diameter
0.127 m
0.2286 m
Inner diameter
0.095 m
0.0762 m
Material density
7850.0 kg/m3
7850.0 kg/m3
Elastic modulus
210
×
1
0
9 N/m2
210
×
1
0
9 N/m2
Shear modulus
7.6923
×
1
0
10 N/m2
7.6923
×
1
0
10 N/m2Table 2
Other parameters used in simulation.
Parameter
Value
a
2
×
1
0
4 Nm
W
4
×
1
0
5 Nm
μ
1
0.3
S
0
15
β
0.03
N
3
μ
s
b
0.3
b
5
×
1
0
3 Nm
r
b
0.22 m
e
0.01 rad/s
α
0.03
Δ
t
0.0002
μ
c
b
0.21


### 4.1. Results from Deterministic Case

It has been well known that stick-slip happens when the rotary speed is lower than a threshold value [27]. This phenomenon is also observed in the simulation of this paper. Under the simulation parameters given above and a rotation speed of 10 rad/s (95 rpm), the response of the deterministic part is shown in Figure 3. In this figure, the time history of the relative rotation speed between the bit and the rock is given in Figure 3(a). It can be clearly seen that the stick-slip motion exists. Two representative points, the leaving point and the entering point to the stick-slip, are marked in this figure. In the deterministic case, these points (time instants) within a motion period are uniquely determined. In Figure 3(b) the phase plane is shown; the stick-slip is characterized by the flat straight line. Increasing the rotary speed or decreasing the WOB can effectively mitigate the stick-slip [11]. However, decreasing the WOB should be done cautiously because it may have negative influences on drilling efficiency. A simulation with the same parameters but a higher rotation speed of 35 rad/s (332 rpm) is also conducted; the results are shown in Figure 4. Clearly, stick-slip is gone. The response is analogous to a sinusoidal curve, and no horizontal line exists in the phase plane any more.Figure 3
Stick-slip when the ground speed is 10 rad/s: (a) time history; (b) phase plane.
(a)
(b)Figure 4
No stick-slip when the ground speed is 35 rad/s: (a) time history; (b) phase plane.
(a)
(b)


### 4.2.  Results from Random Case

Due to the random component of the friction, it is expected that the response will become random as well. Correspondingly, the phase plane will become diffused around a single line in the deterministic case. To capture the response features of the motion in the random case, MC simulation is conducted with a sample size of 200 and noise intensity ofS0=15.Although field test data have long recognized the random nature of downhole vibration, no data have been reported on the strength of the randomness in public literature. The power spectral density ofS0=15 may not be realistic in real drilling operation; however, this does not affect the significance of the study.The first simulation case in random case is at the relatively high rotation speed of 35 rad/s which is the same as the deterministic case in Figure4. Figure 5 gives the mean response of the case in which no stick-slip is observed. Compared with Figure 4, the two cases are very similar to each other. Here, for saving of space, we do not give the diffused phase plane for the random case.Figure 5
Mean response without stick-slip: (a) time history; (b) phase plane.
(a)
(b)The stick-slip case corresponding to the deterministic case in Figure3, but with a random excitation component, is also simulated. The mean response of this case is shown in Figure 6. No obvious difference exists between this figure and Figure 3. However, if we draw all the simulation samples together in the phase plane, the single phase plane curve in Figure 3 will become diffused as shown in Figure 7. For any time instant, the points will be scattered in a range with varying probability. Using (14) and (15), the probabilistic distribution density at each time point can be obtained, as is shown in Figure 8. If stick-slip exists, the two representative points, entering and leaving the stick stage, have special interests, based on which some researchers developed control algorithms in order to mitigate the stick-slip. So knowing the times of these two points is crucial for the control strategy to be successful. In the deterministic case, the time instants for these two points are deterministic. For instance, the two points in Figure 3 aret = 184.2 s andt = 186.7 s, respectively. The statistics of these two points (marginal PDF) are shown in Figure 9. Clearly, both points are scattered in a range. In addition, the distribution is distorted from normal. The reason for that is the nonlinearity caused by the friction. In theory, for a nonlinear system, excitation with normal distribution may not lead to responses with normal distribution.Figure 6
Mean response with stick-slip: (a) time history; (b) phase plane.
(a)
(b)Figure 7
Phase plane for all MC samples between two representative points (stick-slip).Figure 8
Probability density evolution in stick-slip.Figure 9
Probability density of the switching points: (a) leaving stick; (b) leaving slip.
(a)
(b)


## 4.1. Results from Deterministic Case

It has been well known that stick-slip happens when the rotary speed is lower than a threshold value [27]. This phenomenon is also observed in the simulation of this paper. Under the simulation parameters given above and a rotation speed of 10 rad/s (95 rpm), the response of the deterministic part is shown in Figure 3. In this figure, the time history of the relative rotation speed between the bit and the rock is given in Figure 3(a). It can be clearly seen that the stick-slip motion exists. Two representative points, the leaving point and the entering point to the stick-slip, are marked in this figure. In the deterministic case, these points (time instants) within a motion period are uniquely determined. In Figure 3(b) the phase plane is shown; the stick-slip is characterized by the flat straight line. Increasing the rotary speed or decreasing the WOB can effectively mitigate the stick-slip [11]. However, decreasing the WOB should be done cautiously because it may have negative influences on drilling efficiency. A simulation with the same parameters but a higher rotation speed of 35 rad/s (332 rpm) is also conducted; the results are shown in Figure 4. Clearly, stick-slip is gone. The response is analogous to a sinusoidal curve, and no horizontal line exists in the phase plane any more.Figure 3
Stick-slip when the ground speed is 10 rad/s: (a) time history; (b) phase plane.
(a)
(b)Figure 4
No stick-slip when the ground speed is 35 rad/s: (a) time history; (b) phase plane.
(a)
(b)


## 4.2.  Results from Random Case

Due to the random component of the friction, it is expected that the response will become random as well. Correspondingly, the phase plane will become diffused around a single line in the deterministic case. To capture the response features of the motion in the random case, MC simulation is conducted with a sample size of 200 and noise intensity ofS0=15.Although field test data have long recognized the random nature of downhole vibration, no data have been reported on the strength of the randomness in public literature. The power spectral density ofS0=15 may not be realistic in real drilling operation; however, this does not affect the significance of the study.The first simulation case in random case is at the relatively high rotation speed of 35 rad/s which is the same as the deterministic case in Figure4. Figure 5 gives the mean response of the case in which no stick-slip is observed. Compared with Figure 4, the two cases are very similar to each other. Here, for saving of space, we do not give the diffused phase plane for the random case.Figure 5
Mean response without stick-slip: (a) time history; (b) phase plane.
(a)
(b)The stick-slip case corresponding to the deterministic case in Figure3, but with a random excitation component, is also simulated. The mean response of this case is shown in Figure 6. No obvious difference exists between this figure and Figure 3. However, if we draw all the simulation samples together in the phase plane, the single phase plane curve in Figure 3 will become diffused as shown in Figure 7. For any time instant, the points will be scattered in a range with varying probability. Using (14) and (15), the probabilistic distribution density at each time point can be obtained, as is shown in Figure 8. If stick-slip exists, the two representative points, entering and leaving the stick stage, have special interests, based on which some researchers developed control algorithms in order to mitigate the stick-slip. So knowing the times of these two points is crucial for the control strategy to be successful. In the deterministic case, the time instants for these two points are deterministic. For instance, the two points in Figure 3 aret = 184.2 s andt = 186.7 s, respectively. The statistics of these two points (marginal PDF) are shown in Figure 9. Clearly, both points are scattered in a range. In addition, the distribution is distorted from normal. The reason for that is the nonlinearity caused by the friction. In theory, for a nonlinear system, excitation with normal distribution may not lead to responses with normal distribution.Figure 6
Mean response with stick-slip: (a) time history; (b) phase plane.
(a)
(b)Figure 7
Phase plane for all MC samples between two representative points (stick-slip).Figure 8
Probability density evolution in stick-slip.Figure 9
Probability density of the switching points: (a) leaving stick; (b) leaving slip.
(a)
(b)


## 5. Conclusions

Using a finite element model, this paper investigates the stick-slip behavior of a drill string subjected to both deterministic excitation and random excitation. Stick-slip behaviors in both deterministic and random cases are obtained through simulation and compared. The PDF of the stick-slip response in the phase plane is also obtained. The single curve of the phase plane in the deterministic case becomes diffused in a range. Correspondingly, the two points, entering and leaving the stick stage, become diffused too. Due to the nonlinearity caused by the friction, the PDF of these two points are not normal even when the excitation is assumed as Gaussian white noise.



---

*Source: 9168747-2016-05-29.xml*