Numerical Investigation of the Nonlinear Drill String Dynamics Under Stick–Slip Vibration

Abstract

This paper presents a comprehensive analysis of the influence of rotary table velocity, weight-on-bit, and viscous damping on the drill string stick–slip vibration. The analysis allows for studying the qualitative and quantitative variation of the dynamic response of the drill pipes and drill collars/bit. To achieve this goal, a robust and practical finite element (FE) model of the full-scaled drill string was developed based on a velocity-weakening formulation of the nonlinear bit–rock interaction. A detailed investigation of damping parameters was carried out. The performance of the developed model was verified through comparisons with a lumped-parameter model and a field test example. Parametric studies on the stick–slip response of the entire drill string under different field operational conditions were conducted. The dynamical time series of the system response were analyzed in terms of the phase planes, response spectra, and descriptive statistics of the drill pipes and drill collars. The findings of the study revealed that for a realistic drill string geometry, the angular velocity (i.e., mean, peak-to-peak amplitude, and standard deviation) and dominant frequency of self-excited torsional stick–slip oscillations along the drill pipes and drill collars/bit are mainly governed by the rotary table velocity. Furthermore, it was shown that the contribution of higher harmonics in the torsional stick–slip response of the drill pipes is more substantial than the drill collars/bit.

1. Introduction

Drill string is an essential component of rotary drilling systems employed in onshore/offshore oil and gas extraction. A significant length of drill string is made of thin-walled drill pipes, which are loaded in tension and transfer the top rotary motion to the drill bit [1]. The lower end of the drill string is composed of various components, such as the thick-walled drill collars and stabilizers, collectively called the bottom-hole assembly (BHA) [2]. Drill collars are loaded in compression and provide the required axial force, known as weight-on-bit (WOB), to facilitate straight drilling [3]. Figure 1 shows an offshore drilling vessel and the basic components of the drill string.

Drill string is subjected to axial (i.e., bit bounce), lateral (i.e., whirl), and torsional (i.e., stick–slip) dynamic vibrations [4]. These vibrations can exist separately or coupled together in a linear or nonlinear manner [5]. Stick–slip vibration is a severe type of torsional oscillation, detrimental to the fatigue life of drill string and downhole equipment [6]. The primary cause of stick–slip vibration is the nonlinear friction at the bit–formation interface [1,5,7]. This type of vibration involves cyclic fluctuations in the bit angular velocity, ranging from zero up to more than twice the rotary table velocity [8]. The high angular velocities of the drill bit during the slip phase can induce severe axial and lateral vibrations in the BHA, leading to excessive wear and even failure of the bit [9]. Furthermore, the high levels of torque-on-bit during the stick phase, varying between 500 and 10,000 Nm, can twist the drill string for several turns [1]. The consequent large torsional deformation during the twisting can cause failure of the drill pipes or threaded connections [2].

Torsional vibration of drill string, either uncoupled or coupled with other vibration modes, has been extensively studied using (i) lumped-parameter models (e.g., [1,4,5,10,11,12,13,14,15,16,17]), (ii) continuous-parameter models (Cosserat theory of rods) (e.g., [18,19]), and (iii) finite element method (FEM) models (e.g., [9,14,20,21,22]). The influence of system parameters on the drill string vibrations has also been analyzed in a number of studies: Dareing and Livesay (1968) investigated the effect of rotary velocity, friction, and shock sub on the longitudinal vibration of the drill string [23]. The authors assumed a prescribed periodic motion at the bit in the absence of bit–rock excitation and contact forces. Liao et al. (2012) studied the influence of unbalanced mass, friction coefficient, and drive speed on the vibrations of a rotor–stator system [24]. Germay et al. (2009) analyzed the effect of weight-on-bit and rotational speed on the axial and torsional vibrations at the drill bit [22].

Although extensive literature has been dedicated to the dynamics of the drill string, the effects of field operating parameters on the torsional vibrations have not been extensively studied. Analyzing the vibrations under different operational conditions is essential for predicting the behavior of the entire system. It may also provide some valuable information about the influential external parameters in order to avoid harmful vibration modes and achieve a reliable and efficient drilling operation. Furthermore, the majority of the previous studies have analyzed the dynamics of the drill string through vibrations of the drill collars. Even though the drill collars have a significant influence on the overall dynamic response of the system, the behavior of the drill pipes as a slender structure has been less explored and is yet to be understood [19]. Therefore, the present paper aimed at comprehensive analyses of the stick–slip dynamics of the drilling assembly as an integrated system using an efficient integrated FEM model.

The objective of this paper was twofold. First, an integrated nonlinear numerical model was developed to obtain the dynamic response of the entire drill string under stick–slip vibration. A computationally efficient approach was used to model the rate-dependent bit–rock contact interface and capture the cutting and friction effects. The nonlinear Timoshenko beam element was used to model the drill string, accounting for the shear deformation and the axial and torsional stiffness of the drill pipes and drill collars. The axial and torsional vibration modes were geometrically coupled, and the nonlinear effects of large rotation were taken into account. A linear perturbation frequency analysis was performed to estimate the eigenfrequencies of the drill string. The effect of energy dissipation due to the presence of the drill mud was incorporated along the drill pipes and drill collars using the Rayleigh viscous damping. A procedure was proposed for proper quantification of the Rayleigh damping coefficients to maintain the effect of higher vibration modes. The performance of the numerical model was compared with a five-degree-of-freedom lumped-parameter model and verified against the field test example of stick–slip.

Second, comprehensive dynamic analyses of the entire drill string under different operating conditions were carried out. The effects of rotary table velocity, weight-on-bit, and damping ratio on the torsional stick–slip vibration of the drill pipes and drill collars/bit were examined quantitatively and qualitatively. Due to the vital role of viscous damping in the overall response of the system, special attention was paid to the detailed analysis of damping. The spectral analyses were carried out to determine the effects of operating parameters on the stick–slip frequency of the system. The descriptive statistics, i.e., mean, amplitude, and standard deviation of angular velocities along the drill string, were computed for varying operational conditions.

The results showed that the amplitude, mean, and standard deviation of the angular velocities of the drill pipes and drill collars/bit under stick–slip occurrence were largely related to the rotary velocity and its threshold value. Also, the dominant frequency of stick–slip vibration was dependent on the rotary velocity, i.e., decreased with decreasing rotary velocity. The results indicated that the drill pipes can be excited at frequencies higher than the second natural vibration mode of the system with noticeable amplitudes.

The findings of the work were in good qualitative agreement with the field and experimental observations under stick–slip oscillation. The outcomes revealed new features of the overall torsional dynamics of the drill string as an integrated system.

2. Developing the FEM Model

A global three-dimensional FEM model of the drill string was developed in ABAQUS/Explicit. The overall configuration and material properties were taken from a drill string structure, as given in

Table 1. Geometry and material properties of the drill string [1].

Parameter	Variable	Value	Unit
Drill pipe specification
Drill pipe length	Lp	2000	m
Drill pipe outer diameter	Dp	0.1270	m
Drill pipe inner diameter	dp	0.1086	m
Drill collar specification
Drill collar length	Lc	150	m
Drill collar outer diameter	Dc	0.2286	m
Drill collar inner diameter	dc	0.0762	m
Material specification
Steel density	ρ	7850	kg/m3
Young’s modulus	E	2 × 1011	N/m2
Shear modulus	G	7.96 × 1010	N/m2

. The slenderness ratio (diameter divided by length) of less than 0.1 implies that the three-dimensional continuum drill string can be abstracted to a one-dimensional beam [25]. Therefore, the two-node linear interpolation beam element B31 (Timoshenko beam) was selected from the ABAQUS element library to model both the drill pipes and drill collars. The B31 is a one-dimensional line element in three-dimensional space that has stiffness associated with axial deformation, bending, and torsion [26]. Each node has six degrees of freedom: three translational and three rotational. Using the beam element in modeling the drill string can significantly reduce the computational cost yet provide a good approximation of the model response. As the tool joints have a negligible effect on the axial and torsional vibrations [27], the drill pipes were modeled with a uniform cross-sectional area. The effect of geometric nonlinearities, i.e., large rotational displacements and geometric coupling between the axial and torsional vibration modes, were taken into account. The consideration of geometric nonlinearities is vital in dynamics of structures and drill string systems [20,28,29]. Figure 2 shows the overall scheme of the modeled drill string, where x, y, and z are the translational degrees of freedom, whereas φ_x, φ_y, and φ_z are the rotational degrees of freedom around x, y, and z axes, respectively.

2.1. Bit–Rock Contact Boundary Condition

Classically, for a static weight-on-bit, the bit–rock contact interface is characterized by a frictional torque decreasing with the bit angular velocity [30]. The nonlinear velocity-weakening behavior of the torque-on-bit has been observed in field and laboratory drilling experiments and is often considered the primary cause of the drill bit stick–slip oscillations [5,31,32].

ABAQUS/Explicit employs the kinematic and penalty contact models, in which the contact surfaces, i.e., the bit and the rock bodies, should be modeled through two- or three-dimensional finite elements. The inherent nonlinearity of the contact surfaces often leads to computational challenges and convergence issues. As an alternative approach to avoid the difficulties in contact analysis and efficiently employ the beam elements for minimizing the computational cost, the frictional contact at the bit–rock interface was modeled through an equivalent frictional torque-on-bit T_b as follows:

$$T_b=\left\{\begin{array}{c}M, if \left|M\right|<T_0 (\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{k} \mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}),\hfill \\ M_{sl}, otherwise (\mathrm{s}\mathrm{l}\mathrm{i}\mathrm{p} \mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}).\hfill \end{array}\right.$$

(1)

where M is the magnitude of all the frictional moments acting on the bit (drill collar’s end) about the longitudinal axis (tangent to the bit–rock contact surface), T₀ = µ_sR_bW_b is the break-away moment (maximum static friction torque), M_sl = µ(${\dot{\phi }}_b$)R_bW_b is the sliding friction moment about the longitudinal axis, R_b is the bit radius, W_b > 0 is the weight-on-bit, and $\mu \left({\dot{\phi }}_b\right)$ is the Static–Kinetic Exponential Decay friction formulation defined as follows:

$$\mu \left({\dot{\phi }}_b\right)=\left[{\mu }_c+\left({\mu }_s-{\mu }_c\right)e^{-\gamma \left|{\dot{\phi }}_b\right|}\right]$$

(2)

where ${\dot{\phi }}_b$ = Ω_b is the bit angular velocity, µ_c and µ_s are the kinetic (Coulomb) and static friction coefficients, respectively, and γ > 0 is the constant decay coefficient. It is noted that µ_s can be greater than one depending on the surfaces, while µ_c ∈ (0, 1) and µ_s > µ_c.

The drill string was uniformly discretized with 64 elements on the pipe section and 32 elements on the collar section. The number of elements was selected based on a mesh sensitivity analysis, which verified convergence of the numerical results for the entire drill string.

2.2. Eigenfrequency Extraction

The stationary behavior of the rotating drill string is equivalent to the dynamics of the same structure fixed at the surface and free at the bit [33]. Therefore, the boundary condition at the upper end (rotary table) was completely fixed, while the lower end (bit) was only free to rotate about the longitudinal axis of the string.

The first 500 undamped natural frequencies of the drill string with the specified boundary conditions were extracted using the Lanczos linear perturbation method [26]. The criterion for extracting 500 natural modes was to obtain at least ten natural torsional frequencies of the drill string. The obtained natural frequencies are given in

Table 2. Torsional natural frequencies of the drill string.

Mode Number	Natural Frequency	Natural Frequency
n	fn (Hz)	ωn (rad/s)
1	0.19	1.21
2	0.84	5.26
3	1.60	10.02
4	2.37	14.88
5	3.14	19.75
6	3.92	24.62
7	4.69	29.48
8	5.46	34.33
9	6.23	39.16
10	7.04	44.22

, where ω = 2πf.

2.3. Viscous Damping Effect

The viscous damping due to the drill mud is an important aspect of the drill string dynamics. Increasing the torsional damping can disfavor the stick–slip oscillation [1,34,35] and extend the upper boundary of the stable operating region of the drill string [36,37,38,39]. Therefore, the inclusion of viscous damping effects is essential for realistic modeling and satisfactory results.

In the present paper, the frequency-dependent Rayleigh damping was used to model the mud viscous damping in the FEM simulations. The standard form of Rayleigh damping is defined as a linear combination of the mass M and stiffness K matrices as follows [40,41]:

$$C=\alpha M+\beta K$$

(3)

where α and β are the mass and stiffness proportional coefficients, respectively. The damping ratio for mode n of the system is given by the following [41]:

$${\xi }_n={\displaystyle \frac{\alpha }{2{\omega }_n}}+{\displaystyle \frac{\beta {\omega }_n}{2}}$$

(4)

where ω_n (rad/s) is the natural frequency of mode n. The coefficients α and β can be determined from two specified modal damping ratios, i.e., ξ_i and ξ_j. Rewriting Equation (4) in matrix form for the specified modes leads to a set of two algebraic equations as follows:

$${\displaystyle \frac{1}{2}}\left[\begin{array}{cc}1/{\omega }_i& {\omega }_i\\ 1/{\omega }_j& {\omega }_j\end{array}\right]\left\{\begin{array}{c}\alpha \\ \beta \end{array}\right\}=\left\{\begin{array}{c}{\xi }_i\\ {\xi }_j\end{array}\right\}$$

(5)

Assuming ξ_i = ξ_j = ξ in Equation (5), the coefficients α and β can be computed as follows:

$$\alpha =\xi {\displaystyle \frac{2{\omega }_i{\omega }_j}{{\omega }_i+{\omega }_j}}, \beta =\xi {\displaystyle \frac{2}{{\omega }_i+{\omega }_j}}$$

(6)

For a drill string with many modes over a large range of natural frequencies, the Rayleigh method can result in damping ratios significantly different than the target modal damping. Therefore, the natural frequencies ω_i and ω_j should be chosen to ensure that the damping ratios over the intermediate modes are close to the target damping. A preliminary investigation was conducted to determine ω_i and ω_j for which the variation of the damping ratio over the intermediate modes was minimized. To achieve ξ = 0.02, for example, ω_i was set to be the first natural frequency ω₁, whereas ω_j was set to be about the second, the third, the fifth, and the tenth natural frequency of the system. Figure 3 shows the variation of Rayleigh damping over the selected frequency bands.

The results indicated that under a narrow frequency band, i.e., ω_i =1.2 rad/s = ω₁ and ω_j = 5 rad/s ≈ ω₂, the damping ratios for higher modes were immensely larger than the target damping, the so-called de-amplification effect. However, under the wider frequency bands, i.e., ω_i = 1.2 rad/s = ω₁ and ω_j = 40 rad/s ≈ ω₁₀, the damping ratios for higher modes were smaller than the target damping. Therefore, in the rest of this paper, α and β were calculated using ω_i = ω₁ and ω_j = ω₁₀ (see Tabe 2) in order to minimize the frequency-dependent damping variations and maintain the high-frequency torsional vibrations reported in the previous studies (e.g., [22,42]). This broad frequency band gave a conservative underdamped system for intermediate modes, which is more favorable than a non-conservative overdamped system, particularly in fatigue assessment of the drill string system.

3. Comparison and Verification of the Models

3.1. Comparison of the FEM Model with a Lumped-Parameter Model

A five-degree-of-freedom lumped-parameter model was developed based on the generic lumped-parameter model of Navarro-Lopez and Cortes (2007) [11] to make comparisons with the FEM model; see Figure 4. The key assumptions of the model are as follows: (i) the drill pipes behave as torsional spring without structural damping; the drill collars behave as a rigid body, (ii) the frictional contact in the pipe threaded connections and at the pipe–borehole interface is ignored, (iii) the energy dissipation effect of the drilling mud is modeled by lumped viscous damping, (iv) the motor dynamic is not considered, (v) the rotary table is massless; the rotary velocity is constant, and (vi) the weight-on-bit is constant. The corresponding set of motion equations are as follows:

$$\left\{\begin{array}{c}{k}_p\left({\phi }_r-{\phi }_1\right)-T_r=0\\ J_p{\ddot{\phi }}_1+c_p{\dot{\phi }}_1+k_p\left({\phi }_1-{\phi }_r\right)+k_p\left({\phi }_1-{\phi }_2\right)=0\\ J_p{\ddot{\phi }}_2+c_p{\dot{\phi }}_2+k_p\left({\phi }_2-{\phi }_1\right)+k_p\left({\phi }_2-{\phi }_3\right)=0\\ J_p{\ddot{\phi }}_3+c_p{\dot{\phi }}_3+k_p\left({\phi }_3-{\phi }_2\right)+k_p\left({\phi }_3-{\phi }_b\right)=0\\ J_b{\ddot{\phi }}_b+c_b{\phi }_b+k_p\left({\phi }_b-{\phi }_3\right)+T_b=0\end{array}\right.$$

(7)

where φ_r, φ₁, φ₂, φ₃, and φ_b are the angular displacements of the rotary table, drill pipe 1, drill pipe 2, drill pipe 3, and drill collars/bit, respectively, $T_r=k_p({\varphi }_r-{\varphi }_1)$ is the torque delivered by the rotary table to the drill string (T_r > 0), and T_b is a nonlinear function representing the frictional torque-on-bit. The model parameters were calculated as follows:

$$J_p=\rho I_p(L_p/4), J_c=\rho I_c{L}_c, J_b=J_c+J_p/2, k_p=G{I}_p/(L_p/4)$$

(8)

with

$$I_p=(\pi /32)(D_p^4-d_p^4), I_c=(\pi /32)(D_c^4-d_c^4)$$

(9)

where J_p, J_c, and J_b, are the mass moment of inertia of the drill pipe, drill collar, and drill bit, respectively, k_p is the torsional stiffness of the drill pipe, and I_p and I_c are the polar moment of inertia of the pipe and the collar sections, respectively. The equivalent viscous damping coefficient of the drill pipe (c_p) and drill collar/bit (c_b) were calculated as follows:

$$c_p=\widehat{\mathrm{c}}(L_p/4), c_b=\widehat{\mathrm{c}}(L_p/8)$$

(10)

where $\widehat{\mathrm{c}}$ is the viscous damping coefficient per unit length of the drill pipe.

The frictional torque-on-bit T_b was modeled as a combination of the Karnopp (1985) model [43] and the Stribeck effect (negatively sloped friction at the slip phase) through a variable structure proposed by Sepehri et al. (1996) [44], Leine et al. (1998) [45], and Navarro-Lopez and Cortes (2007) [11], as illustrated in Figure 5. Thus,

$$T_b=\left\{\begin{array}{c}{T}_{st} ,if \left|{\dot{\phi }}_b\right|<V_d and \left|T_{st}\right|\le T_0 (\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{k} \mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}),\\ T_{0}sgn\left(T_{st}\right) ,if \left|{\dot{\phi }}_b\right|<V_d and \left|T_{st}\right|>T_0 (\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} \mathrm{f}\mathrm{r}\mathrm{o}\mathrm{m} \mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{k} \mathrm{t}\mathrm{o} \mathrm{s}\mathrm{l}\mathrm{i}\mathrm{p}),\\ T_{sl}sgn\left({\dot{\varphi }}_b\right), if \left|{\dot{\phi }}_b\right|\ge V_d (\mathrm{s}\mathrm{l}\mathrm{i}\mathrm{p} \mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}).\end{array}\right.$$

(11)

with

$$T_{st}=-c_b{\dot{\phi }}_b-k_p\left({\phi }_b-{\phi }_3\right), T_0={\mu }_s{R}_b{W}_b, T_{sl}=\mu \left({\dot{\phi }}_b\right)R_b{W}_b$$

(12)

where V_d > 0 is the threshold velocity, T_st is the reaction torque during the stick phase that must overcome the break-away torque to make the bit move, T₀ is the break-away torque (maximum static friction torque), T_sl is the sliding friction torque, R_b is the bit radius, W_b > 0 is the weight-on-the bit, and $\mu \left({\dot{\phi }}_b\right)$ is the bit exponential-decay friction coefficient defined in Equation (2). The exponential decaying behavior of the torque-on-bit (Stribeck effect) is in agreement with the field and experimental observations [5,31,32].

The drill string vibration obtained from the numerical model was compared with the lumped-parameter model to have a comparative evaluation of the models’ performance. The simulations were carried out for Ω_r = 7.33 rad/s (7.33 × 60/2π = 70 rpm), which is within the common field operating range, W_b = 40 kN, and ξ = 0.03.

The model parameters are given in

Table 3. Model parameters computed for the full drill string configuration given in Table 1.

Quantity	Variable	Value	Unit
Mass moment of inertia of drill pipe	Jp	46.6438	kgm2
Mass moment of inertia of bit	Jb	335.1180	kgm2
Torsional stiffness of drill pipe	kp	1891.8971	Nm/rad
Equivalent viscous damping coefficient for drill pipe	cp	18.9	Nms/rad
Equivalent viscous damping coefficient for drill collars/bit	cb	9.45	Nms/rad
Bit–rock interface parameters [11]
Bit outer diameter	Db	0.311 (12¼)	m (inch)
Bit radius	Rb	0.1555	m
Static friction coefficient	µs	0.8	–
Kinetic (Coulomb) friction coefficient	µc	0.5	–
Decay coefficient	γ	0.9	–
Threshold velocity	Vd	1 × 10−6	–

. The equations of motion, Equations (6) and (11), were numerically integrated using the 4th–5th order Runge–Kutta method with variable time steps through the built-in function ‘ode45’ in MATLAB. The initial condition for all of the state variables (except Ω_r) was taken as zero. To achieve ξ = 0.03, the viscous damping coefficient per unit length of drill pipes ĉ = 0.0378 Ns/rad adapted from Jansen and van den Steen (1995) [1] was used.

Figure 6 shows the stick–slip response of the drill string obtained from the FEM and lumped-parameter models. There is a good agreement between the results. Also, the limit cycling behavior is in agreement with field and experimental measurements [1,46].

3.2. Field Test Verification

The FEM and lumped-parameter models developed in this paper were verified against field test data published in Kyllingstad and Nessjøen (2009) [47]. A realistic drill string of about 3200 m long consisting of 3100 m drill pipe and 100 m BHA was modeled based on the numerical and analytical procedures presented in Section 2 and Section 3.1, respectively. Figure 7 shows the result of the comparison. Although differences were observed, a fairly good match was achieved between the simulation results and the field test, which verifies the performance of the developed models under stick–slip vibration. It is noted that the difference between simulation results and the field test may be related to the bit–rock interface parameters [48]. Due to the lack of availability of drilling parameters at the bit–rock interface, these parameters were obtained from other studies (e.g., [11]). A better agreement between the simulation results and the field test can be obtained if the bit–rock-related parameters are validated based on downhole field measurements.

4. Results and Discussion

The numerical model developed in Section 2 was used for parametric studies. The values of rotary table velocity Ω_r, weight-on-bit W_b, and damping ratio ξ were varied to investigate their effects on the dynamics of the drill string and the existence of stick–slip. Only one variable was changed at each simulation to assess its individual influence. The time series of angular velocities were obtained for observation nodes along the drill string over 110 s for which the drill string exhibits steady-state stick–slip vibrations. The first 10 s of the analysis was the ramping-up of the rotary velocity from zero to the desired level, which is not displayed in the time series.

For each simulation, the peak-to-peak angular velocities of the observation nodes were calculated from t = 70 s to t = 90 s. Figure 8 demonstrates the response amplitude of the drill pipe mid-point (1150 m above the bit) and the drill bit under stick–slip oscillations. The observed torsional oscillations of the drill string under stick–slip conditions were limit cycles (bounded periodic vibrations). Therefore, the considered time duration covered at least one full cycle of the velocity oscillation for each node, and the calculated amplitudes represented the absolute amplitude of the whole simulation for that node. Furthermore, the apparent oscillations of the peak-to-peak velocities under stick–slip situation were negligible and did not significantly affect the calculated amplitudes; see Figure 8. The statistical properties of the system response, such as mean and standard deviation of the angular velocities, were then evaluated at each node throughout the time series from t = 10 s to t = 110 s. Adequate quantification of the statistical characteristics can be further used in fatigue analysis of the drill string system.

To gain more in-depth insight into the behavior of the drill string under steady-state stick–slip vibration, the frequency analyses of the numerically obtained angular velocities were carried out. The time series of angular velocities along the drill string were treated using the fast Fourier transform (FFT). The obtained spectra indicated the discrete frequencies and amplitudes of the angular velocities that contributed to the stick–slip response of the system. The time series from t = 10 s to t = 110 s were sampled at a frequency of 80 Hz (Nyquist frequency 40 Hz) to achieve high-resolution spectra. The resolution of the obtained spectra was 0.01 Hz; therefore, the results were quantitative representations of the contributing frequencies. It is worth noting that the amplitudes of zero frequency components of the response spectra correspond to twice the mean angular velocities, which are not represented in the plots. Also, it is noted that the ordinate in each spectrum plot represents the normalized amplitude; that is, the amplitude of each frequency divided by the amplitude of the dominant frequency.

4.1. Effect of Rotary Velocity Ω_r

Four different levels of Ω_r were considered: 2.09, 5.24, 7.33, and 11.52 rad/s. The W_b and ξ were fixed at 40 kN and 0.03, respectively, with α = 0.070416 1/s and β = 0.00132142 s.

Figure 9 and Figure 10 show the time series and the corresponding phase planes of angular velocity at the bit and a drill pipe located 1900 m above the bit, respectively, under varying Ω_r. For Ω_r smaller than a particular threshold value, i.e., 2.09, 5.24, and 7.33 rad/s, the drill bit undergoes periodic stick–slip vibration, whereas the drill pipe experiences regular periodic torsional vibration; see plots (a)–(c) in Figure 9 and Figure 10 This situation can be identified by the closed trajectories of the phase planes representing the ceaseless torsional oscillations. As Ω_r increases up to a certain threshold value, the duration of the stick phase at the bit is decreased, while the amplitude of the vibration at both the drill pipe and bit is increased. This agrees with field and experimental observations (e.g., [47,49]). For higher values of Ω_r, i.e., 11.52 rad/s, which exceed the threshold value, the torsional vibrations at both the bit and drill pipe gradually damp out, and the drill string performs at stable condition with a velocity close to the rotary velocity; see plot (d) in Figure 9 and Figure 10. The suppression of stick–slip oscillation by increasing the rotary velocity beyond the threshold value is in agreement with field observations (e.g., [1,34,47]).

Table 4. Mean angular velocity µ (rad/s) at different points along the drill string for W_b = 40 kN, ξ = 0.03, and varying Ω_r.

Distance Above the Bit (m)	Ωr = 2.09 rad/s	Ωr = 5.24 rad/s	Ωr = 7.33 rad/s	Ωr = 11.52 rad/s
2150	2.09	5.24	7.33	11.52
1900	2.10	5.23	7.33	11.52
1650	2.10	5.21	7.32	11.51
1400	2.11	5.20	7.32	11.51
1150	2.11	5.19	7.31	11.50
900	2.12	5.18	7.31	11.50
650	2.12	5.16	7.31	11.50
400	2.13	5.15	7.30	11.49
150	2.14	5.15	7.30	11.49
112.5	2.14	5.15	7.30	11.49
75	2.14	5.15	7.30	11.49
37.5	2.14	5.15	7.30	11.49
0	2.14	5.15	7.30	11.49

presents the calculated mean angular velocities along the drill string for varying Ω_r. Under each Ω_r, the mean angular velocities of the drill pipes and drill collars along the drill string are in the vicinity of the rotary velocity. This indicates that the mean velocity of the drill string is directly related to the rotary velocity at the surface.

Figure 11 presents the amplitude spectra of angular velocities at the drill bit and a drill pipe located 1900 m above the bit under stick–slip vibration for different Ω_r. For low levels of Ω_r, i.e., 2.09 rad/s, the dominant frequencies of both the bit and drill pipe f_D = 0.14 Hz are considerably smaller than the first torsional natural frequency of the system f₁ = 0.19 Hz, f_D/f₁ = 0.73; see Figure 11a. This indicates that stick–slip vibration can occur at frequencies lower than the first natural frequency of the drill string. This is an interesting finding that agrees with published field observations [47]. For higher Ω_r, i.e., 5.24 and 7.33 rad/s, the dominant frequencies of both the bit and drill pipe f_D = 0.18~0.19 Hz are slightly smaller than the first torsional natural frequency of the system; see Figure 11b,c. The difference between the dominant frequencies and the first natural frequency of the system is associated with the nonlinear frictional torque due to the bit–rock interaction.

According to Figure 11a–c, the amplitude decreases with frequency at the drill bit; that is, the higher harmonics have a smaller contribution to the stick–slip vibration of the drill collars/bit, especially for higher Ω_r. In contrast, the contribution of higher harmonics is significant in the drill pipes, especially for lower Ω_r. For instance, in Figure 11a, corresponding to the drill pipe at 1900 m above the bit, the peak frequencies are 0.29, 0.86, 1.58, and 2.44 Hz, which the latter three are related to the second (0.838 Hz), third (1.595 Hz), and fourth (2.368) natural torsional frequency of the system, respectively; see Table 2. From these results, it can be concluded that under a stick–slip vibration, the contribution of higher frequencies in the response of the drill pipes is more noticeable than that in the drill collars/bit. This can be related to the smaller cross-sectional area and longer length of the drill pipe section compared to the drill collar section.

4.2. Effect of Weight-on-Bit W_b

Simulations were carried out using four different values of W_b: 10, 30, 40, and 50 kN. The Ω_r and ξ were considered to be fixed at 6.3 rad/s and 0.03, respectively, with α = 0.070416 1/s and β = 0.00132142 s.

Figure 12 and Figure 13 show the time series and the corresponding phase planes of angular velocity at the bit and a drill pipe located 1650 m above the bit, respectively, for varying W_b. For levels of W_b smaller than a certain threshold value, i.e., 10 kN, no stick–slip vibration is observed, and the torsional vibrations at both the bit and drill pipe gradually disappear; see plot (a) in Figure 12 and Figure 13. As the W_b increases and exceeds the threshold value, i.e., 30, 40, and 50 kN, the drill bit exhibits periodic stick–slip vibration, resulting in periodic torsional vibration at the drill pipe; see plots (b)–(d) in Figure 12 and Figure 13. This indicates that increasing the W_b can give rise to more severe stick–slip vibration with a larger amplitude, whereas decreasing the W_b can diminish the stick–slip condition. This is in accordance with field measurement data (e.g., [5]).

The increase in W_b is associated with a rise in the break-away torque, T₀ = µ_sR_bW_b, which is accompanied by a prolonged stick phase and a larger slip amplitude. This can be described in terms of the physical phenomenon as follows: if the weight-on-bit is increased, it will produce a greater frictional torque that needs to be overcome in order to enter the slip phase. This, in turn, requires a larger twisting of the drill string to accumulate reaction torque in the bit until it reaches the frictional torque. Therefore, the duration of the stick phase under fully developed stick–slip condition is more extended for larger W_b; see Figure 12b–d.

Table 5. Mean angular velocity µ (rad/s) at different points along the drill string for Ω_r = 6.3 rad/s, ξ = 0.03, and varying W_b.

Distance Above the Bit (m)	Wb = 10 kN	Wb = 30 kN	Wb = 40 kN	Wb = 50 kN
2150	6.3000	6.3000	6.3000	6.3000
1900	6.3010	6.3024	6.2816	6.2807
1650	6.3020	6.3048	6.2631	6.2619
1400	6.3029	6.3072	6.2450	6.2441
1150	6.3038	6.3096	6.2282	6.2274
900	6.3047	6.3119	6.2129	6.2117
650	6.3055	6.3142	6.1986	6.1968
400	6.3062	6.3164	6.1853	6.1826
150	6.3069	6.3185	6.1728	6.1694
112.5	6.3069	6.3185	6.1727	6.1693
75	6.3069	6.3185	6.1726	6.1692
37.5	6.3069	6.3185	6.1726	6.1692
0	6.3069	6.3185	6.1726	6.1693

shows the calculated mean angular velocity along the drill string. For small levels of W_b with no stick–slip vibration, i.e., W_b = 10 kN, the mean angular velocity of the entire drill string is equal to the rotary velocity. As the W_b increases and exceeds the threshold value, the mean angular velocity along the drill string becomes different from the rotary velocity. The difference between the rotary velocity and the mean angular velocity along the drill string indicates a continuous twist and, thus, torsional instability in the drill string during stick–slip vibration.

Figure 14 illustrates the response spectra of the angular velocities for varying W_b. According to this figure, for different levels of W_b, the stick–slip dominant frequency f_D~0.19 Hz remains close to the first torsional natural frequency of the system 0.19188 Hz. This implies that the dominant frequency is independent of W_b. It is noted that increasing W_b has been shown to increase the frequency of the stick–slip limit cycle depending on the operating parameter [22]. Moreover, it is observed that larger W_b can excite higher frequency components in the drill pipes; see Figure 14.

4.3. Effect of Damping Ratio ξ

To conduct a detailed investigation of the effect of mud damping on stick–slip vibration of the drill string, nine values of ξ were considered, as given in

Table 6. ξ and the corresponding values for α and β.

ξ	α	β
	(1/s)	(s)
0.03	0.070416	0.00132142
0.05	0.117361	0.00220237
0.07	0.164305	0.00308332
0.10	0.234721	0.00440474
0.15	0.352082	0.00660711
0.20	0.469442	0.00880948
0.25	0.586803	0.01101185
0.27	0.633747	0.01189280
0.30	0.704164	0.01321422

. For each ξ, the corresponding α and β were calculated using Equation (6) with ω_i = 1.2056 rad/s and ω_j = 44.2 rad/s. The Ω_r and W_b were fixed at 5.24 rad/s and 40 kN, respectively.

Figure 12 and Figure 13 illustrate the time series and phase planes of the angular velocity at the bit and a drill pipe located 1400 m above the bit, respectively, for different damping ratios. For ξ smaller than a certain value, i.e., 0.05, 0.10, and 0.20, the drill string responds with periodic torsional vibrations at the drill pipe and stick–slip at the drill bit; see plots (a)–(c) in Figure 15 and Figure 16. It can be seen that increasing ξ decreases both the duration of the stick phase at the bit and the amplitude of the torsional vibration at the drill pipe and bit. For ξ higher than a particular threshold value, i.e., 0.27, the stick–slip vibration is removed, and thus, the angular velocities of the drill pipe and drill collars/bit come close to the rotary velocity; see plot (d) in Figure 15 and Figure 16. The loss of stick–slip vibration above the threshold damping is in agreement with previous observations (e.g., [1,34,47]).

Table 7. Mean angular velocity (µ) at different points along the drill string for Ω_r = 5.24 rad/s, W_b = 40 kN, and varying ξ.

Distance Above the Bit (m)	ξ = 0.05	ξ = 0.10	ξ = 0.20	ξ = 0.27
2150	5.24	5.24	5.24	5.24
1900	5.23	5.23	5.23	5.24
1650	5.22	5.22	5.21	5.23
1400	5.21	5.21	5.20	5.23
1150	5.20	5.20	5.19	5.22
900	5.19	5.19	5.18	5.22
650	5.18	5.18	5.17	5.21
400	5.18	5.17	5.16	5.21
150	5.17	5.15	5.15	5.20
112.5	5.17	5.15	5.15	5.20
75	5.17	5.15	5.15	5.20
37.5	5.17	5.15	5.15	5.20
0	5.17	5.15	5.15	5.20

shows the mean angular velocity along the drill string for different ξ. By increasing ξ, the mean angular velocity along the drill string does not change significantly and remains in the vicinity of Ω_r.

Figure 17 shows the response spectra of the angular velocities for varying ξ. According to this figure, the dominant frequency of the stick–slip vibration f_D = 0.18~0.19 Hz for varying ξ remains close to, i.e., slightly smaller than, the first natural torsional frequency of the system 0.19188 Hz. Although the influence of viscous mud damping on the dominant frequency might be negligible, it introduces a ragged behavior in the spectra.

The comparison between the stick–slip response predicted by the analytical and the numerical model when varying damping ratios showed a noticeable difference between their threshold damping values. Figure 18 demonstrates the result of the comparison. For ξ up to 0.10, both the lumped-parameter and FEM models predict a steady-state stick–slip vibration; see Figure 18a. By reaching ξ = 0.13, however, the stick–slip is eliminated in the lumped-parameter model, whereas the FEM model shows a continuous stick–slip oscillation; see Figure 18b. The difference in the threshold damping between the numerical and analytical models may be related to the Rayleigh damping model, which does not guarantee an equal value of the damping coefficient for all the vibration modes.

4.4. System Behavior

To gain a better understanding of the dynamical behavior of the entire drill string, the variations of the amplitude and standard deviation of angular velocities along the drill string with Ω_r, W_b, and ξ are shown in Figure 19 and Figure 20, respectively.

Under each stick–slip condition, the vibration amplitude A and standard deviation σ along the drill string increase rapidly from zero at the rotary table to the maximum at the top of the drill collars and remain almost constant along the collar section. This shows that the torsional stick–slip oscillation of the drill collars/bit is more crucial than the drill pipes with regard to the intensity. The almost constant-amplitude torsional vibration along the drill collars is mainly due to the high torsional stiffness of the collar section, which makes it behave as a rigid body.

The amplitude and standard deviation of the angular velocity of the drill pipes and drill collars/bit along the drill string under stick–slip oscillation increase with Ω_r and W_b, as depicted in plots (a)–(b) in Figure 19 and Figure 20, respectively, but decrease with ξ, as depicted in plot (c) in Figure 19 and Figure 20. A Ω_r or ξ higher than a threshold value or a W_b lower than a threshold value can remove the stick–slip oscillations. As the stick–slip is removed, the amplitudes and standard deviations of the angular velocity of both the drill pipes and drill collars/bit reach zero. It is noted that although stick–slip vibration can be reduced or even eliminated by increasing the rotary velocity above the threshold value, the vibration amplitudes and standard deviations along the drill string can be intensified before reaching the threshold velocity; see plot (a) in Figure 19 and Figure 20. This, in turn, can result in larger cyclic stresses in the drill pipes and drill collars, detrimental to the life of the system.

Varying the ξ does not significantly alter the stick–slip vibration; see plot (c) in Figure 19 and Figure 20. For example, increasing the ξ from 0.07 to 0.10, which is an increase of 42.86%, decreases the standard deviation of the drill collar by only 3%. Similarly, increasing the ξ from 0.03 to 0.25, which is an increase of 833.33%, decreases the standard deviation of the drill collars/bit by 17.89%. Therefore, if Ω_r and W_b are to be kept constant, and the stick–slip is treated only by altering the viscous damping, then a large damping should be considered for mitigation. It is noted that although viscous damping can reduce the amplitude of torsional vibration along the drill string or even remove the stick–slip vibration, the realistic damping value may be too small to eliminate the oscillations at field operating conditions.

The growth of the amplitude of the torsional vibration at each point along the drill string with regard to the stick–slip is significantly affected by the rotary velocity rather than the weight-on-bit and viscous damping.

5. Summary and Conclusions

In this paper, based on mathematical and mechanical theories, a detailed numerical investigation of the drill string stick–slip motion was presented. Firstly, an efficient and robust nonlinear FE model was developed to predict the vibrations of the full drill string. The stick–slip phenomenon was efficiently modeled by incorporation of the highly nonlinear bit–rock interaction and the viscous mud-damping effect along the entire drill string. An alternative approach was proposed based on the rate-dependent behavior of the torque-on-bit that has been observed in the published data. The nonlinear Timoshenko beam (shear deformation) was implemented. The nonlinear effects of large rotational displacements, the geometrically nonlinear axial–torsional coupling effects, and the axial and torsional stiffness of the drill pipes and drill collars were considered in the model. A modal analysis was conducted to estimate the eigenfrequencies of the drill string system through a linear perturbation frequency analysis. A procedure was proposed to properly select the Rayleigh damping coefficients in order to minimize the variation of the frequency-dependent damping and retain the effect of higher vibration modes of the system. A nonlinear torsional lumped-parameter model of the full drill string was developed considering the dynamics of the drill pipes, drill collars, and drill bit, accounting for the bit–rock interaction and the equivalent mud damping. The accuracy of the developed FEM and lumped-parameter models was verified against the field test data.

Secondly, the influence of rotary velocity, weight-on-bit, and viscous mud damping on the system dynamics was investigated. Particular attention was paid to the torsional behavior of both the drill pipes and drill collars/bit rather than focusing only on the BHA response. The time series of angular velocities were obtained for different points along the drill string under each operating condition. The mean, peak-to-peak amplitude, and standard deviation of the obtained velocity time series under different operating conditions were computed and compared. Spectral analyses were carried out to determine the frequency components of the torsional vibration of the drill pipes and drill collars/bit.

The results showed that the mean angular velocity of each point along the drill string was directly related to the rotational velocity at the surface but independent of the weight-on-bit and damping ratio. Variation of rotary velocity was more influential with regard to the growth of the amplitude and standard deviation of stick–slip vibration at each point along the drill string compared to the weight-on-bit and damping ratio. The stick–slip dominant frequency was dependent mainly on the angular velocity of the rotary table and was decreased with decreasing the rotary velocity. Nevertheless, the dominant frequency appeared to be independent of the weight-on-bit and damping for the considered drill string and operating conditions. The contribution of higher frequencies to the torsional stick–slip response of the drill collars/bit was negligible. However, noticeable peaks were observed at frequencies close to the higher natural vibration modes of the system in the response spectra of the drill pipes under low rotary velocities. The intensity of the higher frequency components in the response of the drill pipes increased with decreasing the rotary velocity.

In summary, the study demonstrated new characteristics of the stick–slip motion of the entire drill string and provided the significant influence of the rotary table velocity on the self-excited torsional motions of the drill string.