# Analysis of the Longitudinal-Bending-Torsional Coupled Vibration Mechanism of the Drilling of a Roof Bolter for Mine Support System

Authors: Youyu Liu; Liteng Ma; Renjun Li; Bo Chen; Yongcun Guo
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The condition of the drilling of a roof bolter for mine support system is complex, and the existing analytical approaches for the vibration of the drilling for oil wells cannot be fully applicable to it. In order to reveal the vibration mechanism of bolt drilling, some models for longitudinal, bending, and torsional coupled vibration were established based on the energy method and finite element method. A global matrix was assembled in the order of units. Taking the drilling condition as the boundary conditions, a reduced-order model was constructed. The simulations for drilling mudstone show that the order of mode shape is positively correlated with the corresponding natural frequency, and the low-order mode shape plays a decisive role in the vibration characteristics of a drill string. The order of causing drill string displacement from large to small is as follows: longitudinal vibration, torsional vibration, and bending vibration. Compared with the bit, the displacements of longitudinal and torsional vibration of the tail of the drill string are decreased by 93.2% and 84.7%, respectively. The perturbed force and reaction torque are the main reasons for inducing the coupled vibration of it. This research provides a theoretical basis for analysing the vibration characteristics of the drilling of roof bolters for mine support systems and developing their absorbers.


---

## Body

## 1. Introduction

Roof bolters are key mechanical equipment for roadway supports for coal mines, which can realize fragmentizing coal rock and drilling by the high-speed impact and rotation of drill string. There is a complex vibration in drill string in operation under the combined excitation of different geotechnical parameters, axial thrust, and torque. The main vibration forms including bending vibration, longitudinal vibration, and torsional vibration, and various vibrations are nonlinear coupled, which possibly lead to the stick-slip and bit-bounce of drill string [1]. These are the main reasons for the premature failure of drill string and affecting the support effect.Currently, there are transfer matrix method [2], lumped mass method [3], generalized coordinate method [4], hypothetical mode method [5], finite element method [6, 7], elastic line method [8], and other methods to analyse the vibration mechanism of drill string. For a drill string system with multiple pipe sections, Han et al. calculated the total vibration transfer matrix by multiplying all individual matrices; each is obtained for an individual pipe section [2]. This method simplifies the modelling process and vibration analysis. However, its essence is to solve higher-order algebraic equations. It is difficult to obtain accurate solutions while dealing with multi-DOF vibration systems. Chen et al. concentrated the masses along the drill string and proposed a kinematic coupled approach to correlate the axial displacement of the bit to the rotation angle of the bit [3]. However, this approach only discusses torsional and longitudinal coupled vibration, and the natural frequency error increases with the increase of the order of the mode shape. Cheng et al. discretized the drill string into some absolute nodal coordinate beam units and established its kinetic models of coupled vibration, which can describe the coupled vibration phenomenon of the drill string [4]. The key to this approach lies in the physical coordinates of the drill string. Moreover, and the selected mode shape must meet all the boundary conditions given by the system, which is difficult to operate and harsh to use. Jeong and Yoo proposed a hypothetical mode method to conduct the static or dynamic analysis of a flexible beam undergoing large deflection, and simplify the expressions of strain energy and geometric constraints among deformation variables [5]. However, this method uses a finite number of hypothetical modal vibrations to approximately simulate it, which requires many nodal coordinates and cumbersome calculation. Ren et al. presented a coupled model for axial/torsional/lateral vibrations of the drill string, in which the nonlinear dynamics and qualitative analysis method are employed to find out the key factors and sensitive zone for coupled vibration [6]. Deng et al. obtained the time domain response of the drill string based on the shape and strength of the geotechnical surface [9]. Ritto et al. proposed a stochastic computational model by the nonlinear Timoshenko beam theory to model uncertainties in the bit-rock interaction model [7]. However, the 3-order elastic matrix in 2D plane is used will reduce the accuracy of the whole drill string model while calculating the strain energy of it. Tucker and Wang analysed the stability of a drill string under torsional, axial, and bending coupled disturbances by the elastic line method and transformed an exact motion equation into an ordinary differential-difference equation [8]. The space does not need to be discretized by this method; thus, it is suitable for higher frequency vibration.Most of the research objects in the above work are the drill pipes for oil wells, which are essentially different from the drill strings for mine support systems, which are mainly reflected as follows: (1) there are large differences in structure and size between them. A drill pipe system for oil wells is composed of a drilling rig, a top drive system, drill pipes, a drill collar, and a drill bit. The drill pipes are connected by bolts and have a length of hundreds of meters. However, the drill string of the roof bolter is mostly a whole steel pipe wound with spiral blades, with a length of 0.5 to 3 m. (2) Their stress conditions are very different. The force loading on the drill string of the roof bolter during construction is much less than that of another.The existing analytical methods for the drilling vibration of the drill pipes for oil wells are not completely suitable for the roof bolter for the mine support system. In this paper, a dynamic model for 6 degree of freedom (DOF) bolt drilling will be established; a set of mass, stiffness, damping, and force loading matrixes will be assembled; and the order of the finite element models will be reduced. These are aimed at revealing the longitudinal-bending-torsional coupled vibration mechanism of the bolt drilling. Through the modal analysis and vibration response analysis for drilling mudstone, the law and change trend of the coupled vibration of the drill string will be explored. It is expected to provide a theoretical basis for analysing the vibration characteristics of drill strings and developing absorbers.


## 2. Kinetic Equation

### 2.1. Underlying Assumption

For the nonlinear coupled of the 3 types of vibration forms of a drill string, there are 6 assumptions as follows: (1) The drill string is a thin straight rod with equal cross section, and the cross section is circular; (2) The drill bit is regarded as a rigid body, and the body of the drill string as an elastomer with small deformation; (3) Ignoring the influence of temperature rise caused by friction between drill pipe and coal rock during the process of fragmentizing coal rock [10, 11]; (4) The direction of axial thrust is along the axis of the drill string; (5) The displacement of the bending vibration of the drill string is decomposed into x and y axes; (6) The drill string is regarded as an Euler–Bernoulli beam, namely, the effects of shear deformation and rotational inertia are ignored.


### 2.2. Kinetic Energy Equation

As shown in Figure1, the drill string system is discretized into 1 to n units, and each unit has 6-DOF. The displacement of longitudinal vibration is u; that of bending vibration are v and w; the rotation angle around x, y, and z axis is θx, θy, and θz, respectively.Figure 1
Discrete system of drill string.The translational kinetic energy of the unit of drill string (UDS) is expressed as [12](1)Tp=12∫0lρAsTsdz.According to Figure1, the translational velocity of the unit is(2)s=∂u∂t∂v∂t∂w∂tT.Substitute (2) into (1),(3)Tp=∫0lAρ2∂u∂t2+∂v∂t2+∂w∂t2dz.As shown in Figure2, the rotation velocity of the unit is expressed by Euler angle. Firstly, rotating the global coordinate system o­xyz around the y-axis by an angle ψ to get the coordinate system o­x1y1z1; secondly, rotating it again around the z1-axis by an angle θ to get the coordinate system o­x2y2z2; finally, rotating it once more around the x2-axis by an angle φ to get the coordinate system o­x3y3z3.Figure 2
Transforming the coordinate system of the unit of drill string (a) rotating aroundy-axis; (b) rotating around z1-axis; and (c) rotating around x2-axis.
(a)(b)(c)Using Euler kinematics equation [13], the rotational velocity is shown as,(4)Rs=RsxRsyRsz=∂ψ∂tsinθ+∂ϕ∂t,∂ψ∂tcosθcosϕ+∂θ∂tsinϕ,∂θ∂tcosϕ−∂ψ∂tcosθsinϕ.The rotational kinetic energy of the UDS is expressed as [12](5)Tz=12Rs⋅J⋅Rs.Substitute (4) into (5),(6)Tz=Iz2∂ϕ∂t2+2θ∂ϕ∂t∂ψ∂t+Ixy2∂ψ∂t2+∂θ∂t2.Due to the rotation angle of the drill string restricted by coal rock, there is an assumption as(7)ϕ=θx,ψ=θy,θ=θz.Substitute (7) into (6),(8)Tz=∫0lρIz2∂θz∂t2+2θz∂θy∂t∂θx∂t+ρIyx2∂θy∂t2+∂θx∂t2dz.According to simultaneous (3) and (8), the total kinetic energy of the UDS is shown as(9)T=∫0lρIz2∂θz∂t2+2θz∂θy∂t∂θx∂t+ρIyx2∂θy∂t2+∂θx∂t2dz+∫0lAρ2∂u∂t2+∂v∂t2+∂w∂t2dz.


### 2.3. Strain Energy Equation

The strain energy of the UDS is described as follows [14]:(10)H=12∫VεTσdV=12∫VεTDεdV.Due to the length of the drill string is much larger than its diameter, the stress and strain are mainly concentrated in the normal direction of the cross-section of the drill string. Its components can be simplified as(11)σy=σx=σyx=0,εy=εx=−μεz,εyx=0.Substituting (11) into (10), the strain energy of the UDS is(12)H=12∫VλEεz2+Gεzy2+εzx2dV,where λ=2μ−1μ−1/μ+1.The strains of the UDS are as follows [15]:(13)εz=∂u∂z+12∂u∂z2+∂v∂z2+∂w∂z2,εzy=12∂u∂y+∂w∂z+∂u∂y∂u∂z+∂v∂y∂v∂z+∂w∂y∂w∂z,εzx=12∂u∂x+∂v∂z+∂u∂z∂u∂x+∂v∂z∂v∂x+∂w∂z∂w∂x.According to the assumption of Euler–Bernoulli beam [16], compared with bending deformation, shear deformation can be ignored, thus,(14)θx=∂v∂z,θy=−∂w∂z.Substituting (14) into (13), and then substituting the result into (12), the linear strain energy of the UDS is as follows:(15)HL=λEA2∫0l∂u∂z2dz+GIx2∫0l∂θz∂z2dz+λEIy2∫0l∂θy∂z2dz+λEIx2∫0l∂θx∂z2dz.The nonlinear strain energy of the UDS is generated by longitudinal vibration, longitudinal-torsional coupled vibration, longitudinal-bending coupled vibration, and bending-torsional coupled vibration. The expressions of them are as follows:(16)HN1=λEA2∫0l∂u∂z3dz,(17)HN2=λEIz2∫0l∂u∂z∂θz∂z2dz,(18)HN3=λEA2∫0l∂u∂zθy2+∂u∂zθz2dz+3λEIx2∫0l∂u∂z∂θx∂z2dz+3λEIy2∫0l∂u∂z∂θy∂z2dz,(19)HN4=λE−GIz2∫0lθy∂θz∂z∂θx∂z−θx∂θz∂z∂θy∂zdz.According to simultaneous (16) and (19), the total strain energy of the UDS is shown as(20)H=HL+HN1+HN2+HN3+HN4.


## 2.1. Underlying Assumption

For the nonlinear coupled of the 3 types of vibration forms of a drill string, there are 6 assumptions as follows: (1) The drill string is a thin straight rod with equal cross section, and the cross section is circular; (2) The drill bit is regarded as a rigid body, and the body of the drill string as an elastomer with small deformation; (3) Ignoring the influence of temperature rise caused by friction between drill pipe and coal rock during the process of fragmentizing coal rock [10, 11]; (4) The direction of axial thrust is along the axis of the drill string; (5) The displacement of the bending vibration of the drill string is decomposed into x and y axes; (6) The drill string is regarded as an Euler–Bernoulli beam, namely, the effects of shear deformation and rotational inertia are ignored.


## 2.2. Kinetic Energy Equation

As shown in Figure1, the drill string system is discretized into 1 to n units, and each unit has 6-DOF. The displacement of longitudinal vibration is u; that of bending vibration are v and w; the rotation angle around x, y, and z axis is θx, θy, and θz, respectively.Figure 1
Discrete system of drill string.The translational kinetic energy of the unit of drill string (UDS) is expressed as [12](1)Tp=12∫0lρAsTsdz.According to Figure1, the translational velocity of the unit is(2)s=∂u∂t∂v∂t∂w∂tT.Substitute (2) into (1),(3)Tp=∫0lAρ2∂u∂t2+∂v∂t2+∂w∂t2dz.As shown in Figure2, the rotation velocity of the unit is expressed by Euler angle. Firstly, rotating the global coordinate system o­xyz around the y-axis by an angle ψ to get the coordinate system o­x1y1z1; secondly, rotating it again around the z1-axis by an angle θ to get the coordinate system o­x2y2z2; finally, rotating it once more around the x2-axis by an angle φ to get the coordinate system o­x3y3z3.Figure 2
Transforming the coordinate system of the unit of drill string (a) rotating aroundy-axis; (b) rotating around z1-axis; and (c) rotating around x2-axis.
(a)(b)(c)Using Euler kinematics equation [13], the rotational velocity is shown as,(4)Rs=RsxRsyRsz=∂ψ∂tsinθ+∂ϕ∂t,∂ψ∂tcosθcosϕ+∂θ∂tsinϕ,∂θ∂tcosϕ−∂ψ∂tcosθsinϕ.The rotational kinetic energy of the UDS is expressed as [12](5)Tz=12Rs⋅J⋅Rs.Substitute (4) into (5),(6)Tz=Iz2∂ϕ∂t2+2θ∂ϕ∂t∂ψ∂t+Ixy2∂ψ∂t2+∂θ∂t2.Due to the rotation angle of the drill string restricted by coal rock, there is an assumption as(7)ϕ=θx,ψ=θy,θ=θz.Substitute (7) into (6),(8)Tz=∫0lρIz2∂θz∂t2+2θz∂θy∂t∂θx∂t+ρIyx2∂θy∂t2+∂θx∂t2dz.According to simultaneous (3) and (8), the total kinetic energy of the UDS is shown as(9)T=∫0lρIz2∂θz∂t2+2θz∂θy∂t∂θx∂t+ρIyx2∂θy∂t2+∂θx∂t2dz+∫0lAρ2∂u∂t2+∂v∂t2+∂w∂t2dz.


## 2.3. Strain Energy Equation

The strain energy of the UDS is described as follows [14]:(10)H=12∫VεTσdV=12∫VεTDεdV.Due to the length of the drill string is much larger than its diameter, the stress and strain are mainly concentrated in the normal direction of the cross-section of the drill string. Its components can be simplified as(11)σy=σx=σyx=0,εy=εx=−μεz,εyx=0.Substituting (11) into (10), the strain energy of the UDS is(12)H=12∫VλEεz2+Gεzy2+εzx2dV,where λ=2μ−1μ−1/μ+1.The strains of the UDS are as follows [15]:(13)εz=∂u∂z+12∂u∂z2+∂v∂z2+∂w∂z2,εzy=12∂u∂y+∂w∂z+∂u∂y∂u∂z+∂v∂y∂v∂z+∂w∂y∂w∂z,εzx=12∂u∂x+∂v∂z+∂u∂z∂u∂x+∂v∂z∂v∂x+∂w∂z∂w∂x.According to the assumption of Euler–Bernoulli beam [16], compared with bending deformation, shear deformation can be ignored, thus,(14)θx=∂v∂z,θy=−∂w∂z.Substituting (14) into (13), and then substituting the result into (12), the linear strain energy of the UDS is as follows:(15)HL=λEA2∫0l∂u∂z2dz+GIx2∫0l∂θz∂z2dz+λEIy2∫0l∂θy∂z2dz+λEIx2∫0l∂θx∂z2dz.The nonlinear strain energy of the UDS is generated by longitudinal vibration, longitudinal-torsional coupled vibration, longitudinal-bending coupled vibration, and bending-torsional coupled vibration. The expressions of them are as follows:(16)HN1=λEA2∫0l∂u∂z3dz,(17)HN2=λEIz2∫0l∂u∂z∂θz∂z2dz,(18)HN3=λEA2∫0l∂u∂zθy2+∂u∂zθz2dz+3λEIx2∫0l∂u∂z∂θx∂z2dz+3λEIy2∫0l∂u∂z∂θy∂z2dz,(19)HN4=λE−GIz2∫0lθy∂θz∂z∂θx∂z−θx∂θz∂z∂θy∂zdz.According to simultaneous (16) and (19), the total strain energy of the UDS is shown as(20)H=HL+HN1+HN2+HN3+HN4.


## 3. Finite Element Model with 6-DOF

### 3.1. Mass Matrix

Taking any unit on the drill string to analyse, the upper and lower nodes of it are represented byj and i, respectively. The unit coordinate systems j-xyz and i-xyz are established respectively, as shown in Figure 3.Figure 3
Generalized coordinates of the nodal displacement of UDS.There are 6-DOF for each node of the discrete model of the UDS. The displacement matrix of the unit is as follows.(21)qe=uiviwiθxiθyiθziujvjwjθxjθyjθzjT.According to the displacement shape function [17], each displacement matrix of the UDS is as follows:(22)uvwθz=N100000N2000000N3N40000N5N6000000N7N80000N9N10000000N1100000N12qe,where N1=1−x/l;N2=x/l;N3=1−3x2/l2+2x3/l3;N4=x−2x2/l+x3/l2; N5=3x2/l2−2x3/l3;N6=−x2/l+x3/l2;N7=1−3x2/l2+2x3/l3;N8=2x2/l−x3/l2−x; N9=3x2/l2−2x3/l3;N10=x2/l−x3/l2;N11=1−x/l;andN12=x/l.For the convenience of subsequent calculation, orderu=Nuqu,v=Nvqv,w=Nwqw,θz=Nθzqθ. Substituting (22) into (9) by Hamiltonian principle [18, 19], ignoring infinitesimal quantity, the mass matrix Me of the UDS is obtained as (23).


### 3.2. Stiffness Matrix

The linear strain energy of the UDS is shown in (15). In the same way as solving the unit mass matrix, substituting (22) into (15), the linear stiffness matrix of the UDS is shown in (24). Substituting (22) into equations (16) to (19), the nonlinear stiffness of the UDS is shown in equations (25) to (28). The stiffness matrix of the UDS, combined with the linear and nonlinear ones of it, is shown in (29).(23)Me=∫0lAρNuTNu+NvTNv+NwTdz+ρIz∫0lNθzTNθzdz+ρIxy∫0lNw′TNw′+Nv′TNv′dz,(24)KL=∫0lλEANu′TNu′+GIzNθz′TNθz′+λEIyNu′′TNu′′+λEIxNv′′TNv′′dz,(25)KN1=∫0l3EλA2Nu′TNu′qeNu′dz,(26)KN2=∫0lEλIx12Nu′TNθz′qeNθz′+Nθz′TNu′qeNθz′dz,(27)KN3=∫0LEλA12Nu′TNw′qeNw′+Nw′TNu′qeNw′+12Nu′TNv′qeNv′+Nv′TNu′qeNv′dz+∫0lEλIz32Nu′TNv′′qeNv′+3Nv′′TNu′qeNv′′+3Nw′′TNu′qeNw′′dz+∫0l3EλIy2Nu′TNw′′qeNw′′dz,(28)KN4=λE−GIz2∫0lNw′′TNθz′qeNv′+Nv'TNθz′qeNw′′+Nθz′TNv′qeNw′′−Nθz′TNw′qeNv′′−Nw′TNθz'qeNv′′−Nv′′TNθz′qeNw′dz,(29)Ke=KL+KN1+KN2+KN3+KN4.


### 3.3. Damping Matrix

According to the Rayleigh damping [20], a unit damping matrix [21] is built by combing the mass matrix and linear stiffness matrix of the UDS with a certain proportion coefficient.(30)Ce=αMe+βKL,where α=2ξωaωb/ωa+ωbandβ=2ξ/ωa+ωb.The damping matrix changes with the natural frequency of the drill string. The frequency band that is easy to induce vibration can be obtained according to the kinetic characteristics of the drill string structure and the excitation of force loading.


### 3.4. Force Loading Matrix

There is mainly gravity (Fg), axial thrust (Fa), perturbed force (Fr), impact force (Fc), frictional force (Fm), reaction torque (Mo), and torque (Mt) on the drill string. According to the displacement matrix shown in (21), the directions of gravity and axial thrust are along the axis of the drill string, and their equivalent nodal forces are as equations (31) and (32). The torque rotates along the axial direction of the UDS, and its equivalent nodal force is shown in (33). If the lateral displacement is greater than the distance between the drill string and coal rock, it is judged that an impact collision occurs. The impact forces in the x-axis and y-axis are, respectively, as equations (34) and (35) [22], and their equivalent nodal force is shown in (36). The frictional force between the drill string and coal rock is shown in (37), and its equivalent nodal force is shown in (38) [23].(31)Fg=Algρ200000Algρ200000,(32)Fa=000000Fa00000,(33)Mt=000000000Mt00,(34)Fcx=−kvt−Ωsgnvt,vt≥Ω,0,else,(35)Fcy=−kwt−Ωsgnwt,wt≥Ω,0,else,(36)Fc=0FcyiFcxi0−dFcyi2dFcxi20FcyjFcxj0dFcyj2−dFcxj2,(37)Fm=−μdk2vt2+wt2−Ωsgnϖ,(38)Fm=0FmyiFmxi0−dFmyi2dFmxi20FmyjFmxj0dFmyj2−dFmxj2.


### 3.5. Assembly of Unit Matrixes

According to the order of units, the global force loading matrix is assembled, as shown in equation (39).Since the local coordinate system of each unit is consistent with the global one, coordinate transformation is not required. Equations (23) and (29) are assembled according to the order of the units by using the direct stiffness method [24, 25]. In this paper, the number n of the units is 7 (8 nodes). The characteristics of the global mass and stiffness matrix of the drill string are obtained by MATLAB, as shown in Figure 4.(39)F=fj1fi1+fj2fi2+fj3fi3+fj4⋯finT,where f=Fg+Fc+Fa+Mt+Fm−Fr−Mo.Figure 4
Distribution diagram of global mass and stiffness matrix for nonzero units. (a) Global mass matrix; (b) Global stiffness matrix.
(a)(b)As shown in Figure4, the global mass matrix and stiffness matrix are of order 48. The blue dot is the nonzero units in the matrixes; and the corresponding position of the zero units are blank. The number of nonzero units in the global mass matrix is 196, and that in the global stiffness matrix is 534. The nonzero units in Figures 4(a) and 4(b) are concentrated in a strip area centered on the diagonal. The zero units account for the majority, showing banded and sparse characteristics. There is obvious symmetry in Figure 4(a), while asymmetric in Figure 4(b). Consequently, the global stiffness matrix is singular and irreversible.After the UDS matrix is assembled into a global matrix, the kinetic equation of the drill string is shown as(40)MQ¨+CQ˙+KQ=F.


## 3.1. Mass Matrix

Taking any unit on the drill string to analyse, the upper and lower nodes of it are represented byj and i, respectively. The unit coordinate systems j-xyz and i-xyz are established respectively, as shown in Figure 3.Figure 3
Generalized coordinates of the nodal displacement of UDS.There are 6-DOF for each node of the discrete model of the UDS. The displacement matrix of the unit is as follows.(21)qe=uiviwiθxiθyiθziujvjwjθxjθyjθzjT.According to the displacement shape function [17], each displacement matrix of the UDS is as follows:(22)uvwθz=N100000N2000000N3N40000N5N6000000N7N80000N9N10000000N1100000N12qe,where N1=1−x/l;N2=x/l;N3=1−3x2/l2+2x3/l3;N4=x−2x2/l+x3/l2; N5=3x2/l2−2x3/l3;N6=−x2/l+x3/l2;N7=1−3x2/l2+2x3/l3;N8=2x2/l−x3/l2−x; N9=3x2/l2−2x3/l3;N10=x2/l−x3/l2;N11=1−x/l;andN12=x/l.For the convenience of subsequent calculation, orderu=Nuqu,v=Nvqv,w=Nwqw,θz=Nθzqθ. Substituting (22) into (9) by Hamiltonian principle [18, 19], ignoring infinitesimal quantity, the mass matrix Me of the UDS is obtained as (23).


## 3.2. Stiffness Matrix

The linear strain energy of the UDS is shown in (15). In the same way as solving the unit mass matrix, substituting (22) into (15), the linear stiffness matrix of the UDS is shown in (24). Substituting (22) into equations (16) to (19), the nonlinear stiffness of the UDS is shown in equations (25) to (28). The stiffness matrix of the UDS, combined with the linear and nonlinear ones of it, is shown in (29).(23)Me=∫0lAρNuTNu+NvTNv+NwTdz+ρIz∫0lNθzTNθzdz+ρIxy∫0lNw′TNw′+Nv′TNv′dz,(24)KL=∫0lλEANu′TNu′+GIzNθz′TNθz′+λEIyNu′′TNu′′+λEIxNv′′TNv′′dz,(25)KN1=∫0l3EλA2Nu′TNu′qeNu′dz,(26)KN2=∫0lEλIx12Nu′TNθz′qeNθz′+Nθz′TNu′qeNθz′dz,(27)KN3=∫0LEλA12Nu′TNw′qeNw′+Nw′TNu′qeNw′+12Nu′TNv′qeNv′+Nv′TNu′qeNv′dz+∫0lEλIz32Nu′TNv′′qeNv′+3Nv′′TNu′qeNv′′+3Nw′′TNu′qeNw′′dz+∫0l3EλIy2Nu′TNw′′qeNw′′dz,(28)KN4=λE−GIz2∫0lNw′′TNθz′qeNv′+Nv'TNθz′qeNw′′+Nθz′TNv′qeNw′′−Nθz′TNw′qeNv′′−Nw′TNθz'qeNv′′−Nv′′TNθz′qeNw′dz,(29)Ke=KL+KN1+KN2+KN3+KN4.


## 3.3. Damping Matrix

According to the Rayleigh damping [20], a unit damping matrix [21] is built by combing the mass matrix and linear stiffness matrix of the UDS with a certain proportion coefficient.(30)Ce=αMe+βKL,where α=2ξωaωb/ωa+ωbandβ=2ξ/ωa+ωb.The damping matrix changes with the natural frequency of the drill string. The frequency band that is easy to induce vibration can be obtained according to the kinetic characteristics of the drill string structure and the excitation of force loading.


## 3.4. Force Loading Matrix

There is mainly gravity (Fg), axial thrust (Fa), perturbed force (Fr), impact force (Fc), frictional force (Fm), reaction torque (Mo), and torque (Mt) on the drill string. According to the displacement matrix shown in (21), the directions of gravity and axial thrust are along the axis of the drill string, and their equivalent nodal forces are as equations (31) and (32). The torque rotates along the axial direction of the UDS, and its equivalent nodal force is shown in (33). If the lateral displacement is greater than the distance between the drill string and coal rock, it is judged that an impact collision occurs. The impact forces in the x-axis and y-axis are, respectively, as equations (34) and (35) [22], and their equivalent nodal force is shown in (36). The frictional force between the drill string and coal rock is shown in (37), and its equivalent nodal force is shown in (38) [23].(31)Fg=Algρ200000Algρ200000,(32)Fa=000000Fa00000,(33)Mt=000000000Mt00,(34)Fcx=−kvt−Ωsgnvt,vt≥Ω,0,else,(35)Fcy=−kwt−Ωsgnwt,wt≥Ω,0,else,(36)Fc=0FcyiFcxi0−dFcyi2dFcxi20FcyjFcxj0dFcyj2−dFcxj2,(37)Fm=−μdk2vt2+wt2−Ωsgnϖ,(38)Fm=0FmyiFmxi0−dFmyi2dFmxi20FmyjFmxj0dFmyj2−dFmxj2.


## 3.5. Assembly of Unit Matrixes

According to the order of units, the global force loading matrix is assembled, as shown in equation (39).Since the local coordinate system of each unit is consistent with the global one, coordinate transformation is not required. Equations (23) and (29) are assembled according to the order of the units by using the direct stiffness method [24, 25]. In this paper, the number n of the units is 7 (8 nodes). The characteristics of the global mass and stiffness matrix of the drill string are obtained by MATLAB, as shown in Figure 4.(39)F=fj1fi1+fj2fi2+fj3fi3+fj4⋯finT,where f=Fg+Fc+Fa+Mt+Fm−Fr−Mo.Figure 4
Distribution diagram of global mass and stiffness matrix for nonzero units. (a) Global mass matrix; (b) Global stiffness matrix.
(a)(b)As shown in Figure4, the global mass matrix and stiffness matrix are of order 48. The blue dot is the nonzero units in the matrixes; and the corresponding position of the zero units are blank. The number of nonzero units in the global mass matrix is 196, and that in the global stiffness matrix is 534. The nonzero units in Figures 4(a) and 4(b) are concentrated in a strip area centered on the diagonal. The zero units account for the majority, showing banded and sparse characteristics. There is obvious symmetry in Figure 4(a), while asymmetric in Figure 4(b). Consequently, the global stiffness matrix is singular and irreversible.After the UDS matrix is assembled into a global matrix, the kinetic equation of the drill string is shown as(40)MQ¨+CQ˙+KQ=F.


## 4. Boundary Condition

Working conditions of the drill string are as follows: (1) the axial thrust and torque are applied to the tail of the drill string by the roof bolter; (2) the drill string drill along its axis. Consequently, the boundary conditions of the drill string system are as follows: (1) all DOF of the first node of the drill string are constrained, in which the node displacements of the element are all zero. (2) the lateral displacements of the 8th node is constrained, namely, in which the lateral displacements in x and y axes are eliminated. The boundary conditions include global mass, stiffness, damping, and force loading and can be expressed as(41)Mr×1=M1×c=0Kr×1=K1×c=0Cr×1=C1×c=0F1×1=0Mr×2=M2×c=0Kr×2=K2×c=0Cr×2=C2×c=0F2×1=0⋯⋯⋯⋯Mr×6=M6×c=0Kr×6=K6×c=0Cr×6=C6×c=0F6×1=0Mr×44=M44×c=0Kr×44=K44×c=0Cr×44=C44×c=0F44×1=0Mr×45=M45×c=0Kr×45=K45×c=0Cr×45=C45×c=0F45×1=0.There are 7 units and 8 nodes, and 6-DOF in each node. (41) is a 6 × 8 order matrix, in which Mr×1=M1×c=0 indicates that the elements in row r, column 1 of Mr×1 and in row 1, column c of M1×c are all zero, and so on for the rest.Equation (41) is applied to (40) to eliminate the DOF of the displacement of the corresponding node so as to eliminate the singularity of the global stiffness matrix of the drill string.


## 5. Simulated Analysis for Drilling Mudstone

### 5.1. Modal Analysis

This paper takes mudstone as the drilling object. The geotechnical mechanics parameters [26] are as follows: ρ = 2600 kg/m3, R = 38 MPa, Rm = 0.34 MPa, E = 12 GPa, μ = 0.25, Fa = 6000 N, Mt = 130 N·m, Ω = 0 m, k = 109N·m−1μ = 0.25, and ϖ = 220 r/min. The parameters of the drill string are shown in Table 1.Table 1
Parameters of drill string.
Parameters1st unit2nd unit3rd unit4th unit5th unit6th unit7th unit8th unitElasticity modulus (GPa)210210210210210210210210Density (kg/m−3)78507850785078507850785078507850Length (m)0.0860.0860.0860.0860.0860.0860.0860.035Outer diameter (m)0.0380.0380.0380.0380.0380.0380.0380.032Inner diameter (m)0.0340.0340.0340.0340.0340.0340.0340.006A drilling model is constructed according to the above drill string, mudstone parameters, and boundary conditions. As shown in Figure5, point “RP” is fixed, and the drilling direction is along the shaft of the drill string. A drilling simulation is implemented by ABAQUS [26], and then the perturbed force and reaction torque of mudstone on the drill string are obtained.Figure 5
Drilling model.The first six modal shapes of the drill string are obtained by ABAQUS used for the modal analysis of it, as shown in Figure6. The first to third order mode shapes of the drill string are oscillating; and the fourth to sixth ones are second-order bending. The natural frequencies of the first six modal shapes are 0.051, 0.224, 0.314, 0.715, 0.870, and 1.094 Hz, respectively. The order of the mode shape is positively correlated with the corresponding natural frequency. The higher the natural frequency, the smaller the vibration period and the faster the attenuation, and vice versa. It means that the low-order mode shape plays a decisive role in the vibration characteristics of the drill string.Figure 6
First six mode shapes of the drill string. (a) 1st order; (b) 2nd order; (c) 3rd order; (d) 4th order; (e) 5th order; and (f) 6th order.
(a)(b)(c)(d)(e)(f)Consequently, the natural frequencies of the first and second mode shapes are employed, and their corresponding angular frequencies are as follows:ωa = 0.320 rad/s and ωb = 0.224 rad/s. Taking the ξ = 0.0048, substituting the parameters above into (30), the scale coefficients are obtained as follows: α = 1.26 × 10−3 and β = 1.76 × 10−2.


### 5.2. Interpolation of Perturbed Force and Reaction Torque

To substitute the perturbed force and reaction torque into (39), it is necessary to interpolate them. The perturbed force and reaction torque of the mudstone on the drill string are imported into the “ppval” function of MATLAB for cubic interpolation. As shown in Figure 7, “∗” indicates data points; the sampling interval is 0.001 s; and the sampling time t is set to 3 s. It can be seen from Figure 7(a) that the disturbing force changes periodically. The perturbed force on the drill string is zero while t at [0.49, 0.87], [1.10, 1.38], [1.59, 1.83], [1.99, 2.16], [2.32, 2.48], [2.71, 2.82], and the drill string rebounds. As can be seen from Figure 7(b), when t is at [0.63, 1.05], [1.60, 1.87], [2.09, 2.21], and [2.39, 2.87], the reaction torque on the drill string is zero. At this time, the drill string is stuck and stops suddenly, and the actual revolving speed is zero.Figure 7
Perturbed force and reaction torque on the drill pipe: (a) perturbed force and (b) reaction torque.
(a)(b)


### 5.3. Displacement of Unit Vibration

Employing Simulink [27] in MATLAB, the responses to the coupled vibration on the drill pipe are output to 7 oscilloscopes, respectively, as shown in Figure 8. Substitute (41) into (40), the modified global matrix, including mass, stiffness, damping, and force load, are substituted into the “finite element model.” The calculation process can be described as follows.Input: FOutput: Q(1) Substituting to (33), (36), (38), (31), and the interpolation results of Fr and Mo into (39), assemble global force loading matrix F;(2) Assembling (23) into M, solve M−1;(3) According to (40), solve the equation Q¨=M−1F−CQ˙−KQ by the “ode23tb” decoder, and then obtain Q¨, Q˙, and Q.Figure 8
Simulink model.The coupled vibration displacement responses to the drilling of roof bolter are solved by using the “ode23tb” decoder in Simulink, as shown in Figure9.Figure 9
Vibration displacement responses to the drilling of roof bolter. (a) 2nd unit; (b) 3rd unit; (c) 4th unit; (d) 5th unit; (e) 6th unit; (f) 7th unit; and (g) 8th unit.
(a)(b)(c)(d)(e)(f)(g)As shown in Figure9, un, vn, wn, and θxn are the responses to the longitudinal, x-axis, y-axis, and torsional vibration displacement of the UDS, respectively. The positive or negative values of the vibration displacement curves only represent their directions. The variation intervals of u8 to u2 shown in Figures 9(a) to 9(f) are [−2.76 × 10−3, 5.21 × 10−3], [−6.57 × 10−4, 3.03 × 10−4], [0, 6.57 × 10−4], [0, 5.86 × 10−4], [0, 4.86 × 10−4], [0, 3.54 × 10−4], and [0, 1.92 × 10−4], respectively. Their vibration responses curves are similar, and the amplitudes of them decrease gradually. Compared with the bit, the longitudinal vibration displacement of the tail of the drill string is decreased by 93.2%. There is zero in the curves of the longitudinal vibration displacements during this period, which is due to the bouncing movement of the drill string that is consistent with the rebound time of the drill string in Figure 7(a). Compared with Figures 9(a) to 9(f), the amplitude of u8 increased significantly, which shows that during process of fragmentizing coal rock, the longitudinal vibration is first transmitted to the drill bit, and then attenuated to each unit of the drill string in turn. It means that the perturbed force is the main cause of inducing the longitudinal vibration.Thev8∼v2 and w8∼w2 curves in Figures 9(a) to 9(f) are shown as mirror symmetry, but the amplitudes of them is not consistent, which shows that the characteristics of the impact forces caused by the bending vibration are similar, but their magnitudes are different. In the time intervals [0.61, 0.99], [1.28, 1.396], [1.62, 1.88], and [2.39, 2.87], the displacement caused by it is zero, which means that the contact collision between the drill string and the mudstone does not occur continuously.From Figures9(a) to 9(g), the trends of θx8 to θx2 are completely consistent. Their amplitudes decrease by 79.1%, 16.3%, 19.9%, 25.2%, 33.8%, and 49.9% in order, which shows a gradual decreasing trend. Compared with the bit, the torsional vibration displacement of the tail of the drill string decreased by 84.7%. By observing Figures 9(b) and 9(a) to 9(g), it is found that the curve of reaction torque is linearly related to the curve of the torsional vibration displacement of the drill string. Their vibration displacements are intermittently zero, indicating that the drill strings are stuck, which is consistent with the sticking time in Figure 7(b). This means that the uneven change of the reaction torque is the main cause of the torsional vibration.


## 5.1. Modal Analysis

This paper takes mudstone as the drilling object. The geotechnical mechanics parameters [26] are as follows: ρ = 2600 kg/m3, R = 38 MPa, Rm = 0.34 MPa, E = 12 GPa, μ = 0.25, Fa = 6000 N, Mt = 130 N·m, Ω = 0 m, k = 109N·m−1μ = 0.25, and ϖ = 220 r/min. The parameters of the drill string are shown in Table 1.Table 1
Parameters of drill string.
Parameters1st unit2nd unit3rd unit4th unit5th unit6th unit7th unit8th unitElasticity modulus (GPa)210210210210210210210210Density (kg/m−3)78507850785078507850785078507850Length (m)0.0860.0860.0860.0860.0860.0860.0860.035Outer diameter (m)0.0380.0380.0380.0380.0380.0380.0380.032Inner diameter (m)0.0340.0340.0340.0340.0340.0340.0340.006A drilling model is constructed according to the above drill string, mudstone parameters, and boundary conditions. As shown in Figure5, point “RP” is fixed, and the drilling direction is along the shaft of the drill string. A drilling simulation is implemented by ABAQUS [26], and then the perturbed force and reaction torque of mudstone on the drill string are obtained.Figure 5
Drilling model.The first six modal shapes of the drill string are obtained by ABAQUS used for the modal analysis of it, as shown in Figure6. The first to third order mode shapes of the drill string are oscillating; and the fourth to sixth ones are second-order bending. The natural frequencies of the first six modal shapes are 0.051, 0.224, 0.314, 0.715, 0.870, and 1.094 Hz, respectively. The order of the mode shape is positively correlated with the corresponding natural frequency. The higher the natural frequency, the smaller the vibration period and the faster the attenuation, and vice versa. It means that the low-order mode shape plays a decisive role in the vibration characteristics of the drill string.Figure 6
First six mode shapes of the drill string. (a) 1st order; (b) 2nd order; (c) 3rd order; (d) 4th order; (e) 5th order; and (f) 6th order.
(a)(b)(c)(d)(e)(f)Consequently, the natural frequencies of the first and second mode shapes are employed, and their corresponding angular frequencies are as follows:ωa = 0.320 rad/s and ωb = 0.224 rad/s. Taking the ξ = 0.0048, substituting the parameters above into (30), the scale coefficients are obtained as follows: α = 1.26 × 10−3 and β = 1.76 × 10−2.


## 5.2. Interpolation of Perturbed Force and Reaction Torque

To substitute the perturbed force and reaction torque into (39), it is necessary to interpolate them. The perturbed force and reaction torque of the mudstone on the drill string are imported into the “ppval” function of MATLAB for cubic interpolation. As shown in Figure 7, “∗” indicates data points; the sampling interval is 0.001 s; and the sampling time t is set to 3 s. It can be seen from Figure 7(a) that the disturbing force changes periodically. The perturbed force on the drill string is zero while t at [0.49, 0.87], [1.10, 1.38], [1.59, 1.83], [1.99, 2.16], [2.32, 2.48], [2.71, 2.82], and the drill string rebounds. As can be seen from Figure 7(b), when t is at [0.63, 1.05], [1.60, 1.87], [2.09, 2.21], and [2.39, 2.87], the reaction torque on the drill string is zero. At this time, the drill string is stuck and stops suddenly, and the actual revolving speed is zero.Figure 7
Perturbed force and reaction torque on the drill pipe: (a) perturbed force and (b) reaction torque.
(a)(b)


## 5.3. Displacement of Unit Vibration

Employing Simulink [27] in MATLAB, the responses to the coupled vibration on the drill pipe are output to 7 oscilloscopes, respectively, as shown in Figure 8. Substitute (41) into (40), the modified global matrix, including mass, stiffness, damping, and force load, are substituted into the “finite element model.” The calculation process can be described as follows.Input: FOutput: Q(1) Substituting to (33), (36), (38), (31), and the interpolation results of Fr and Mo into (39), assemble global force loading matrix F;(2) Assembling (23) into M, solve M−1;(3) According to (40), solve the equation Q¨=M−1F−CQ˙−KQ by the “ode23tb” decoder, and then obtain Q¨, Q˙, and Q.Figure 8
Simulink model.The coupled vibration displacement responses to the drilling of roof bolter are solved by using the “ode23tb” decoder in Simulink, as shown in Figure9.Figure 9
Vibration displacement responses to the drilling of roof bolter. (a) 2nd unit; (b) 3rd unit; (c) 4th unit; (d) 5th unit; (e) 6th unit; (f) 7th unit; and (g) 8th unit.
(a)(b)(c)(d)(e)(f)(g)As shown in Figure9, un, vn, wn, and θxn are the responses to the longitudinal, x-axis, y-axis, and torsional vibration displacement of the UDS, respectively. The positive or negative values of the vibration displacement curves only represent their directions. The variation intervals of u8 to u2 shown in Figures 9(a) to 9(f) are [−2.76 × 10−3, 5.21 × 10−3], [−6.57 × 10−4, 3.03 × 10−4], [0, 6.57 × 10−4], [0, 5.86 × 10−4], [0, 4.86 × 10−4], [0, 3.54 × 10−4], and [0, 1.92 × 10−4], respectively. Their vibration responses curves are similar, and the amplitudes of them decrease gradually. Compared with the bit, the longitudinal vibration displacement of the tail of the drill string is decreased by 93.2%. There is zero in the curves of the longitudinal vibration displacements during this period, which is due to the bouncing movement of the drill string that is consistent with the rebound time of the drill string in Figure 7(a). Compared with Figures 9(a) to 9(f), the amplitude of u8 increased significantly, which shows that during process of fragmentizing coal rock, the longitudinal vibration is first transmitted to the drill bit, and then attenuated to each unit of the drill string in turn. It means that the perturbed force is the main cause of inducing the longitudinal vibration.Thev8∼v2 and w8∼w2 curves in Figures 9(a) to 9(f) are shown as mirror symmetry, but the amplitudes of them is not consistent, which shows that the characteristics of the impact forces caused by the bending vibration are similar, but their magnitudes are different. In the time intervals [0.61, 0.99], [1.28, 1.396], [1.62, 1.88], and [2.39, 2.87], the displacement caused by it is zero, which means that the contact collision between the drill string and the mudstone does not occur continuously.From Figures9(a) to 9(g), the trends of θx8 to θx2 are completely consistent. Their amplitudes decrease by 79.1%, 16.3%, 19.9%, 25.2%, 33.8%, and 49.9% in order, which shows a gradual decreasing trend. Compared with the bit, the torsional vibration displacement of the tail of the drill string decreased by 84.7%. By observing Figures 9(b) and 9(a) to 9(g), it is found that the curve of reaction torque is linearly related to the curve of the torsional vibration displacement of the drill string. Their vibration displacements are intermittently zero, indicating that the drill strings are stuck, which is consistent with the sticking time in Figure 7(b). This means that the uneven change of the reaction torque is the main cause of the torsional vibration.


## 6. Conclusions

(1)
Based on the energy method and finite element method, combined with the boundary conditions of the cantilever beam, and drilling along the axis of the roof bolter, the reduced-order kinetic models of the bending, longitudinal, and torsional coupled vibrations of the drill string are constructed. The coupled vibration displacement responses of each unit of the drill pipe are obtained, in which the coupled vibration mechanism of the drilling of a roof bolter for a mine support system has been revealed.(2)
According to the simulated analysis of mudstone drilling, the order of causing drill string displacement from large to small is as follows: longitudinal vibration, torsional vibration, and bending vibration. The curve of reaction torque is linearly related to the curve of the torsional vibration displacement of the drill string. The displacements of the longitudinal vibration and torsional vibration decrease gradually from the bit to the tail of the drill string. In the instance of this paper, the longitudinal and torsional vibration displacements of the tail are decreased by 93.2% and 84.7%, respectively, compared with those of the tip. The dynamic characteristics of the tail of the drill string under coupled vibration can be predicted according to the longitudinal and torsional vibration displacement of the bit.(3)
The bending vibration displacement curves along thex and y-axes are shown as mirror symmetry, but the amplitudes of them are not consistent. Moreover, they are intermittently zero, indicating that the contact collision between the drill string and the mudstone does not occur continuously. The perturbed force and reaction torque are the main causes of inducing the coupled vibration. To reduce the vibration amplitude of drill string and avoid tripping and sticking, the axial thrust and torque of drill string can be reduced, which can reduce its disturbing force and reaction torque.The 6-DOF coupled vibration mechanism analysed in this paper can be extended to various drill strings and different kinds of coal rock, which provides a theoretical basis for analysing the vibration characteristics of the drilling of roof bolters for mine support systems and developing their absorbers.



---

*Source: 4279959-2022-05-17.xml*