Theoretical and Numerical Analysis of Coupled Axial-Torsional Nonlinear Vibration of Drill Strings

Based on a lumped parameter model with two degrees of freedom, the periodic response of the coupled axial-torsional nonlinear vibration of drill strings is studied by HB-AFT (harmonic balance and alternating frequency/time domain) method and numerical simulations. *e amplitude-frequency characteristic curves of axial relative displacement and torsional relative angular velocity are given to reveal the mechanism of bit bounce and stick-slip motion. *e stability of periodic response is analyzed by Floquet theory, and the boundary conditions of bifurcation of periodic response are given when parameters such as nominal drilling pressure, angular velocity of turntable, and formation stiffness are varied. *e results show that the amplitude of the periodic response of the system precipitates a spontaneous jump andHopf bifurcationmay occur when the angular velocity of the turntable is varied. *e variation of parameters may lead to the complex dynamic behavior of the system, such as period-doubling motion, quasiperiodic motion, and chaos. Bit bounce and stick-slip phenomenon can be effectively suppressed by varying the angular velocity of turntable and nominal drilling pressure.


Introduction
Rotary drilling systems using drag bits are widely used for the exploration and production of oil and gas. ey often suffer from severe vibrations which can be categorized as lateral, axial, and torsional vibrations. ese vibrations can lead to fatigue problems, reduction of bit life, unexpected changes in drilling direction, and even failure of the drill string.
As we all know, the bit bounce phenomenon in the axial vibration of the drill string and the stick-slip phenomenon in the torsional vibration are very harmful to the drilling operation. e former seriously affects the drilling efficiency, while the latter seriously affects the service life of the drill string.
erefore, the coupled axial-torsional vibration of drill string has been one of the focuses in the field of drill string vibration research in recent years. Although the composition of the drill string system is very complex, the theory and numerical analysis method based on the lumped parameter model is still an effective and practical choice.
Elsayed and Raymond [1] investigated the axial-torsional coupling of the drill string with PDC bits. e measurement of chatter showed evidence of stick-slip as well as coupling between axial and torsional modes. e sidebands in the frequency spectrum near the torsional natural frequency were considered to be an index of axial torsional coupling. Yigit et al. [2,3] used a simple dynamic model to simulate the effects of different operating conditions on the interaction of stick-slip and drill bit rebound. e equations of motion of such a system were developed by using a simplified lumped parameter model with only one compliance.
is model did not account for the effect of the flow inside and outside the collars and drill pipe or complicated cutting and friction conditions at the bit/formation interface. Chi et al. [4] studied the fatigue life of the drill string under coupled axial-torsional vibration, and the coupling was attributed to the bit-rock interaction. omas Richard et al. [5] used an alternative approach based on a discrete model of the drilling system and proposed a rate-independent description of the bit-rock interaction. e model takes into consideration the axial and torsional vibration modes and their coupling through bit-rock interaction laws, which account for both frictional contact and cutting processes at the bit-rock interface. Sampaio et al. [6] utilized nonlinear finite element analysis, considering nonlinear strain displacement and geometrical stiffening to analyze the coupled axial-torsional vibration of the drill string. Assuming linear and nonlinear strain energy, the energy variation methods were used to generate the equation. e nonlinear model was more predictive of stick-slip compared to a linear model. Zamanian et al. [7] used a discrete model with two degrees of torsional freedom and one degree of axial freedom. Frictional contact and nonlinear bit-rock interaction were considered, and this interaction was said to be the cause of axial-torsional coupling. ey used Euler-forward finite difference technique to solve numerical equations and studied the effect of mud damping on stick-slip phenomena.
e effect of torsional vibrations on axial vibrations was studied by Voronov et al. [8] with a two-degree-of-freedom model. ey took into account the nonlinearity caused by the cutting action of the drill bit and studied the influence of axial and torsional dynamics on the chip formation process. Germay et al. [9] also developed a coupled axial-torsional vibration model of the drill string to investigate the effect of axial vibrations on the stick-slip phenomenon. Gulyayev et al. [10] studied the quasistatic stability of a slender drill string under axial load, torque, rotational inertia force, and internal mud flow. For the model that includes both axial and torsional drill string dynamics, which are coupled via a rate-independent bit-rock interaction law. Besselink et al. [11] used the time scale separation between fast axial dynamics and slow torsional dynamics to study the dynamic analysis of drill string model. is semianalytical approach allows for fast and accurate computation of the axial limit cycles, leading to insight into the phenomena leading to torsional vibrations. Divenyi et al. [12] defined the forces and torques according to contact or noncontact scenarios to establish a nonsmooth system for coupled axial-torsional vibration of drill string. e nonsmoothes were modeled by adopting a smoothening procedure which is advantageous in terms of mathematical description and numerical analysis. Nandakumar and Wiercigroch [13] established a fully coupled two-degree-of-freedom model which assumes a state-dependent time delay and viscous damping for both the axial and torsional motions. A detailed stability analysis of the resultant mathematical model was conducted without making any asymptotic assumptions.
Considering the state-dependent time delay and nonlinear motion caused by dry friction and contact loss, Liu et al. [14,15] established a dynamic model of drill string system with two degrees of freedom and axial-torsional coupling.
rough simulation, the time history curves of axial vibration and torsional vibration of the system were simulated, and the parameter intervals of bit bounce and stick-slip phenomena were obtained, and the semidiscrete method was used for stability analysis.
Gupta and Wahi [16] considered the drill bit bounce and stick-slip caused by the axial regeneration effect under variable cutting depth and used the parameter-concentrated drill string axial torsion model and linear cutting force model to study the overall dynamics of the system. In 2018, Gupta and Wahi [17,18] effectively simplified the twodegree-of-freedom model of drill string to a single-degreeof-freedom system, combining the regeneration effect caused by axial vibration in cutting operation and the delay term determined by torsional vibration in cutting force. e method of multiple scales was used to study the bifurcation characteristics of a drill string system with 1 : 1 internal resonance in axial and torsional modes. It was found that, under different operation parameters and system parameters, there are supercritical and subcritical bifurcations. At the same time, Gupta and Wahi [19] systematically studied the influence of various parameters on drilling stability and made an accurate stability analysis of the power system of ideal rotary drilling. ey pointed out that torsional damping is very important, and the 1 : 1 internal resonance can significantly improve the stability of the system.
Ren et al. [20] presented a nonlinear dynamics model which is described as a simplified, equivalent, flexible shell under axial rotation and qualitative analysis are developed to study the key effective factors influencing the coupled axial/ torsional vibrations of a drill string by the method of multiple scales and the Runge-Kutta-Fehlberg method. Hosseinzadeh and Bakhtiari-Nejad [21] gave a new model including the effect of axial compliance along with the complete axial rigid-body dynamics to study the coupled axial-torsional vibration of the drill string. e effect of drill string upper and lower part lengths is studied on the stability of the system. Bakhtiari-Nejad and Hosseinzadeh [22] regarded the drill string as a continuous system and studied the dynamic stability of the torsional coupling motion of the drill string using a semidiscrete method. Sarker et al. [23] developed a bond graph dynamic model of a horizontal oil well drill string to predict longitudinal motion, torsional motion, and coupling between them excited by bit-rock interaction. A lumped segment modeling approach of vertical drill string dynamics was extended to include dynamic wellbore friction in the "build" and "horizontal" sections of the drill string. Liu et al. [24] proposed an integrated model for the axial-torsion dynamics of drill string with spatially distributed inertia and general boundary conditions. Based on the cutting geometry of the bit-rock interface, a notion called a "complex-delay" is introduced. In the proposed complex-delay model, both the state-dependent delay from the torsional vibration and the multiple regenerative effects due to the loss of contact of drill bit are taken into account. Chen et al. [25] proposed a lumped mass model with multiple degrees of freedom for torsional and axial vibration of drill string to investigate the bit bounce and stick-slip behavior. Instead of force coupling on the bit-rock interaction boundary, the kinematic coupling on axial-torsional vibration was used, and the stochastic approach on the bitrock interaction was employed.
Based on a 2-DOF coupled axial-torsional model, this paper focuses on the influence of operating parameters and stiffness of rock on vibration of drill string, which is interesting and meaningful. e structure of this paper is as follows: In Section 2, the axial-torsional coupled two-degree-of-freedom nonlinear dynamic governing equation of the drill string based on the complex bit-formation interaction model is given. In Section 3, which is also the novelty work of this paper, the dynamic governing equations of the drill string system without Taylor expansion were solved by the HB-AFT (harmonic balance-alternating frequency/time domain) method. e periodic response of the system is obtained, and the stability of the periodic solution of the system is studied. Section 4 is based on the numerical results of the Runge-Kutta method to verify the results of the theoretical analysis. Finally, Section 5 is the summary and conclusion of this paper.

Dynamic Model of Axial-Torsional Coupled Vibration of Drill String
Consider the lumped system shown in Figure 1, in which the axial stiffness k 1 and torsional stiffness k 2 are used to describe the flexibility of the drill string; the axial damping c 1 and torsional damping c 2 are caused by drilling fluid; the lumped mass M is composed of the mass of the BHA (bottom-hole-assembly) and one third of the mass of the drill string; the moment of inertia J is also composed of the moment of inertia of BHA and one third of the moment of inertia of the drill string. e motion equation of the two-degree-of-freedom lumped parameter model of the axial-torsional coupling vibration of the drill string can be written as follows [16]: where T 0 is the input torque and W 0 is the nominal drilling pressure (equal to the sum of the weight of the drill string and BHA minus the cable tension and the buoyancy of the drilling fluid). W b (weight-on-bit, WOB) and T b (torque-onbit, TOB) are the resultant force and moment of rock acting on the bit, respectively. For axial vibration, the excitation comes from the lobeshaped formation surface formed on the formation by the rolling of the bit, which will lead to the change of the cutting depth. According to [26], the definition of cutting depth is as follows: where s 0 is the initial cutting depth, that is, the amplitude of the cutting depth s. Figure 2 describes three states of axial degrees of freedom of the system. Figure 2(a) is the condition in which the drill string is suspended without contact with the ground. Figure 2(b) is in BHA static equilibrium state, and axial drilling pressure W 0 is applied, but the drill string does not rotate. W 0 is called nominal drilling pressure since drilling pressure is varied due to axial vibration of drilling when drill-bit penetrates. When the drilling pressure WOB is equal to W 0 , the static displacement U 0 is related to the penetration depth of the bit, which is regarded as the origin of the coordinate axis, that is, U � 0 (as the position of the red line shown in the following figure). e downward cutting of the bit is defined as a positive movement and U is positive. Figure 2(c) shows the process of dynamic response. e forces include the instantaneous force exerted by the formation on the bit W b , the elastic force of the drill string, and the linear damping force. erefore, W b can be given as where k c is the stratum stiffness, T b is the expression of w b , r b is the bit radius, and δ c is the cutting depth.
where the expression of the friction behavior between the bit and the formation μ(ϕ • ) in the formula is where where R p (rate of penetration, ROP) is a function of nominal bit pressure w o , turntable angular velocity Ω, and bit-formation characteristics.

Theoretical Analysis Based on the HB-AFT Method
In the dynamic control equation of drill string system, U and ϕ are the axial displacement and angular displacement of bit, Shock and Vibration 3 respectively. Based on the theoretical analysis of the HB-AFT method [27,28], the displacement fluctuation of bit in axial vibration and torsional vibration is studied. By setting . In order to facilitate analysis, we   Shock and Vibration e solution of the governing equation (8) can be given as By substituting (11) into (8), the Fourier transform method is used to obtain the harmonic coefficients. e following algebraic equation about amplitude can be obtained: Equation (12) is solved by the homotopy extension method. In order to facilitate the explanation of the solution process, we set z � [z 1 , z 2 , . . . , z h ], where h � 4n + 3, , and z h is the system parameter of the change to be studied. en, (12) can be written as where H(z) � [H 1 , H 2 , . . . , H h− 1 ] T . Set z j (z j ∈ R h ) as a regular point on the solution curve of the equation (where z 0 can be obtained by the Newton-Raphson method with given system parameters). e prediction is made along the tangent direction of this point, and the next prediction point can be given as where δ j is the step size of step j , and u j is the prediction direction, which can be determined by the following equation: where 0 − is a h − 1-dimensional zero vector. en, starting from y 0 , it is corrected by Newton iterative method along the direction perpendicular to the tangent of z j point. e algorithm is as follows: Iterate from p � 1 to stop as ‖y p+1 − y p ‖ < ε, and set z j+1 � y p+1 . eoretically, the variation of each harmonic amplitude with the system parameters can be obtained by using the previously mentioned method. However, in the calculation, in order to solve H(y p ) and H • (y 0 ), it is necessary to solve the integral in (12). Because there are trigonometric functions on the right side of (10), it is very difficult to solve analytically. As the harmonic number n is large, it is impossible to solve analytically. erefore, the numerical method is used for the calculation. Divide a period (0 through 2π) into q parts. en, τ becomes an isometric sequence [0, 2π/q, 4π/q, 6π/q, 2(q − 1)π/q]. By giving the value of y p and substituting it into (12), both x 1 and x 2 become time series. en, by substituting x 1 and x 2 into (8), it becomes the sequence of (h − 1) × q. e value of H(y p ) can be obtained by discrete Fourier transform of the sequence. e form of the derivative matrix H • (y 0 ) can be written as follows: Shock and Vibration e items in the matrix can be obtained by the following expression: e previously mentioned method is used to solve the equation. eoretically, the larger the number of harmonics is, the higher the accuracy of the solution is. However, if the value is too large, the calculation efficiency will be affected. A large number of numerical calculations have been carried out on the periodic solution of the system, and it is found that the main frequency component of the periodic solution of the drill string system is less than 10, so the harmonic number of the solution is set as n � 10.
en, the stability of periodic solution is studied by using Floquet theory, and (8) is written in the form of state equation: where 10 • , x 20 , x 20 • ] is the steady periodic solution of the system, and p − is a small disturbance. By substituting it in (10) and ignoring the higher order term of p − , there are where G � DF(x 0 − ). Obviously, G is a time-varying parameter matrix, which can be solved by the Hsu method. Taking the unit matrix as the initial value, the Floquet multiplier matrix can be obtained by numerical integration in one period, and the eigenvalues is Floquet multipliers. When the Floquet multiplier passes through the unit circle from 1, the static bifurcation occurs in the system. When the Floquet multiplier passes through the unit circle from −1, the period-doubling bifurcation occurs in the system. When the Floquet multiplier passes through the unit circle from other positions, the Hopf bifurcation occurs in the system [19].

Figures 3 and 4 are amplitude-frequency curves given by the HB-AFT method.
e red line represents the stable solution, the blue dotted line represents the unstable solution, and the cyan point is the bifurcation point. e amplitude is the maximum value of the relative displacement of axial vibration and the relative angular displacement of torsional vibration obtained by the superposition of harmonics. Table 1 shows the bifurcation form of each bifurcation point in Figure 3 and the change of the Floquet multiplier. It can be seen that, with the change of the driving angular velocity of the turntable, the variation trend of the amplitude in the axial direction and the torsion direction is relatively consistent. ere are two peaks: the first peak in the formula is small, and the system has a stable periodic solution, and the second peak deviates to the left. e system has two stable periodic solutions and one unstable periodic solution, and the bistable phenomenon appears. With the change of the angular velocity of the turntable drive, two static bifurcations occur at first, and then, with the increase of the angular velocity of the turntable drive, two Hopf bifurcations occur in the system. When the value is between the two static bifurcation points, the periodic solution of the system is unstable, and the system has two stable periodic solutions when the value is between the two static bifurcation points. Figure 5 indicates the boundary curve of static bifurcation, Hopf bifurcation, and period doubling bifurcation with different parameters. e red line is the Hopf bifurcation boundary, the blue line is the static bifurcation boundary, and the green line is the period doubling bifurcation boundary.

Numerical Results
In order to verify the results of the theoretical analysis, the numerical results based on Matlab are given this section. Figure 6 is the bifurcation diagram with the turntable driving angular velocity as the bifurcation parameter. e red and blue curves represent the positive increase and reverse decrease of the turntable driving angular velocity, respectively.
When the nominal drilling pressure is 100 kN in the softer formation, the bifurcation diagram is shown in Figure 6(b). It can be seen that the harmonic periodic response of the system bifurcates and leads to chaos under some parameters when the driving angular velocity of the turntable is taken at [1, 10] (rad/s). At the same time, the peak amplitude of increasing Ω 0 is later and larger than that of decreasing Ω 0 , which means that the parameter range of chaos occurs when the angular velocity of turntable drive increases positively. When the driving angular velocity of the turntable is at [16,32] (rad/s), the drill string system makes a simple periodic motion, and the amplitude related to the decrease is greater than that related to the increase. ere is a jump phenomenon in the steady-state response of the axialtorsional coupling vibration of the drill string, and the jump occurs before the peak amplitude with the increase of the angular velocity of the turntable drive.
ese results are consistent with the results of theoretical analysis. 6 Shock and Vibration     frequency characteristic curve of torsional vibration will jump twice compared with axial vibration.
When the angular velocity of the turntable is 38 rpm and 68 rpm and the nominal bit pressure is 50 kN, 100 kN, and 150 kN respectively, the time process of the relative displacement of the axial vibration and the angular velocity of the torsional vibration is shown in Figure 8. Figure 9 is the bit weight corresponding to Figure 8.       time, it can be observed that when the drill string system makes a simple periodic motion, the vibration frequency is basically the same as the angular velocity driven by the turntable, and the vibration frequency in the axial direction and torsion direction is the same, which has nothing to do with the nominal drilling pressure. A supplementary explanation is given to the axial phase diagram of the previously mentioned chaotic vibration and simple periodic motion (the transverse and longitudinal coordinates are the relative displacement and relative velocity of the axial vibration, respectively) and the torsional phase diagram (the transverse and longitudinal coordinates are the relative angular displacement and relative angular velocity of the torsional vibration, respectively). In the softer formation, when the nominal drilling pressure is 100 kN, the axial and torsional phase diagrams of 38 rpm and 68 rpm are shown in Figure 10. e time process is made for the angular velocity of the turntable when the jump phenomenon occurs in different nominal drilling pressure, as shown in Figure 11. Figure 12 is the corresponding drilling weight diagram, and Figure 13 is the corresponding axial phase diagram and torsion phase diagram.

Discussion and Conclusion
In this paper, considering the two states of bit contact and noncontact, the axial-torsional coupled nonlinear vibration of drill string is studied based on a two-degree-of-freedom lumped parameter model, by the HB-AFT method and numerical method, respectively. In terms of the HB-AFT method, not only are the amplitude-frequency characteristic curves of axial motion and torsional motion given, but also the stability of periodic response is analyzed by using Floquet theory, and the boundaries of static bifurcation, Hopf bifurcation, and period doubling bifurcation in different parameter planes are given. e results show that, in addition to the operating parameters such as turntable angular velocity and nominal drilling pressure, the formation stiffness also has a great influence on the coupling vibration of the drill string. For the same operating parameters, the response amplitude of the system with larger formation stiffness is larger than that with smaller formation stiffness, and the greater the formation stiffness is, the easier bit bounce occurs. For the determined formation stiffness, the greater the nominal drilling pressure is, the more likely the stick-slip phenomenon occurs in the system. Meanwhile, the greater the drilling pressure is, the greater the angular velocity is for the occurrence of bit bounce. In the numerical verification, it is found that the simulation results have a good agreement with the theoretical analysis for soft and softer formation compared with hard formation. e harmful phenomena such as bit bounce and stick-slip can be effectively suppressed by changing angular velocity of turntable and nominal drilling pressure.
Data Availability e data that support the findings of this study are available from the corresponding author upon reasonable request. 14 Shock and Vibration