The dynamics of eccentric block in a fully mechanical vertical drilling tool under the effect of torsional vibration

Abstract

In present, the drilling industry have to face the reality of low oil and gas prices, fully mechanical vertical drilling tool is likely to become the best choice for vertical drilling in the future. Without electric measurement and control system, the mechanical vertical drilling tool has a mechanically stable platform with gravity sensitivity. The stable platform included an eccentric block which is not only sensitive to gravity but also severely affected by torsional vibration. In this article, we reveal the dynamics of eccentric block under the effect of torsional vibration, which is the key element of the well deviation control accuracy in mechanical vertical drilling tool. The dynamic model of eccentric block is established using Lagrangian equation and Coulomb friction, and the law of dynamic characteristics under the effect of different amplitudes and period of torsional vibration is found. The results show that, with the deviation angle and eccentric torsion moment of eccentric block is increased, the resistance to torsional vibration increases. While the frequency of torsional vibration decreases, a subordinate period gradually appears after the rotation speed value of ω is less than 1 rad/s, the torsional vibration with large amplitude and high frequency is detrimental to the stability of the eccentric block. In addition, the largest Lyapunov exponent is estimated using the time series of simulation, the dynamics of eccentric block under the effect of torsional vibration can be switched four types, and the maximum and minimum values of largest Lyapunov exponent are 3.99 × 10⁻² and −2.98 × 10⁻², which means with the state of torsional vibration changes and the response of eccentric block would be transformed to chaos.

Keywords

Mechanical vertical drilling tool torsional vibration eccentric block gravity sensitivity non-linear analysis

Introduction

Inclining prevention, while fast drilling in high and steep structures, is a challenge in oil and gas drilling engineering. Automatic vertical drilling technology is advanced in vertical well operation based on downhole vertical drilling tools, which has been successfully applied to well trajectory control in blocks with high and steep formations. Unlike traditional inclining prevention technology, such as pendulum assembly and packed hole assembly, active deviation control is realized by automatic vertical drilling technology, the well deviation could be kept at less than 1°, the drill pipe friction is reduced obviously, the efficient of drilling is significantly improved, and the drilling cost is effectively reduced.^1,2

The vertical drilling tool (VDT) is designed for the first time in KTB borehole-Germany’s superdeep drilling program.³ After three decades development, the VDT with electronic components has been widely and commercially applied, such as Power-V,⁴ SDD,⁵ VDS,⁶ the well quality, and drilling efficiency are further improved.

The operation and maintenance costs of electric-controlled VDT were high due to its high-cost and vibration-sensitive electronics. However, the drilling industry have to face the reality of low oil and gas prices in present, drilling costs needed to be reduced further.⁷ The development of fully mechanical vertical drilling tool (MVDT) was perfectly adapted to the current economic situation. MVDT has the following advantages: (a) low costs due to its fully mechanical structure, (b) high temperature resistance, (c) high suitability and reliability without electric components, and (d) high bearing capacity for larger weight on bit (WOB) and rate of penetration (ROP). Thus, MVDT will become an optimal choice for vertical drilling in the near future.

The MVDT can rotate with drill bit and drillstring during drilling process. Figure 1 shows the working principle. There are three main components in the tool for the mechanical control system. (a) Retractable piston—apply force of inclination correction through pushout to borehole wall. (b) Plate valve—direct the mud to apply pushing force to piston at the higher side of borehole. (c) Eccentric block—locate the mud control valve through holding itself on lower side of borehole due to gravity.⁷ The eccentric block is mounted on outer cylinder using bearings, by which the eccentric block could freely rotate around its axle. An upper plate valve is rigid connected with the eccentric block, and there is an opening opposite to the eccentric block. Under the upper plate valve, there is another piece of lower plate valve with three mud flow channels rotating synchronously with bottom hole assembly (BHA).⁸

Figure 1.

Schematic of working principle of MVDT.

When the inclination reaches a certain value, the eccentric block will rotate to the lower side of borehole under its eccentric torsion moment due to gravity, the opening of upper plate is rotated to the higher side of borehole, the mud flow channel on the higher side of borehole is unfolded, and then the position is pushed out. In the drilling process, the lower plate valve rotates with BHA, the three mud flow channels alternately connected, then the pistons will be pushed out on the borehole wall by turns under the pressure difference of mud between inner tool and annulus, which is approximately equal to pressure-drop of bit nozzles. This process will produce a continuous lateral force from the high side to the low side of the direction on the drill bit, thus the well trajectory will gradually restore to vertical.

During vertical drilling process, field observations based on downhole and surface vibration measurements have indicated that the drillstring and BHA exhibit severe vibration,^9,10 which is characterized by (a) sticking phases where the bit stops and (b) slipping phases where the angular velocity of the tool increases up to two times of the imposed angular velocity. This inhomogeneity rotation speed, also called stick-slip, was first founded by Cunningham¹¹ through analyzing data of downhole measurement as early as 1968, and the research field of its mechanism was initiated by Belokobyl’skii and Prokopov¹² in 1982, by whom proposed that this inhomogeneous phenomenon was caused by the non-linear friction of the negative damping generated on the rupture surface when the drill bit broke rock. Since then, a lot of researchers concluded that the mechanism of drillstring self-excited vibration is mainly caused by drill negative damping torque^13,14 and the non-linear friction between drillstring and borehole wall, and until nowadays, it still attracts great interests of researchers as a very common phenomenon due to its great dangers and inevitability.^15–17 Torsional vibration greatly influences the performance of bottom drilling tools, and the detriment of it must be taken seriously.^18–21

For MVDT, especially using the model of push-the-bit, its retractable pistons constantly push against the borehole wall and generate cycle torque and drag to the BHA, which leads the movement of MVDT to more severely torsional vibration. As a result, the torsional vibration presents obvious feature of non-linear kinetics and exerts a huge impact on inclination preventing behavior of MVDT. Without electric measurement and control system, the performance of deviation control of MVDT is directly determined by the dynamics of eccentric block, which is severely impacted by torsional vibration. Once the torsional vibration of BHA occurs, the rotation speed of the lower plate valve will be vibrating along with the BHA, then the dynamics of eccentric block must be altered simultaneously through the interaction between the upper and lower plate valve by (a) the friction coefficient switches between static and dynamic; (b) the instantaneous power of plate friction changes over time. Therefore, the behavior of eccentric block must be presented complicated dynamics, by which the inclination preventing precision of MVDT is adversely affected. The dynamics of eccentric block under the effect of torsional vibration must be further studied to discover its effect law, which is a matter of crucial and imperative significance for optimizing and obtaining advanced MVDT with high performance.

Mathematical model

Mechanical structure of stable platform with gravity sensitivity

The mechanical model of stable platform with gravity sensitivity of MVDT is shown in Figure 2. There is an upper joint and a lower plate valve on each side of the outer cylinder, respectively, connecting with the upper drill pipe and the piston chamber. In the inner space of the out cylinder, an eccentric block is mounted by bearings between the upper flange and lower flange, the eccentric block can freely rotate to the lower side of the borehole due to gravity, while the outer cylinder, upper joint, upper flange, lower flange, support ring, and lower plate valve synchronously rotate with the upper drill pipe. On the lower side of the eccentric block, an upper plate valve with its spool opening opposite to the eccentric block is assembled by key joint, thus the upper plate valve could rotate together with the eccentric block to control the piston pushout. Under the action of the positive force formed by inner and outer pressure difference of mud, the upper plate valve can keep contacting closely with the lower plate valve all the time.

Figure 2.

Schematic of mechanical structure of gravity-sensitive stable platform.

Dynamic model of eccentric block

The cross section drawing of eccentric block is designed as sector appearance, as shown in Figure 3 because of this axial symmetry shape, and the coordinate of the mass center is set as (0, Y_C), then Y_C can be calculated according to its physics principle

Y_{C} = = \frac{\int \int_{D} yd σ}{\int \int_{D} d σ} = \frac{\int \int_{D} r^{2} \sin α drd α}{\int \int_{D} rdrd α} = \frac{\int_{α_{1}}^{α_{2}} \sin α d α \int_{r_{1}}^{r_{2}} r^{2} dr}{\int_{α_{1}}^{α_{2}} d α \int_{r_{1}}^{r_{2}} rdr}

(1)

wherein $α_{1}$ and $α_{2}$ are refer to the include angle between two boundaries to X-axis; $r_{1}$ and $r_{2}$ are refer to the equivalent inner and outer radius of the eccentric block, respectively.

Figure 3.

Cross section drawing of eccentric block.

Assume that the deviation angle of borehole is $β$ , and the deflection angle between the center line of eccentric block and the lower side of borehole is $φ$ , as shown in Figure 4. The eccentric torsion moment of eccentric block $T_{E}$ can be given as following

\begin{array}{l} T_{E} = G_{e} \sin β \sin φ Y_{C} \\ = \frac{1}{3} ρ g l (r_{2}^{3} - r_{1}^{3}) (\cos α_{1} - \cos α_{2}) \sin β \sin φ \end{array}

(2)

wherein ρ and l are, respectively, referring to the density and the length of eccentric block; and g is the gravitational acceleration. When $α_{1} = 0$ and $α_{2} = π$ , the maximum value of eccentric torque $T_{Emax}$ will be obtained

T_{Emax} = \frac{2}{3} ρ gl (r_{2}^{3} - r_{1}^{3}) \sin β \sin φ

(3)

Figure 4.

Rotation schematic of eccentric block.

In order to ensure that the eccentric block can hold its position on lower side of borehole, as shown in Figure 5, the torque equation must abide by

T_{E} \geq T_{B} + T_{P}

(4)

wherein $T_{B}$ and $T_{P}$ , respectively, refer to the friction torque of bearings and plate valve.

Figure 5.

Principle of stable platform with gravity sensitivity.

Define $μ_{b}$ is the friction coefficient of bearings lubricated by drilling mud, N_L is the overall loading on bearings, the friction torque of bearings $T_{B}$ is determined by

\begin{array}{l} T_{B} = μ_{b} N_{L} \frac{d}{2} \approx \frac{1}{2} μ_{b} d G_{e} \cos β \\ = \frac{1}{4} π μ_{b} ρ g l d_{b} (r_{2}^{2} - r_{I}^{2}) \cos β \end{array}

(5)

Figure 6 shows the contact surface schematic of upper plate valve, with d, D, $α$ , $R_{1}$ , and $R_{2}$ denoting, respectively, the inner diameter, outer diameter of upper plate valve, opening angle, inner radius, and outer radius of spool opening of the valve. Then, the friction torque of plate valve $T_{P}$ can be obtained by

T_{P} = T_{P 1} + T_{P 2} + T_{P 3}

(6)

wherein $T_{P 1}$ , $T_{P 2}$ , and $T_{P 3}$ refers to three parts friction torque according to the radius $d / 2 < R < R_{1}$ , $R_{1} \leq R \leq R_{2}$ , and $R_{2} < r < D / 2$ , respectively. For convenient processing, the shaded area $S_{1}$ and $S_{2}$ in Figure 6 is ignored, and these three friction torques can be calculated as following

{\begin{matrix} T_{P 1} = \int_{0}^{R_{1}} μ_{P} P 2 π R^{2} dR = \frac{2}{3} π μ_{p} P {R_{1}}^{3} \\ T_{P 2} = \frac{α}{2 π} \frac{2}{3} π μ_{P} P ({R_{2}}^{3} - {R_{1}}^{3}) = \frac{α}{3} μ_{P} P ({R_{2}}^{3} - {R_{1}}^{3}) \\ T_{P 3} = \frac{2}{3} π μ_{P} P ({(\frac{D}{2})}^{3} - {R_{2}}^{3}) \end{matrix}

(7)

wherein P is pressure difference of mud between inner tool and annulus.

Figure 6.

Contact surface schematic of upper plate valve.

Due to the special characteristic of bearings, the friction torque $T_{B}$ can be considered as a relatively small constant value when compared with the friction torque of plate valve $T_{P}$ . Thus, for effectively researching the torsional vibration influence on the response behavior of the gravity-sensitive stable platform, the friction torque $T_{B}$ is ignored, the bearings using in this model are all ideal bearings $(μ_{b} = 0)$ , then the torque equation (4) is simplified as $T_{E} \geq T_{P}$ . To express the deflection state of eccentric block when $T_{E}$ is equal to the friction torque $T_{P}$ , here we set a critical deflection angle $φ_{t}$ , which can be calculated as

φ_{t} = \arcsin (μ_{P} P (2 π ({R_{1}}^{3} + {(\frac{D}{2})}^{3} - {R_{2}}^{3}) + α ({R_{2}}^{3} - {R_{1}}^{3})) - ρ gl (r_{2}^{3} - r_{1}^{3}) \sin β)

(8)

Mathematical model of torsional vibration

The rotary drilling structure consists essentially of a rig, drillstrings, and BHA. The BHA in this analysis refers to the lower part of the drillstring, which includes the drill collars, stabilizers, MVDT, and the drill bit. We established the simplified drillstring torsional vibration model as shown in Figure 8. Some assumptions are made as following: (a) the drillstring is homogeneous along its entire length and simply considered as a single linear torsional spring with stiffness coefficient K.²² The borehole and the drillstring are both vertical and straight, no lateral hit motion presents; (b) the boundary conditions applied by the rig to the top drive of the drillstring correspond to a constant angular velocity Ω₀; (c) the influence of drilling mud is neglected; (d) the bottom surface of borehole is composed of a large number of asperities with the shape of ellipsoidal, the height distribution of the asperities is uniform, and the contacting asperities have the same ovality ratio; (e) each of asperities follows the Hertz theory for elliptical contacts, by which the torsional vibration of BHA can be seen as a kind of compelling vibration excited by an additional harmonic vibrating motion of the borehole rock and takes the form $T = T_{0} \sin ω t$ , as shown in Figure 7(a), in which $T_{0}$ and $ω$ are, respectively, the amplitude and inherent frequency of the exciting harmonic wave exerted by borehole bottom.^23–25

Figure 7.

(a) Compelling vibration excited by the wave shape of the bottom and (b) simplified model of the drilling system.

The BHA motion as well as the forces acting on the BHA can in principle calculated by a mechanical model show in Figure 7(b), and we only consider the discrete model of the drilling system characterized by one degree of freedom $δ$ , corresponding to the angular displacement of BHA, and by three mechanical elements, define a mass moment of the inertia of BHA $J$ , an equivalent viscous damping coefficient of C, and a spring of torsional stiffness K. Then, the angular velocity of BHA Ω can be calculated

Ω = Ω_{0} + Ω_{r}

(9)

wherein $Ω_{0}$ is the angular velocity defined by hoisting and turning mechanisms with a constant value, and $Ω_{r}$ is the appendix angular velocity induced by torsional vibration. During the drilling process, the rotation of BHA induced by torsional vibration can be defined as a kind of compelling vibration excited by harmonic wave. Thus, the equation of appendix rotation motion of BHA is given as²⁶

J \frac{d^{2} δ}{d t^{2}} + C \frac{d δ}{dt} + K δ = T_{0} \sin ω t

(10)

wherein C is equivalent viscous damping coefficient of BHA, and $δ$ is angular displacement of BHA. To solve the second-order non-homogeneous linear differential equation above, defined a critical damping coefficient $C_{C}$ and obtained

ω_{p}^{2} = \frac{K}{J}, ε = \frac{C}{C_{C}} = \frac{C}{2 J ω_{p}}, ω_{q} = \sqrt{1 - ε^{2}} ω_{p}, γ = \frac{ω}{ω_{p}}

(11)

Then, the general solution of the rotation motion equation can be solved by the numerical solution method²⁷

\begin{array}{l} δ = δ_{1} + δ_{2} = A e^{- ε ω_{p}} \sin (ω_{q} t + σ_{1}) \\ + \frac{T_{0} / K}{\sqrt{{(1 - γ^{2})}^{2} + γ^{2}}} \sin (ω t + σ_{2}) \end{array}

(12)

wherein $δ_{1} = A e^{- ε ω_{p}} \sin (ω_{q} t + σ_{1})$ is one part of solution induced by damped free vibration, due to the existence of damping; this part of motion will be vanished after a period of time and can be ignored, then the particular solution $δ_{2}$ is the solution of the rotary motion. Take derivative with respect to $δ_{2}$ , the angular velocity excited by the harmonic wave $Ω_{r}$ can be inferred as follow

\begin{array}{l} Ω_{r} = \frac{d δ_{2}}{d t} = \frac{ω T_{0} / K}{\sqrt{{(1 - γ^{2})}^{2} + γ^{2}}} \cos (ω t + σ_{2}) \\ = B \cos (ω t + σ_{2}) \end{array}

(13)

Then, the angular velocity of BHA Ω can be described as

Ω = Ω_{0} + Ω_{r} = Ω_{0} + B \cos (ω t + σ_{2})

(14)

Friction model

As a friction model widely used in engineering, the classic Coulomb model is selected to describe the friction coefficient change with slid velocity in this article. Considering the difficulty in using the discontinuous friction model during simulation, in order to obtain convergence results, the curve of Coulomb model is modified to a kind of continuous curve with the characteristic of gradually changing,^28,29 as shown in Figure 8.

Figure 8.

Friction coefficients with different slid velocities.

During the mechanical model construction, the pressure difference between upper and lower plate valve $P$ can be transformed into positive pressure force $F_{p}$ on the contact surface, then the friction torque of plate valve is determined

F_{f} = μ (v) F_{p} = μ (v) P S_{P}

(15)

wherein $μ (v)$ is the friction coefficient with slid velocity, and $S_{P}$ is the area of the contact surface of the plate valve. The equation of $μ (v)$ is obtained as following, can be illustrated as Figure 8

μ (v) = {\begin{matrix} - \frac{μ_{s}}{v_{s}} v & (| v | < v_{s}) \\ \frac{μ_{d} - μ_{s}}{v_{s} - v_{d}} v + \frac{μ_{s} v_{d} - v_{s} μ_{d}}{v_{s} - v_{d}} & (v_{s} < | v | < v_{d}) \\ μ_{d} & (| v | > v_{d}) \end{matrix}

(16)

wherein $v$ is slip velocity at the friction point; $v_{s}$ is stiction transition velocity; $v_{d}$ is the friction transition velocity; $μ_{s}$ is the static friction coefficient; and $μ_{d}$ is the dynamic friction coefficient.

Parameters determination used in simulation

Parameters used in simulation

The parameters used in simulation are shown in Table 1, including design parameters according to structure design of MVDT, measured parameters according to experimental measurement, and fitting parameters fixed by field drilling data in sections “Field measurement data analysis” and “Parameters determination for torsional vibration” as follows. The drilling process and mechanical processing of MVDT are taken into consideration to determine the design parameters of the eccentric block and the plate valve, the $μ_{s}$ and $μ_{d}$ between plate valve are calculated as 0.1 to 0.05 according to friction testing with the lubricating medium is drilling mud, and the values of $v_{s}$ and $v_{d}$ are relatively set as 10 and 1000 mm/s, respectively.³⁰ However, because of the measured pressure difference $P$ is 2 MPa, the positive pressure force $F_{p}$ can be calculated and approximately set as 1800 N using the design parameters.

Table 1.

Parameters used in simulation.

No.	Design parameters		No.	Measured parameters		No.	Fitting parameters
1	r ₁ (mm)	15
2	r ₂ (mm)	75
3	l (mm)	3000	9	$μ_{s}$	0.1
4	R ₁ (mm)	11	10	$μ_{d}$	0.05
5	R ₂ (mm)	17	11	$v_{s} (mm / s)$	10	15	$Ω_{0} (rad / s)$	10
6	D (mm)	40	12	$v_{d} (mm / s)$	1000	16	B	0, 3, 6, 10
7	d (mm)	16	13	$P (MPa)$	2	17	ω	0.5, 1, 2, 3
8	α (°)	90	14	$F_{p} (N)$	1800
9	ρ (kg/mm³)	7.8 × 10⁻⁶
10	g (N/kg)	9.8

Field measurement data analysis

Field BHA data directly reflect the motion of BHA during drilling, it can be measured and stored by measuring system installed in non-magnetic collar near the drill bit, which is shown in Figure 9. In this measuring system, downhole memories and sensors are installed, in which included a three-axis magnetometers, a three-axis accelerometers, and an angular speed sensor (gyroscope), and the sampling frequency is 100 Hz. The calculation of angular speed of BHA was presented in previous work and no more detailed description is presented here.²¹ One section of downhole angular speed data of BHA are shown in Figure 10. Before data analysis, three kinds of torsional vibration state is defined as follows

{\begin{matrix} Ω_{s 1} = Ω_{0} + B \cos (ω t + σ_{2}) & | B | < 30 % | Ω_{0} | \\ Ω_{s 2} = Ω_{0} + B \cos (ω t + σ_{2}) & 30 % | Ω_{0} | < | B | < 70 % | Ω_{0} | \\ Ω_{s 3} = Ω_{0} + B \cos (ω t + σ_{2}) & | B | > 70 % | Ω_{0} | \end{matrix}

(17)

wherein $Ω_{s 1}$ , $Ω_{s 2}$ , and $Ω_{s 3}$ , respectively, refer to the torsional vibration state of mild, moderate, and severe.

Figure 9.

Field drilling data measuring system installed near the bit.

Figure 10.

Downhole angular speed of field drilling data.

The measurement results show that, at the beginning, the torsional vibration of BHA is switched frequently between the states of mild, moderate, and severe, but gradually the amplitude of drilling torsional vibration generally increased due to the penetration deepens into difficult formation. Approximately at 1200 and 1500 s, the angular speed $Ω$ drops to zero consecutively twice, which means that the drilling system stops rotating at the connection time of three-joint unit. After restoring motivation by top drive, the torsional vibration of BHA is slight at an early stage, which can last 150–200 s. Then, the torsional vibration increased rapidly from a wild stage to a relative steady severe stage lasted for a long time, the angular speed $Ω$ even appears negative value occasionally, which signifies that the rotary-jump happened, and the BHA is bearing severe stick-slip.

Parameters determination for torsional vibration

As shown in Figure 10, the region that emphasized by red rectangle listed from No. 1, No. 2 to No. 3, respectively, refers to three kinds of typical stick-slip, namely, mild torsional vibration, moderate torsional vibration, and severe torsional vibration. For excluding the interference of random factor, these field measured data need to deal with curve fitting using the torsional vibration model of BHA shown in equation (14), the fitting results with 95% confidence bounds of these data are as follows

{\begin{matrix} Ω_{B 1} = 10.4720 - 3.0929 \cos (1.7279 t + 4.4465) \\ Ω_{B 2} = 10.4720 - 5.7611 \cos (1.7279 t + 2.7845) \\ Ω_{B 3} = 10.4720 - 11.1942 \cos (1.7279 t + 2.3721) \end{matrix}

(18)

wherein $Ω_{B 1}$ , $Ω_{B 2}$ , and $Ω_{B 3}$ , respectively, denoting the angular speed of mild, moderate, and severe torsional vibration, and the curve fitting results compared with the actual drilling data are shown in Figure 11. The goodness of fit of three typical stick-slips are shown in Table 2, including the sum of squares due to error (SSE), coefficient of determination (R²), degree-of-freedom-adjusted coefficient of determination (Adj R²), and root mean squared error (RMSE).

Figure 11.

Fixing curves of three typical section parameters corresponding to Figure 9: (a) mild torsional vibration refers to No. 1 region, (b) moderate torsional vibration refers to No. 2 region, and (c) severe torsional vibration refers to No. 3 region.

Table 2.

Goodness of fit compared to measured data.

No. of curve	State of torsional vibration	SSE	R ²	Adj R²	RMSE
1	Mild	6.7364e+03	0.7411	0.7410	1.2980
2	Moderate	9.7569e+03	0.8823	0.8822	1.5622
3	Severe	1.2951e+04	0.9505	0.9505	1.7998

SSE: sum of squares due to error; RMSE: root mean squared error.

To facilitate an effective analysis of the influence of torsional vibration, the angular speed of BHA for simulation could approximate to

{\begin{matrix} Ω_{B 1} \approx 10 - 3 \cos (ω t) \\ Ω_{B 2} \approx 10 - 6 \cos (ω t) \\ Ω_{B 3} \approx 10 - 10 \cos (ω t) \end{matrix}

(19)

where $ω$ is not a constant value during the drilling process; here, we set the possible values of $ω$ as 0.5, 1, 2, and 3 corresponding to different formation parameters and drilling process, which includes parameters of WOB, ROP, revolutions per minute (RPM), and so on.

Methods

Lagrangian approach

The dynamic equation of the gravity sensing platform is established by the combined use of the conservation of momentum and the Lagrangian approach, by which the statics, kinematics, and dynamics of the system can be analyzed. The Lagrange equation can be acquired by expressing a general dynamic equation of a non-free particle system using generalized coordinates, and on this basis, a motion differential equation of a general mechanical system can be obtained by Lagrangian multiplier method. We set the generalized coordinates of a rigid body $j_{r}$ using rectangular coordinates and Euler angles as $q_{j_{r}} = [x, y, z, Ψ, θ_{L}, ϕ]_{j_{r}}^{E}$ , then a system with n rigid bodies can be presented by $q = [q_{1}^{E}, q_{2}^{E}, \dots, q_{n}^{E}]^{E}$ , and the motion differential equation of mechanical system is

\frac{d}{dt} {(\frac{\partial E}{\partial q})}^{E} - {(\frac{\partial E}{\partial q})}^{E} + ρ_{L} ϕ_{q}^{E} + μ_{L} θ_{q}^{E} = Q

(20)

wherein E is kinetic energy of system, q is an array of generalized coordinates, Q is an array of generalized force, $ρ_{L}$ is an array of the Lagrange multiplier with full constraints, and $μ_{L}$ is an array of the Lagrange multiplier with incomplete constraints. We define $ϕ (q, t) = 0$ as the full constraint equation and $θ (q, \overset{\cdot}{q}, t) = 0$ as the non-holonomic constraint equation, and equation (20) can be rewritten by setting $u - \overset{\cdot}{q} = 0$ in general form

{\begin{matrix} F (q, u, \overset{\cdot}{u}, λ, t) = 0 \\ Φ (q, t) = 0 \\ G_{D} (u, \overset{\cdot}{q}) = u - \overset{\cdot}{q} = 0 \end{matrix}

(21)

wherein λ is an array of constraint force and applied force, Φ is an array of algebraic equation describing constraints, and F is the dynamic differential equation of the system. We define a state vector $y_{D} = [q^{E}, u^{E}, λ^{E}]^{E}$ , then equation (21) can be simplified as

g_{D} (y_{D}, {\overset{\cdot}{y}}_{D}, t) = 0

(22)

Equation (22) can be solved step by step through (a) prediction using Taylor series and Gear integral polynomial, (b) transforming the solve of non-linear equations into approximately solve linear equations using iteration correction, (c) integral error analysis, and (d) optimization of the integral step and order of integral polynomial, then the curves of displacement, speed, acceleration, applied forces can be obtained.³¹

Power spectral-density analysis

Power spectral-density (PSD) is one of the most important presentation form of stochastic signal, and the power density spectrum has been obtained using Welch Power method. According to Wiener–Khinchin theorem, auto-PSD of the steady stochastic signal and its auto-correlation function $R_{x} (τ)$ are a pair of Fourier transform, which the PSD function $S_{x} (ω_{D})$ can be obtained by

S_{x} (ω_{D}) = \int_{- \infty}^{\infty} R_{x} (τ) e^{- j ω_{D} τ} d τ

(23)

Nonlinear analysis of simulation time series

Phase-space reconstruction

The dynamics of the time series $x_{0}, x_{1}, \dots, x_{n - 1}$ are fully captured or embedded in the m-dimensional phase space, m ≥ d_a, where d_a is the dimension of the original attractor. A vector $\vec{x_{i}}$ in the reconstructed phase space³² is constructed from the time series as follows

\vec{x_{i}} = [x_{i}, x_{i - τ}, \dots, x_{i - (m - 1) τ}]

(24)

wherein τ is the delay time. Cao’s³³ method computes E₁ and E₂ for the data set of dimension 1 up to a dimension of D_d, which is the largest embedding dimension, used for calculation. E₁ and E₂ are defined as follows

E_{1} (d_{e}) = \frac{1}{n - τ} | \sum_{i = 1}^{n - d_{e} τ} | x_{i + d_{e} τ} - x_{n (i, d_{e}) + d_{e} τ} | |

(25)

E_{2} (d_{e}) = E_{1} (d_{e} + 1) / E_{1} (d_{e})

(26)

wherein d_e is the embedding dimension; N is the number of data points; τ is the embedding delay; $x_{i + d_{e} τ}$ and $x_{n (i, d_{e}) + d_{e} τ}$ are the ith vector in the data sets and its nearest neighbors of d_e-dimensional phase space, respectively.³⁴

Largest Lyapunov exponent

The basic characteristics of chaotic motion are that the movement is extremely sensitive to initial conditions, two very close initial values resulting in orbit over time by separating exponentially, Lyapunov exponent^35,36 that describes the amount of this phenomenon. We use the algorithm of Rosenstein et al. to calculate the Largest Lyapunov exponent (LLE). Consider the representation of the time series data as a trajectory in the embedding space and assume that observe a very close return $s_{n'}$ to a previously visited point s_n. Then, consider the distance $Δ_{0} = s_{n} - s_{n'}$ as a small perturbation, $Δ l_{s} = s_{n + l} - s_{n' + l}$ . If one finds that $| Δ l_{s} | \approx Δ_{0} e^{λ_{l} l_{s}}$ , then λ_l is the LLE. Assuming $S (ξ, m, t)$ exhibits a linear increase with identical slope for all m larger than some m₀ and for a reasonable range of $ξ$ and then this slope can be taken as an estimate of the largest exponent

S (ξ, m, t) = {\ln (\frac{1}{u_{n}} \sum_{s_{n'} \in u_{n}} | s_{n + l} - s_{n' + l} |)}

(27)

Results

Dynamics of eccentric block without torsional vibration

Before analysis the dynamics of eccentric block under the effect of torsional vibration, the dynamics of eccentric block without torsional vibration must be simulated first, and results under different deviation angles from 3° to 90° using the same constant angular speed Ω = 10 rad/s are shown in Figure 12. Because the Lagrangian approach is selected as the simulated method, its numerical solution must be different from analytical solution calculated by equation (8), the results comparison between the simulated and the calculated values of deflection angle are shown in Table 3. When the deviation angle increased from 3° to 90°, the simulated values of $φ_{t}$ decreased from 25.32° to 1.28° and the errors percentage of all contrast values is about 4%, which proves the dynamic simulated process using the models in this article in line with actual drilling situation.

Figure 12.

Simulation results under different deviation angles when the angular speed was 10 rad/s.

Table 3.

Results comparison between the simulated and the calculated values.

β (°)	Calculated value of $φ_{t}$ (°)	Simulated value of $φ_{t}$ (°)	Percentage error (%)
3	24.3259	25.3185	4.0804
20	3.6139	3.7531	3.8512
45	1.7471	1.8133	3.7875
90	1.2653	1.2812	3.7149

Under the force of friction, the motion of eccentric block is a type of damped pendulum action as shown in Figure 12. Here, we define a steady state, after the vibration amplitude damping to less than 0.1°, the steady state of eccentric block starts going on. The time required for the pendulum action damped into the steady state, we defined it as damping time t_d. Values of t_d under different deviation angles are shown in Figure 13. As β increases from 3° to 90°, the value of t_d sharply reduced from 638.5 to 25.8 s and the maximum deflection angle significantly increases from −58.2° to −1.79°, due to the larger stability acquired by larger $T_{Emax}$ of eccentric block. If the deviation angle is reduced lesser than 2°, the $T_{Emax}$ cannot overcome the friction resistance moment of plate valve, then eccentric block will rotate around and will never tend to a steady state, the deflection angle of eccentric block will continuously change between −90° and 90°. Thus, the value change of damping time t_d and maximum deflection angle $φ_{\max}$ presents a characteristic of hyperbolic cotangent. Through curve fitting, the estimation formula along deviation angle can be given as

t_{d} = 0.02785 \times \coth (1.478 e - 05 \times β) + 10.46 \frac{1}{2}

(28)

φ_{\max} = - 1.767 \times \coth (0.0112 \times β) + 0.4419

(29)

Figure 13.

Damping time and maximum deflection angle under different deviation angles.

Dynamics of eccentric block affected by vibratory amplitude

The impact of different vibratory amplitudes on deflection behavior of eccentric block when β is 3° as shown in Figure 14. After torsional vibration occurs, the swing amplitude of deflection angle no longer presenting as single attenuation wave, and the swing of eccentric block gradually transforms into a kind of pendulum motion with a constant period stimulated by the additional harmonic wave induced by borehole bottom. In addition, the swing of eccentric block is independent of the amplitude of torsional vibration. As shown in Figure 14, it is obvious that the larger intensity of torsional vibration, the larger swing amplitude of deflection angle. As the value of B increases from 3, 6 to 10, the relevant swing amplitude of deflection angle is, respectively, 2.5°, 14.4°, and 44.9°, and the time spent before the periodic pendulum motion reduces from about 168, 32 to 20 s. However, the torsional vibration is wild or severe, and the pendulum motion of eccentric block swings around a constant intermediate value of about 24.2° all the time. It is slightly smaller than simulated value of $φ_{t}$ , which is 25.3° without torsional vibration. This means the torsional vibration is beneficial to hold the eccentric block swing around a less deflection angle.

Figure 14.

Impact of different vibratory amplitudes on deflection behavior of eccentric block with deviation angle of 3°.

Figure 15 is the PSD analysis of the impact of different vibratory amplitudes on deflection behavior of eccentric block under the deviation angle of 3°. Each curve showed the same characteristics in the high frequency, because the high-frequency component of the spectrum analysis mainly comes from the random noise, which is neglected in this article. When the BHA rotates without torsional vibration, there is only one significant peak at point A with a characteristic frequency of 9.5 × 10⁻³π rad as shown in Figure 15, which is consistent with the swinging deflection angle of eccentric block. The simple damping motion of the angle will gradually approach to a constant value. After torsional vibration occurs, it is obvious that the phenomenon of fractal appears. There are many peaks along the PSD curves, when the intensity of torsional vibration increases, the general power of characteristic frequencies grows obviously and the number of significant peaks increases as the distribution of effective peaks lasts wider along normalized frequency axis. With the increase in the value of amplitude B gradually from 0, 3, 6 to 10, the number of significant peaks of each curve is 2, 3, 6, and 6, and the maximum power of characteristic peak is, respectively, 13.1, 1.24, 14.17, and 24.1 dB at points A, B, C, and D in Figure 15. The impact of vibratory amplitude on deflection behavior can be seen as a gradual change process along with the amplitude growth. The single damping motion without torsional vibration tends to gradually transform to a periodic motion with the increase in intensity of torsional vibration.

Figure 15.

PSD analysis of the impact of different vibratory amplitudes on deflection behavior of eccentric block under β of 3°.

Dynamics of eccentric block affected by vibratory period

On basis of analysis above, when the value of vibratory amplitude is 10, the intensity of pendulum motion of eccentric block will be the largest, and will fully transform the periodic motion in short time, thus the severe torsional vibration is selected to research the impact of different vibratory amplitudes on deflection behavior of eccentric block, as shown in Figure 16, which are the simulation results acquired under the deviation angle β is 12°. Under this circumstance, a larger eccentric torque can provide a better stability of eccentric block and obtain curves with smaller amplitude. When the value of ω increases, the frequency of vibration becomes high, the curve of swing motion of eccentric block tends to transform to a sine-shape wave, and the transition time before steady period motion turns longer. When ω is 0.5 and 1, respectively, there is another subordinate period in one cycle of the curve line. During one subordinate period, the swing motion of eccentric block is around the deflection angle of 6.25°, the number of peaks in one subordinate period is 6 and 3, respectively, as ω is 0.5 and 1, which is depended on the frequency of the excited wave of borehole bottom. After the value of ω increases larger than 2, the subordinate period disappeared while the amplitude of vibration turns larger. On the scale of one period, the calculated value is 5.95 according to equation (8), set the percentage error is 4, and the simulated value of $φ_{t}$ can be estimated as 6.19 without torsional vibration. When the value ω increased from 0.5 to 3, the average deflection angle is, respectively, 4.73, 4.69, 4.65, and 4.93. All of them are much smaller than the simulated value of $φ_{t}$ when compared with the deviation angle of 3°. This means the benefit of holding the eccentric block swings around a less deflection angle is strengthened due to a larger eccentric torque.

Figure 16.

Impact of different vibratory periods on deflection behavior of eccentric block under β of 12°.

The PSD analysis of these four time series are shown in Figure 17, the difference between them is much larger when compared with Figure 15, which means the impact of vibratory period on deflection behavior is larger than the impact of different vibratory periods. When the value of ω gradually increases from 0.5, 1, 2 to 3, the power value at low normalized frequency presents the phenomenon that first decreases and then increases, while at high normalized frequency presents the different phenomenon that monotonically decreases. It is obvious that with the increase in ω, the peak distribution of characteristic frequencies becomes more and wider due to the existence of subordinate period mentioned above. After the subordinate period disappeared when ω is 2 and 3, at 0.02 πrad along normalized frequency, the power of characteristic peak increases to the maximum value, respectively, is 15.62 and 22.76 dB at point C and D; the number of peaks decreases obviously and the peaks get closer with each other as well as the power grows significantly, which means the difference value between the maximum amplitude and second major amplitude of deflection angle gradually decreases to zero. Thus, an inference can be made here, if the value of ω is large enough, the swing of eccentric block can be eventually involved to a steady pendulum motion with a short period.

Figure 17.

Power-density spectrum analysis of the impact of different vibratory periods on deflection behavior of eccentric block under β of 12°.

Classification of eccentric block dynamics

Through comparative analysis of simulating results, the defection curve of eccentric block under the impact of torsional vibration can be generalized into four typical examples as shown in Figure 18. They can be represented as follows: (a) the work of eccentric torque on eccentric block is less than the work of plate friction, the eccentric block rotates around its axle, and the deflection angle is reciprocated between −90° and 90°. (b) The work of eccentric torque is larger than the friction work, the eccentric block does not rotate anymore and swings around its critical angle with a steady period, which can be seen as a kind of quasi-sine pendulum motion. (c) The motion of eccentric block is still swinging around a critical angle, but the motion curve is no longer along a curve like a single sine shape during one period, while the shape of a period curve is superimposed by two or more sine curves. (d) With the eccentric torque increasing further, if the stability of eccentric block grows enough, a significant subordinate period will appear in one period of motion curve. These four typical motion curves are expressed, respectively, by the position of mass center along X-axis, and the angular speed is shown in Figure 19; in them, Figure 19(a) shows the steady state of the curve shown in Figure 18 for a better exhibition, and other three figures show the overall process including transition section. If the eccentric block rotates around its axle, the motion of it is a kind of varying accelerated motion under the action of eccentric torque and friction torque, the angular speed will accelerate to a relative high level and the MVDT loses the deviation correction ability at this time. When the quasi-sine pendulum motion occurs, the shape of angular speed curve along X-axis position is close to circular, which presents the quasi-sine pendulum motion. With the single sine pendulum motion disappearing, the circular shape will be fractal to two or more ellipses with different orbits. After the significant subordinate period occurs, the shape of angular speed curve along the X-axis position becomes irregular, and as the ratio of subordinate period in one period increases, the irregularity increases.

Figure 18.

Four kinds of deflection behavior expressed by time as X-axis and deflection angle as Y-axis.

Figure 19.

Four kinds of deflection behavior by position as X-axis and angular speed as Y-axis.

Nonlinear analysis of simulation time series

A result can be obtained through the above analysis that the $T_{Emax}$ directly impacts the mechanism of eccentric block response, the value of $T_{Emax}$ is depended on the design parameters and the deviation angle. The value of $T_{Emax}$ increases by increasing deviation angle and length equivalently if the cross section of eccentric block is decided. During the vertical drilling process, if the value of LLE is positive, the swing motion of eccentric block maybe chaos, the motion cannot be predicted. Otherwise, if the value of LLE is small than 0, the swing motion of eccentric block is a simple periodic motion or a damping motion, which we can make an accurate forecast. The non-linear analysis of the quasi-periodic state motion of more than 30 simulating results is carried out through equations (24)–(27), and the curve between LLE, the $T_{Emax}$ , and torsional vibration is shown as follows, in which the impact of torsional vibration is expressed, respectively, by amplitude (see Figure 20) and period (see Figure 21) as well as the four typical motions marked refer to Figure 17. For all simulate results, the maximum and minimum value of LLE is 3.99 × 10⁻² and −2.98 × 10⁻², which means the motion of eccentric block under the impact of torsional vibration has both chaotic and linear characteristics. It is obvious that there is no simple rule between LLE $T_{Emax}$ and torsional vibration, but in Figure 20 when the value of $T_{Emax}$ is larger than about 20 N m, the value of LLE is all smaller than 0. After the $T_{Emax}$ decreases to less than about 20 N m, the LLE is larger than 0 only in two spaces with the value of B less than about 2 or larger than about 5. There are three spaces that the LLE is positive in Figure 21, respectively, when the value of ω is approximately less than 1 or larger than about 2.5 while the $T_{Emax}$ is roughly less than about 20 N m, and when the ω is less than about 0.5 while the $T_{Emax}$ is larger than about 50. The influence law between LLE, $T_{Emax}$ , and torsional vibration is not clear till now, and further study is needed afterward.

Figure 20.

3D surface of LLE affected by amplitude of torsional vibration.

Figure 21.

3D surface of LLE affected by period of torsional vibration.

Discussion and conclusion

Due to the working principle of MVDT, the upper plate valve and the eccentric block are rigid connected and the dynamics of eccentric block presents obvious non-linear behavior under the effect of torsional vibration, by which the control precision and inclining preventing performance of MVDT are directly determined. The effect mechanism can be generalized as follows:

Both amplitude and period of torsional vibration impact the control response behavior of MVDT, with the $T_{Emax}$ of the eccentric block increasing, the extent of these impact decreases because the larger value of $T_{Emax}$ , the better stability of the eccentric block, and the more resistance to torsional vibration. Therefore, to design a MVDT with high performance, the $T_{Emax}$ must be improved using a longer length or a larger radius, as well as selecting materials with larger density, like lead or tungsten, to produce the eccentric block.

With the increase in the amplitude of torsional vibration, the swing amplitude of deflection angle increases. While the frequency of torsional vibration decreases, a subordinate period gradually appears after the value of ω is less than 1. Once the torsional vibration occurs, the average deflection angle in one periodic motion is slightly smaller than the deflection angle without torsional vibration, and transition time is less before entering the steady periodic state. However, a torsional vibration with large amplitude and fast frequency is detrimental to the stability of the eccentric block. If the torsional vibration is controllable, then small amplitude and slow frequency are the ideal state for acquiring a better response of eccentric block.

Through analyzing more than 30 groups of simulation results, the state of torsional vibration is that the response of eccentric block can be generalized in four typical motions, which are rotate motion, quasi-sine pendulum motion, swing motion without subordinate period, and swing motion with subordinate period. When $T_{Emax}$ of eccentric block or the amplitude and period of torsional vibration changed, the deflection motion of eccentric block will be switched between these four typical motions.

The LLE of time series of simulation results are calculated, and the non-linear analysis is carried out. The maximum and minimum values of LLE are 3.99 × 10⁻² and −2.98 × 10⁻², which means the response of eccentric block would be transformed to chaos with the state of torsional vibration changing, especially in spaces of $T_{Emax}$ less than about 20 N m when the amplitude of torsional vibration is too large or the period of torsional vibration is too small. If the deflection motion is chaos, the accurate position of eccentric block cannot be predicted explicitly.

Based on above researches, a suggestion for engineering can be given, just taking the ideal uniform rotation of BHA into consideration is inadequate during the design of MVDT, the deviation control performance of the eccentric block must be tested and evaluated under the effect of torsional vibration. In this article, some laws of the dynamics of eccentric block under the torsional vibration are explored, but there is no clear influence law among LLE, $T_{Emax}$ , and torsional vibration. The intrinsic reason of this complicate non-linear phenomenon has not been completely explained yet, and further studies are needed to explore the interactions of LLE, $T_{Emax}$ , and torsional vibration.

Footnotes

Appendix 1

Handling Editor: Elsa de Sa Caetano

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by Beijing Natural Science Foundation (3174054), Natural Science Foundation of China (51704264), National Key R&D Program of China (2016YFE0202200), and Special scientific research fund for public welfare industry of ministry of land and resources of the People’s Republic of China (No. 201411054).

ORCID iD

Qilong Xue

References

Oppelt

Chur

Feld

et al . New concept for vertical drilling of boreholes (SPE-21905). Richardson, TX: Society of Petroleum Engineers, 1991.

Chur

Oppelt

. Vertical drilling technology: a milestone in directional drilling (SPE-25759). In: Proceedings of the SPE/IADC drilling conference, Amsterdam, 22–25 February 1993. Richardson, TX: Society of Petroleum Engineers.

Bram

Draxler

Hirschmann

et al . The KTB borehole: Germany’s superdeep telescope into the earth’s crust. Oilfield Rev 2004; 7: 4–22.

Ligrone

Oppelt

Calderoni

et al . The fastest way to the bottom: straight hole drilling device—drilling concept, design considerations, and field experience (SPE-36826). Richardson, TX: Society of Petroleum Engineers, 1996.

Calderoni

Savini

Treviranus

. Outstanding economic advantages based on new straight-hole drilling device proven in various oilfield locations (SPE-56444). Richardson, TX: Society of Petroleum Engineers, 1999.

Janwadkar

Fortenberry

Dawkins

et al . Pioneering advanced directional drilling technologies overcome challenges of S and J type wells in N. America (SPE-103198). In: Proceedings of the SPE annual technical conference and exhibition, San Antonio, TX, 24–27 September 2006. Richardson, TX: Society of Petroleum Engineers.

Jones

Feddema

Castro

et al . Fully mechanical vertical drilling system delivers RSS performance in vertical drilling applications while providing an economical alternative to conventional rotary steerable systems set-up for vertical hold mode (SPE-178788). In: Proceedings of the IADC/SPE drilling conference and exhibition, Fort Worth, TX, 1–3 March 2016. Richardson, TX: Society of Petroleum Engineers.

Han

Zhao

et al . Development of mechanical tool for automatic vertical drilling. Acta Petrol Sin 2008; 29: 766–768.

Ghasemloonia

Rideout

Butt

. A review of drillstring vibration modeling and suppression methods. J Petrol Sci Eng 2015; 131: 150–164.

10.

Hohl

Tergeist

Oueslati

et al . Derivation and experimental validation of an analytical criterion for the identification of self-excited modes in drilling systems. J Sound Vib 2015; 342: 290–302.

11.

Cunningham

. Analysis of downhole measurements of drill string forces and motions. J Eng Ind: T ASME 1968; 90: 208–216.

12.

Belokobyl’skii

Prokopov

. Friction-induced self-excited vibrations of drill rig with exponential drag law. Sov Appl Mech 1982; 18: 1134–1138.

13.

Gulyaev

Lugovoi

Borshch

. Self-excited vibrations of a drillstring bit. Int Appl Mech 2013; 49: 350–359.

14.

Brett

. The genesis of bit-induced torsional drillstring vibrations. SPE Drill Eng 1992; 7: 168–174.

15.

Najafov

Sofiyev

Kuruoglu

. Torsional vibration and stability of functionally graded orthotropic cylindrical shells on elastic foundations. Meccanica 2013; 48: 829–840.

16.

Sheikhlou

Rezazadeh

Shabani

. Stability and torsional vibration analysis of a micro-shaft subjected to an electrostatic parametric excitation using variational iteration method. Meccanica 2013; 48: 259–274.

17.

Liu

Gao

. A nonlinear dynamic model for characterizing downhole motions of drill-string in a deviated well. J Nat Gas Sci Eng 2017; 38: 466–474.

18.

Pelfrêne

Sellami

Gerbaud

. Mitigating stick-slip in deep drilling based on optimization of PDC bit design (SPE-139839). In: Proceedings of the SPE/IADC drilling conference and exhibition, Amsterdam, 1–3 March 2011. Richardson, TX: Society of Petroleum Engineers.

19.

Ertas

Bailey

Wang

et al . Drillstring mechanics model for surveillance root cause analysis and mitigation of torsional and axial vibrations. SPE Drill Completion 2013; 29: 1064–1071.

20.

Ahmadi

Altintas

. Stability of lateral, torsional and axial vibrations in drilling. Int J Mach Tool Manu 2013; 68: 63–74.

21.

Wang

Xue

Han

et al . Torsional vibration analysis of push-the-bit rotary steerable drilling system. Meccanica 2014; 49: 1601–1615.

22.

Lian

Zhang

Lin

et al . Experimental and numerical study of drill string dynamics in gas drilling of horizontal wells. J Nat Gas Sci Eng 2015; 27: 1412–1420.

23.

Greenwood

Williamson

JBP

. Contact of nominally flat surfaces. P Roy Soc Lond A Mat 1966; 295: 300–319.

24.

Majumdar

Bhushan

. Fractal model of elastic-plastic contact between rough surfaces. J Tribol: T ASME 1991; 113: 1–11.

25.

Liu

Chen

Wang

et al . Modeling friction-reducing performance of an axial oscillation tool using dynamic friction model. J Nat Gas Sci Eng 2016; 33: 397–404.

26.

Graff

. Wave motion in elastic solids. Columbus, OH: Ohio State University Press, 1975.

27.

Hubbard

West

. Differential equations: a dynamical systems approach. Berlin: Springer-Verlag, 1991.

28.

Martins

JAC

Oden

. A numerical analysis of a class of problems in elastodynamics with friction. Comput Method Appl M 1983; 40: 327–360.

29.

Kikuuwe

Takesue

Sano

et al . Fixed-step friction simulation: from classical coulomb model to modern continuous models. In: Proceedings of the 2005 IEEE/RSJ international conference on intelligent robots and systems, Edmonton, AB, Canada, 2–6 August 2005, pp.1009–1016. New York: IEEE.

30.

Doebelin

. System dynamics: modeling, analysis, simulation, design. New York: Marcel Dekker, Inc., 1998.

31.

Qiu

. Analytical mechanics. Beijing, China: China Railway Publishing House, 1998.

32.

Kennel

Brown

Abarbanel

HDI

. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys Rev A 1992; 45: 3403–3411.

33.

Cao

. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D 1997; 110: 43–50.

34.

Xue

Leung

Wang

et al . The chaotic dynamics of drilling. Nonlinear Dynam 2016; 83: 2003–2018.

35.

Wolf

Swift

Swinney

et al . Determining Lyapunov exponents from a time series. Physica D 1985; 16: 285–317.

36.

Rosenstein

Collins

Luca

CJD

. A practical method for calculating largest Lyapunov exponents (LLE) from small data sets. Physica D 1993; 65: 117–134.