State-Dependent Delay and Drill-String Dynamics

Available online at www.sciencedirect.com

ScienceDirect Procedia IUTAM 22 (2017) 31 – 38

IUTAM Symposium on Nonlinear and Delayed Dynamics of Mechatronic Systems

State-Dependent Delay and Drill-String Dynamics Xie Zhenga , Balakumar Balachandrana,∗ a Department

of Mechanical Engineering, University of Maryland, College Park, MD 20742, USA

Abstract In this work, stability analysis carried out with a reduced-order model of a drill string is presented. This model, which is nonlinear in nature and contains a state-dependent delay, is used to study coupled axial and torsion dynamics of the system. The axial penetration rate and drill spin speed are used as control parameters, and the analysis reveals that the consideration of the state-dependent delay is critical for capturing the influence of the axial penetration rate. As a part of the analysis, the linearized system associated with state-dependent delay is constructed and the D-subdivision method is utilized. The results are presented in the form of stability charts in the plane of spin speed and cutting coefficient and in the plane of spin speed and penetration rate. These results confirm the earlier findings obtained through numerical means with full and reduced-order models. c 2017 © Published by by Elsevier B.V.B.V. This is an open access article under the CC BY-NC-ND license  2017The TheAuthors. Authors. Published Elsevier (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of organizing committee of the IUTAM Symposium on Nonlinear and Delayed Dynamics of Peer-review under responsibility of organizing committee of the IUTAM Symposium on Nonlinear and Delayed Dynamics of Mechatronic Systems.

Mechatronic Systems

Keywords: Drill sting; State-dependent delay; Stability analysis; D-subdivision method

1. Introduction Rotary drilling systems are used in oil and gas exploration operations. Long and slender flexible structures called drill strings are key components of these systems. These rotating, flexible structures often experience stick-slip oscillations (e.g., Liao, Balachandran, Karkoub, and Abdel-Magid 1 ). As shown in earlier work (Richard, Germay, and Detournay 2,3 ), the time-delay effect associated with bit-rock cutting mechanics can play an important role in causing the stick-slip behavior. Drill-string system models may be constructed through a variety of discretization means (e.g., references 1,4,5,6 ). In earlier work 6,7 conducted in the authors’ group, reduced-order models, finite-element based discretization, and the presence of the state-dependent delay have been discussed. In the present work, with the goal of further understanding the effect of this delay, the authors have used the reduced-order model presented in the earlier work of Besselink, Deonoel, and Nijmeijer 4 and Liu, Vlajic, Long, Meng, and Balachandran 7 . The remaining sections of this article are organized as follows. In Section 2, a reduced-order system with two degrees of freedom for describing the drill-string dynamics is presented. With the goal of reducing the number of parameters, the governing equations are rewritten in nondimensionalized form. Next, the linearization is carried out in Section 3 as a step towards stability analysis. In Section 4, the obtained linearized system is analyzed by using ∗

B. Balachandran. Tel.: +1-301-405-4747 ; fax: +1-301-314-9477. E-mail address: [email protected]

2210-9838 © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of organizing committee of the IUTAM Symposium on Nonlinear and Delayed Dynamics of Mechatronic Systems doi:10.1016/j.piutam.2017.08.006

32

Xie Zheng and Balakumar Balachandran / Procedia IUTAM 22 (2017) 31 – 38

the D-subdivision method and the illustrative stability charts obtained are presented. Finally, concluding remarks are collected together and presented at the end. 2. Modeling and Nondimensionalization In Figure 1, an illustrative model of a drill-string system is illustrated. The axial and torsion motions of interest are also shown. At the top end of the drill sting, the system is imposed with a constant axial speed V0 and a rotation speed Ω0 . In terms of the axial displacement Z(t) and the rotational displacement Φ(t), the governing equations take the form ¨ + Ca Z(t) ˙ + Ka (Z(t) − V0 t) = W s − Wb (t) M Z(t) (1) ¨ + Ct Φ(t) ˙ + Kt (Φ(t) − Ω0 t) = −T b (t) I Φ(t) Here, M and I are the respective translational and rotational inertias, Ka and Kt are the respective translational stiffness and torsion stiffness, and Ca and Ct represent the respective translational damping and torsion damping. Furthermore, W s is the sum of the weight of both the drill pipe and drill collar. Wb and T b respectively denote the weight and torque on the bi, and they are both determined by bit-rock interactions. Each of them can be decomposed in terms of cutting and friction components, as follows. Wb (t) = Wbc (t) + Wb f (t) (2) T b (t) = T bc (t) + T b f (t) The subscripts bc and b f are used to denote the cutting and friction components, respectively. Following the earlier work of Detournay and Defourny 8 , those components can be expressed as ˙ Wbc (t) = aζR(d(t))H(Φ(t)) ˙ T bc (t) = 12 a2 R(d(t))H(Φ(t)) ˙ Wb f (t) = σalH(d(t))H(Z(t)) ˙ ˙ T b f (t) = 12 μγa2 σlsgn(Φ)H(d(t))H( Z(t))

(3)

where the R(.) function is the unit ramp function and H(.) is the Heaviside step function. The different parameters along with representative values are listed in Table 1. A drag bit or polycrystalline diamond compact (PDC) bit is often used in drilling operations. This bit is assumed to have N identical blades that are symmetrically distributed. In Figure 2, two successive blades are shown along with the delayed states. For an individual blade, the instantaneous depth of cut can be determined as

Fig. 1. Representative reduced-order model of drill-string system, following prior work 4,7 .

33

Xie Zheng and Balakumar Balachandran / Procedia IUTAM 22 (2017) 31 – 38

Fig. 2. Two successive blades of a drill bit.

Table 1. Parameters used for drilling operations (values adopted from references 3,7 ) Parameter Mass Axial damping Axial stiffness Moment of inertia Torsion damping Torsion stiffness Radius of drill bit Wear flat length Intrinsic specific energy of rock Contact strength Cutter face inclination Friction coefficient Geometry parameter of drill bit Number of blades on drill bit

Symbol

Value

Unit

M Ca Ka I Ct Kt a l  σ ζ μ γ N

3.4 × 104

kg N s/m N/m kg m2 N s m/rad N/m m m MPa Mpa -

1.56 × 104 7.0 × 105 116 32.9 938 0.108 0.0012 0 − 110 60 0.6 0.6 1 4

dn (t) = Z(t) − Z(t − τ)

(4)

Assuming that the cutting action is uniform across the N blades, then, the cutting depth is d(t) = Ndn (t)

(5)

where the delay τ is given by

2π (6) N The state-dependent delay τ(Φ(t)) is the elapsed time for the drill bit to rotate over an angle of 2π N , and this delay depend only on the state Φ. Next, form. Following earlier work 4 , the characteristic time √ the equations of motion are cast into dimensionless 2 t∗ = I/Kt and characteristic length L∗ = 2Kt /a are introduced. Then, one can obtain the nondimensional variables as Φ(t) − Φ(t − τ) =

z=

Z − Z¯ L∗

¯ ϕ=Φ−Φ

tˆ = t/t∗

τˆ = τ/t∗

(7)

34

Xie Zheng and Balakumar Balachandran / Procedia IUTAM 22 (2017) 31 – 38

¯ correspond to the equilibrium solution of Eq.(1), which is a trivial solution in the absence of vibrations. ¯ and Φ Here, Z, The axial state z and angular state φ are functions of dimensionless time tˆ. With the nondimensional variables, the governing equations can be recast as x˙ = f(x(tˆ), x(tˆ − τˆ (x)), τˆ (x))

(8)

where  T x = z ϕ z˙ ϕ˙ ⎞ ⎛ z˙(tˆ) ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ ⎜⎜⎜ ϕ( ˙ tˆ) ⎟⎟ f = ⎜⎜⎜⎜ 2 ⎜⎜⎝−2ξη˙z(tˆ) − η z(tˆ) − Nψ(z(tˆ) − z(tˆ − τˆ )) − Nψv0 ω0 τˆ ⎟⎟⎟⎟⎠ −2κϕ( ˙ tˆ) − ϕ(tˆ) − N(z(tˆ) − z(tˆ − τˆ )) − Nv0 ω0 τˆ

(9)

and xt (s) = x(t + s), s ∈ [-r 0], r ∈ R+ . The dimensionless form of equation (6) takes the form ω0 τˆ + ϕ(tˆ) − ϕ(tˆ − τˆ ) =

2π N

(10)

The dimensionless parameters are defined as Ca Ka 2 Ka I Ct t = κ= √ η2 = √ M ∗ Kt M 2 Ka M 2 Kt I aζI V0 /L∗ a2 ψ= ω0 = Ω0 /t∗ v0 = = V0 Kt M Ω0 2Kt Ω0 ξ=

(11)

3. Construction of Associated Linear System According to the work of Hartung 15 , the true linearization of system with state-dependent delay is not possible due to the fact that the solution of the system is not differentiable with respect to state-dependent delay. Hence, one needs to find a constant delay model which has the same local stability properties as the original system. Making use of the method discussed by Hartung 9 , the linear system associated with the constant solution x(t) ≡ x¯ is written as y˙ = D1 f(¯x, x¯ , τ(¯xt ))y(t) + D2 f(¯x, x¯ , τ(¯xt ))y(t − τ(¯x)) + D3 f(¯x, x¯ , τ(¯xt ))Dτ(¯xt )yt

(12)

where, D1 , D2 , and D3 are the derivatives with respect to the 1st, 2nd, and 3rd argument of f, and Dτ denotes the Fr´echet derivative of time delay τ with respect to yt .The y vector is defined as  T y = u θ u˙ θ˙

(13)

where u = z − z¯ and θ = ϕ − ϕ¯ are the perturbations around the constant solutions (case of no vibrations). ⎛ ⎞ 0 0 1 0 ⎟⎟ ⎜⎜⎜ ⎟ ⎜⎜⎜ 0 0 0 1 ⎟⎟⎟⎟ ⎟ D1 f(¯x, x¯ , τ(¯xt )) = ⎜⎜⎜⎜ 2 ⎜⎜⎝−η − Nψ 0 −2ξη 0 ⎟⎟⎟⎟⎠ −N −1 0 2κ ⎛ ⎜⎜⎜ 0 ⎜⎜⎜ 0 D2 f(¯x, x¯ , τ(¯xt )) = ⎜⎜⎜⎜ ⎜⎜⎝Nψ N

0 0 0 0 0 −2ξη 0 0

⎞ 0⎟⎟ ⎟ 0⎟⎟⎟⎟ ⎟ 0⎟⎟⎟⎟⎠ 0

(14)

(15)

35

Xie Zheng and Balakumar Balachandran / Procedia IUTAM 22 (2017) 31 – 38

and

⎞ ⎛ 0 ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ ⎜⎜⎜ 0 ⎟⎟ ⎜ D3 f(¯x, x¯ , τ(¯xt )) = ⎜⎜⎜ ⎜⎝⎜−Nψv0 ω0 ⎟⎟⎟⎠⎟ −Nv0 ω0

(16)

After taking the Fr´echet derivatives of both sides of Eq.(10), one arrives at Dτ(¯xt )yt = −

θ(tˆ) − θ(tˆ − τ¯ ) ω0

(17)

where τ¯ = 2π/(Nω0 ) is the constant delay for motion without vibrations. On substituting Eqs.(14)-(17) into Eq.(12), the resulting linearized system is u¨ (tˆ) + 2ξη˙u(tˆ) + η2 u(tˆ) = −Nψ(u(tˆ) − u(tˆ − τ¯ )) + Nψv0 (θ(tˆ) − θ(tˆ − τ¯ )) ¨ tˆ) + 2κθ( ˙ tˆ) + θ(tˆ) = −N(u(tˆ) − u(tˆ − τ¯ )) + Nv0 (θ(tˆ) − θ(tˆ − τ¯ )) θ(

(18)

For the previously discussed constant delay case, the linearized equations take the form u¨ (tˆ) + 2ξη˙u(tˆ) + η2 u(tˆ) = −Nψ(u(tˆ) − u(tˆ − τ¯ )) ¨ tˆ) + 2κθ( ˙ tˆ) + θ(tˆ) = −N(u(tˆ) − u(tˆ − τ¯ )) θ(

(19)

Through comparison, the system with the state-dependent delay is seen to have additional terms related to torsion motion, which is an implicit function of the state-dependent delay. These additional terms are dependent on the axial penetration rate v0 . The influence of this delay on the stability of the system will be investigated in the next section. 4. Stability analysis and results The D-subdivision method is a commonly used approach for linear autonomous systems with time delays. The main idea of the D-subdivision method can be described as follows: first, one defines the stability crossing set in the parameter space where characteristic equation has imaginary roots. Then, by studying the crossing direction of imaginary roots of the crossing set, the numbers of right-half-plane roots are determined. For more details about D-subdivision method see references 12,14 . For the current system, after a rather lengthy calculation, the authors arrive at the following characteristic equation: P0 (s) + P1 (s)(1 − e−¯τ s ) = 0

(20)

Here, P0 , and P1 are polynomials in the eigenvalue s and these polynomials take the form P0 (s) = s4 + (2ξη + 2κ)s3 + (η2 + 4κξη + 1)s2 + (2ξη + 2κη2 )s + η2 P1 (s) = (Nψ − Nv0 )s2 + (2κNψ − 2ξηNv0 )s + (Nψ − Nη2 v0 )

(21)

After substituting s = iω and τ¯ = 2π/ω0 into Eqs.(21) and separating the real and imaginary parts, one obtains the stability crossing set in the ω0 -ψ domain: 1 α 2 + β2 + (α(ω2 − η2 ) + 2βξηω)Nv0 ] [ 2 N[α(ω2 − 1) − 2βκω] 2πω , k = 1, 2, ..., = N(Θ1 + (2k − 1)π)

ψS DD = ω0S DD Here,

α = Real(P0 )

β = Imag(P0 )

Θ1 = ∠

−P1 P0 + P1

(22)

(23)

36

Xie Zheng and Balakumar Balachandran / Procedia IUTAM 22 (2017) 31 – 38

Fig. 3. Stability charts in the plane of drive speed ω0 and cutting coefficient ψ for different penetration speeds v0 .

From Eqs. (19), which form the linearized equation for constant delay case, it is clear that the system is independent of the axial penetration rate v0 . Therefore, when v0 = 0, Eqs.(21) should degenerate to the characteristic equation for the constant delay model; that is, P0 (s) + P2 (s)(1 − e−¯τ s ) = 0

(24)

P2 (s) = Nψ(s2 + 2κs + 1)

(25)

Setting v0 in P1 to zero, one obtains

Following the same procedure as before, the stability boundary for constant delay model is given by 1 α 2 + β2 N[α(ω2 − 1) − 2βκω] 2 2πω = , k = 1, 2, ..., N(Θ2 + (2k − 1)π)

ψCD = ω0CD and

Θ2 = ∠

−P2 P0 + P2

(26)

(27)

For the parameters listed in Table 1, the stability chart generated in the ω0 − ψ plane is shown in Figure 3. In this chart, the stability boundaries obtained for v0 = 0, 0.1, 1, and 2 are given respectively. The stability chart for v0 = 0 is equivalent to that for the constant delay case. It can be seen that the left boundaries for both the constant delay (v0 = 0) and state-dependent delay cases are almost identical to each other. With an increase in v0 , the right boundary for the chart of the state-dependent delay model is shifted upwards. By contrast, the right boundary of the stability chart for the constant delay remains the same and independent of the penetration rate v0 . It is clear that the consideration of the state-dependent delay is important to capture the reduction in the stability region in the considered ω0 − ψ plane. Furthermore, from these results, it can be discerned that for stability purposes, a higher cutting coefficient ψ is preferrable. The obtained results are consistent with those obtained in prior work conducted in the authors’ group 6,7 . Here, the linearization construction and analytical construction of the stability boundaries allows for more convenient generation of the stability charts. Next, another stability chart is provided in the parameter space of penetration rate v0 and drive speed ω0 , which are two critical control parameters during drill operations. By using the same approach as before, the stability boundary for state-dependent delay model can be derived as the following:

37

Xie Zheng and Balakumar Balachandran / Procedia IUTAM 22 (2017) 31 – 38

Fig. 4. Stability charts in the plane of the drive speed ω0 and the penetration speed v0 , for different values of ξ and κ.

1 α 2 + β2 + (α(1 − ω2 ) + 2βκω)Nψ] [ 2 2 − η ) + 2βξηω] 2πω = , k = 1, 2, ..., N(Θ1 + (2k − 1)π)

v0S DD = ω0S DD

N[α(ω2

(28)

Unlike the state-dependent delay case, for the constant delay case, one can only get the stability boundary in terms of spin rate ω0 , which is a function of the crossing frequency ωc . The boundary is determined from ω0CD =

2πωc N(Θ2 + (2k − 1)π)

(29)

For the crossing frequency, the obtained expression is

2 +β2 κβ + κ2 β2 − α(α + α2Nψ ) ωc = (30) α Based on the expressions obtained, the stability charts obtained for the state-dependent delay model in the space of the drive speed ω0 and the axial penetration speed v0 are shown in Figure 4, for different damping ratios and stiffness values. The damping parameter is meant to capture the influence of drill mud, which is used in drilling operations to control drill-string dynamics. The upper boundary of the stability crossing set, is generated when the crossing frequency ω is small (typically ω < 1) and this boundary shifts upwards as the damping is increased. As previously noted in the context of the results presented in the ω0 -ψ plane, the left boundaries do not shift with respect to the increase in damping. In fact, the left boundary can be represented by Eq.(29), which corresponds to the constant delay and is associated with a large crossing frequency ω. The results indicate that damping can be influential in addressing instabilities associated with the state-dependent delay. Again, the obtained results agree well with those obtained in prior work conducted in the authors’ group 6,7 . 5. Concluding Remarks A reduced-order system descriptive of axial and torsion motions of a drill string has been considered and stability analyses has been conducted to examine the influence of the state-dependent delay. The stability boundaries for both state-dependent delay and constant delay cases are obtained by using the D-subdivision method. The analyticalnumerical results are helpful in understanding the influence of the state-dependent delay and the importance of considering it for obtaining accurate stability region predictions.

38

Xie Zheng and Balakumar Balachandran / Procedia IUTAM 22 (2017) 31 – 38

Acknowledgements Partial support received for this work from the Qatar National Research Fund for NPRP Project 7-083-2-041 is gratefully acknowledged. References 1. Liao,C., Balachandran,B., Karkoub, M., and Abdel-Magid, Y.L. (2011) Drill-string dynamics: reduced-order models and experimental studies. ASME J. Vib. Acoust. 133(4),p. 041008. 2. Richard,T., Germay, C., and Detournay, E.(2004) Self-excited stick-slip oscillations of drill bits. C.R.Mec.,332(8),pp.619-626. 3. Richard,T., Germay, C., and Detournay, E. (2007) A simplified model to explore the root cause of stickslip vibrations in drilling systems with drag bits. J. Sound Vib. 305(3),432-456. 4. Besselink,B., van de Wouw, N., and Nijmeijer, H.(2011) A semi-analytical study of stick-slip oscillations in drilling systems. ASME J. Comput. Nonlinear Dyn. 6(2),p.021006. 5. Germay,C, Denoel, V., and Detournay, E. (2009) Multiple mode analysis of the self-excited vibrations of rotary drilling systems. J. Sound and Vib. 325(12), pp. 362-381. 6. Liu X., Vlajic N., Long X., Meng G., and Balachandran B.(2014) State-dependent delay influenced drill-string oscillations and stability analysis ASME J. Vibration and Acoustics DOI: 10.11115/1.4027958. 7. Liu X., Vlajic N., Long X., Meng G., and Balachandran B. (2014) Coupled axial-torsional dynamics in rotary drilling with state-dependent delay: stability and control J. Non-linear Dyn 78:1891-1906. 8. Detournay,E., and Defourny, P.(1992) A phenomenological model for the drilling action of drag bits. Int. J. Rock Mech. Min. Sci. Geomech. Abstr.,29(1):13-23. 9. Hartung F. and Turi J. (2000) Linearized stability in functional-differential equations with stat-dependent delays. Proceeding of the conference Dynamical Systems and Differential Equations, added volume of Discrete and Continuous Dynamical Systems: 416-425. 10. Insperger T., Stepan G., and Turi J.(2007) State-dependent delay in regenerative turning processes, Nonlinear Dyn. 47(1-3): 275-283. 11. Insperger T., Barton D. A. W., and Stepan G. (2007) Criticality of Hopf bifurcation in state-dependent delay model of turning processes Int. J. Non-Linear Mechanics 43: 140-149. 12. Gu K. and Zheng X. (2014) An overview of stability crossing set for systems with scalar delay channels 33rd Chinese control conference. 13. Brokate M., and Colonius F. (1990) Linearizing equations with state-dependent delays Appl Math Optim 21:45-52. 14. Stepan, G. (1989) Retarded dynamical systems. Longman, Harlow. 15. Hartung, F. and Turi, J. (1997) On differentiability of solutions with respect to parameters in state-dependent delay equations. J. Differ. Equ. 135(2),192-237.