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.
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 and 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 Figure 1, the drill string system is discretized into 1 to units, and each unit has 6-DOF. The displacement of longitudinal vibration is ; that of bending vibration are and ; the rotation angle around , , and axis is , , and , respectively.

The translational kinetic energy of the unit of drill string (UDS) is expressed as [12]
According to Figure 1, the translational velocity of the unit is
As shown in Figure 2, the rotation velocity of the unit is expressed by Euler angle. Firstly, rotating the global coordinate system around the -axis by an angle to get the coordinate system ; secondly, rotating it again around the -axis by an angle to get the coordinate system ; finally, rotating it once more around the -axis by an angle to get the coordinate system .

(a)

(b)

(c)
Using Euler kinematics equation [13], the rotational velocity is shown as,
The rotational kinetic energy of the UDS is expressed as [12]
Due to the rotation angle of the drill string restricted by coal rock, there is an assumption as
According to simultaneous (3) and (8), the total kinetic energy of the UDS is shown as
2.3. Strain Energy Equation
The strain energy of the UDS is described as follows [14]:
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
Substituting (11) into (10), the strain energy of the UDS iswhere .
The strains of the UDS are as follows [15]:
According to the assumption of Euler–Bernoulli beam [16], compared with bending deformation, shear deformation can be ignored, thus,
Substituting (14) into (13), and then substituting the result into (12), the linear strain energy of the UDS is as follows:
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:
According to simultaneous (16) and (19), the total strain energy of the UDS is shown as
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 by j and i, respectively. The unit coordinate systems j-xyz and i-xyz are established respectively, as shown in Figure 3.

There are 6-DOF for each node of the discrete model of the UDS. The displacement matrix of the unit is as follows.
According to the displacement shape function [17], each displacement matrix of the UDS is as follows:where ; ; .
For the convenience of subsequent calculation, order . Substituting (22) into (9) by Hamiltonian principle [18, 19], ignoring infinitesimal quantity, the mass matrix 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).
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.where .
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 (), axial thrust (), perturbed force (), impact force (), frictional force (), reaction torque (), and torque () 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 -axis and -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].
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 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.where .

(a)

(b)
As shown in Figure 4, 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
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 and axes are eliminated. The boundary conditions include global mass, stiffness, damping, and force loading and can be expressed as
There are 7 units and 8 nodes, and 6-DOF in each node. (41) is a 6 × 8 order matrix, in which indicates that the elements in row , column 1 of and in row 1, column of 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, = 6000 N, = 130 N·m, = 0 m, = 109N·m−1 = 0.25, and = 220 r/min. The parameters of the drill string are shown in Table 1.
A drilling model is constructed according to the above drill string, mudstone parameters, and boundary conditions. As shown in Figure 5, 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.

The first six modal shapes of the drill string are obtained by ABAQUS used for the modal analysis of it, as shown in Figure 6. 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.

(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: = 0.320 rad/s and = 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 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 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 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.

(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: Output: (1) Substituting to (33), (36), (38), (31), and the interpolation results of and into (39), assemble global force loading matrix ; (2) Assembling (23) into , solve ; (3) According to (40), solve the equation by the “ode23tb” decoder, and then obtain , , and .

The coupled vibration displacement responses to the drilling of roof bolter are solved by using the “ode23tb” decoder in Simulink, as shown in Figure 9.

(a)

(b)

(c)

(d)

(e)

(f)

(g)
As shown in Figure 9, , , , and are the responses to the longitudinal, -axis, -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 to 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 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.
The and 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 Figures 9(a) to 9(g), the trends of to 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 the x 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.
Nomenclatures
: | Elastic coefficient of impact force (N·m−1) |
: | Displacement matrix |
: | Translational velocity (m/s) |
: | Drill string outer diameter (m) |
: | Scale factor of gravity (N/kg) |
: | Rotational velocity (m/s) |
: | Total kinetic energy |
: | Translational kinetic |
: | Rotational kinetic |
: | Unit length (m) |
: | Cross sectional area (m2) |
: | Inertia tensor |
: | Moment of inertia of -axis |
: | Moment of inertia of -axis |
: | Moment of inertia of -axis |
: | Global stiffness matrix |
: | Global mass matrix |
: | Global damping matrix |
: | Generalized acceleration |
: | Generalized velocity |
: | Generalized displacement |
: | Global force loading matrix |
: | Elastic matrix |
: | Unit volume |
R: | Compressive strength (MPa) |
Rm: | Tensile strength (MPa) |
E: | Elastic modulus (GPa) |
: | Nonlinear strain energy caused by longitudinal vibration |
: | Nonlinear strain energy caused by coupled longitudinal and torsional vibration |
: | Shape function |
: | Nonlinear strain energy caused by longitudinal and bending coupled vibration |
: | Nonlinear strain energy caused by coupled bending and torsional vibration |
: | Strain energy |
: | Linear strain energy |
: | Linear stiffness matrix |
: | Longitudinal vibration stiffness matrix |
[]: | Longitudinal and torsional coupled vibration stiffness matrix |
[]: | Longitudinal and bending coupled vibration stiffness matrix |
[]: | Transverse and bending coupling vibration stiffness matrix |
: | Impact force in x-axis direction (N) |
: | Impact force in y-axis direction (N) |
: | Frictional force (N) |
: | Frictional force in x-axis direction (N) |
: | Frictional force in y-axis direction (N) |
: | Perturbed force (N) |
: | Reaction torque (N·m) |
: | Torque (N·m) |
: | Material density (kg/m3) |
: | Gap between drill string and coal rock (m) |
: | Poisson’s ratio |
: | Friction coefficient |
: | Drill string speed (r/min) |
: | Stress matrix |
: | Strain matrix |
: | Damping ratio |
: | -order circular frequency (rad/s) |
: | -order circular frequency (rad/s) |
: | Scale factor in Rayleigh damping |
: | Scale factor in Rayleigh damping. |
: | The first derivative with respect to |
: | The second derivative with respect to |
: | The first derivative with respect to t |
: | The second derivative with respect to t |
: | Unit number. |
: | Number of matrix rows, |
: | Number of matrix columns, |
: | Lower node |
: | Upper node |
: | Rotational velocity of x-axis |
: | Rotational velocity of y-axis |
: | Rotational velocity of z-axis. |
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
All authors contributed equally to this article.
Acknowledgments
This work was supported in part by the Special Fund for Collaborative Innovation of Anhui Polytechnic University and Jiujiang District under grant no. 2021CYXTB3 and by the University Synergy Innovation Program of Anhui Province of China (grant no. GXXT-2019-048).