Modeling of chatter vibrations in gun drilling process

Sergey Voronov1 , Ivan Pleshcheev2

1, 2Applied Mechanics Department, Bauman Moscow State Technical University, Moscow, Russia

1Corresponding author

Vibroengineering PROCEDIA, Vol. 38, 2021, p. 19-25. https://doi.org/10.21595/vp.2021.22084
Received 19 May 2021; received in revised form 3 June 2021; accepted 10 June 2021; published 28 June 2021

Copyright © 2021 Sergey Voronov, et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Table of Contents Download PDF References
Cite this article
Views 22
Reads 8
Downloads 118
CrossRef Citations 0
Abstract.

The dynamics of the gun drilling process is analyzed in this paper. The tool shank is modeled as long straight beam vibrating in transverse direction under action of cutting forces. Axial force component is expressed as proportional to cutting thickness, which is determined as nonlinear function of beam transverse deflection with time delay. Nonlinear equations of motion of the drilling shank are derived. The stability diagram of the system dynamics was determined. The bifurcation analysis of nonlinear differential delay equations by means of multiple scale method was performed. The obtained results were verified by numerical integration of nonlinear equations. The influence of cutting conditions on system stability and chatter amplitude was observed.

Keywords: gundrilling, chatter, nonlinear dynamics.

1. Introduction

Long straight holes are usually produced by means of a special drilling tool. This efficient process is widely used in the automotive industry to drill deep holes in cylinder heads, crankshafts, fuel pump housings, turbine blades and etc. Gundrilling is the mostly used method of deep small hole machining. Gun drill has asymmetrical single edge tool design with long straight tube of asymmetrical cross-section with a typical reachable diameter range of 0.5 mm up to 40 mm and length-to-diameter- ratios up to L/D = 400 (in special applications even L/D = 900 [1]). The method is widely applied in machining small deep holes as it provides a good straightness and high quality of machined surface due to its self-guiding action [2]. Optimal drill performance in gundrilling is achieved when the combination of the cutting speed, feed rate, tool geometry, carbide grade, and coolant parameters are selected properly depending upon the work material, deep-hole tool machine conditions, and the quality requirements to the drilled holes. Due to low flexural stiffness of gun drill shank lateral vibrations of high magnitude could be excited during the machining. Excessive vibrations are detrimental to finish surface quality and may damage the tool. Therefore, it is important to predict in advance regimes with chatter vibrations. The regenerative mechanism is the main source of chatter vibrations. And it requires that time-delayed terms in model equations should be taken into account. The same mechanism emerges not only in drilling [3-5], but in milling [6], boring etc. The comprehensive review of present state of deep hole drilling modeling was given in [1]. Most authors modeled drill shank using the reduced single degree of freedom system. In this paper gun drill is considered as flexible continuous beam loaded with eccentrically applied cutting force. The new approach allows considering the influence of lateral vibrations on the dynamics of the gun drilling system. The multiple scale method is applied for nonlinear vibrations analysis. Stability diagram was constructed and bifurcation diagrams were obtained by multi-scale expansion. The nonlinear behavior of system in vicinity of stability borders was analyzed by using numerical integration of nonlinear equation.

2. Model of gun drilling system

Gun drill shank is modeled as long slender beam using Euler-Bernoulli beam theory. The beam axial line is supposed to be straight line in undeformed state. Schematic of gun drilling process is presented on Fig. 1. Usually gun drill rests upon intermediate supports, but we suppose that clearance in these supports is large enough, so that these supports are not taken into consideration. Therefore, for the sake of simplicity we assume, that gun drilling shank as simply supported beam.

Fig. 1. Scheme of gun drilling process.

 Scheme of gun drilling process.

Gun drill tool has asymmetric cutting edge design, components of cutting force are shown in Fig. 2(a). In-plane component of cutting force (Fpl) is balanced by guiding pads, while tangential component of force leads to torque loading upon gun drill.

Axial component of cutting force is applied with eccentricity – e. Therefore, beam is loaded with axial force and bending moment. The mathematical model is presented in Fig. 2(b).

Fig. 2. Loading conditions of gun drill cutting edge and drawing of mathematical model of gun drilling system

 Loading conditions of gun drill cutting edge and drawing  of mathematical model of gun drilling system

a) Loading conditions of gun drill cutting edge

 Loading conditions of gun drill cutting edge and drawing  of mathematical model of gun drilling system

b) Drawing of mathematical model of gun drilling system

Nonlinear partial differential equation of beam vibration in plane of minimal rigidity were derived in [7] as:

(1)
m 0 2 w t 2 + E I 4 w s 4 + F a 2 w s 2 = - E I s 3 w s 3 w s 2 + 2 w s 2 2 w s
                + 1 2 m 0 w s s L 2 t 2 0 s w s 2 d s d s - 3 2 F a w s 2 2 w s 2 + M b δ ' s - L ,

where w – lateral displacement of beam axis, t – time, s – axial coordinate, m0 – mass of unit length, EI – flexural stiffness of the beam, Fa – axial force acting at the beam end, Mb – bending moment, due to cutting force eccentricity, δ – Dirac function, L – length of the beam. Deflection in other plane and torsion are ignored.

Chatter vibrations occur due to regenerative mechanism as the vibrating tool cutting edge interacts with the surface left by edge at its previous pass [8, 9]. To model this phenomenon we assume, that cutting force is proportional to the uncut chip thickness h, which is dependent upon axial position of the cutting edge at instant moment and time delayed by period of workpiece rotation T:

(2)
F a = - K c h ,         M b = - e y K c h ,                     h = a + u L , t - u L , t - T ,

where Kc – cutting coefficient, ey – eccentricity of the applied cutting force, h – chip thickness, a – feed, T – delay equal to work piece rotation period, u – axial displacement.

Axial displacement of the beam end is expressed nonlinearly through lateral displacement as:

(3)
u L , t = - 0 L 1 2 w s 2 d s .

Under such loading transverse vibrations arise and it leads to axial deflection of tool tip and therefore varying uncut chip thickness and forces [10]. Various models of tool vibrations under drilling were analyzed by authors in [8, 11, 12]. Chatter vibrations in these papers are explained due to axial or axial-torsional deformation of tool. However, for gun drill with high length-to-diameter ratio flexural vibrations are more significant and are analyzed in the paper.

The equations of flexural-torsional tool vibration under finite tool shank deflection were analyzed in [9]. It was shown that chatter vibration occurs in vicinity of lower eigenfrequency with eigenmode corresponding to this frequency.

Equations of motion in nondimensional form were obtained [7], following nondimensional variables and parameters were introduced:

(4)
w ~ = w L ,             s ~ = s L e y ~ =   e y L ,             a ~ =   a L ,             t ~ = ω 0 t ,             T ~ = T ω 0 ,          
ω 0 = E I ρ A L 4 ,             K c ~ = K c L 3 E I .

Here and after we omit all tildes and use nondimensional variables only.

Next, we discretize system equations by means of Galerkin method using first flexural mode in plane of minimal flexural rigidity. We also add viscous damping to take into account dissipation of energy [7]:

(5)
w s , t = x t sin π s ,
x ¨ + 2 d x ˙ + p 0 2 x = - b 0 x 2 - x d 2 x - b 1 x ˙ 2 + x x ¨ x - b 2 x 3 + b 3 - b 4 x 2 - x d 2 ,

where:

(6)
p 0 2 = π 4 1 - K c a π 2 ,           b 0 = 1 4 K c π 4 ,             b 1 = π 3 2 π 3 - 3 8 π ,           b 2 = π 6 2 - 3 8 π 4 K c a ,
b 3 = 2 π K c e y a ,           b 4 =   π 3 K c e y 2 ,

x d = x ( t - T ) – time-delayed variable value, d – damping coefficient.

Further we will study the effect of nonlinear regenerative mechanism for chatter vibrations excitation, thus we simplify equations by neglecting terms corresponding the geometrical nonlinearity of the system, except for terms, which are associated with determining axial displacement and loading containing the delayed functions. This could be done due to the fact, that usually cutting coefficient is much higher compared to other terms and are highly responsible for chatter vibration excitation. The geometrical nonlinearity consideration leads to nonsignificant variation of chatter stability borders.

After such assumptions we get:

(7)
x ¨ + 2 d x ˙ + π 4 - π 2 K c ( a - π 2 4 x 2 + π 2 4 x d 2 ) x = 2 π K c e y a - π 2 4 x 2 + π 2 4 x d 2 .

Term 2πKceya represents constant component of bending moment, which defines the static deflection xst. To study dynamic displacement near the equilibrium position, we rewrite equation of motion in vicinity of bent equilibrium state for dynamic component of beam deflection x as follows:

(8)
x ¨ + 2 d x ˙ + p 0 2 x + k 1 x - x d + a 1 x 2 - x d 2 + a 2 x x - x d + a 3 x x 2 - x d 2 = 0 ,

where:

(9)
k 1 = π 4 2 K c x s t 2 + π 3 K c e y x s t ,     a 1 = π 3 2 K c e y + π 4 4 K c x s t ,     a 2 = π 4 2 K c x s t ,
a 3 = π 4 4 K c ,             x s t = 2 π K c e y a p 0 2 .

The Eq. (8) feature concludes in existence of time delay variables in linear and nonlinear terms, that requires stability analysis and investigation of stationary regimes of linearized equation, and self-vibrations analysis near the border of stability.

3. Numerical analysis

3.1. Stability of linearized equation

First, we perform linear stability analysis via D-decomposition technique in domain of system parameters (coordinates T-Kc). Linearized Eq. (8) reads as followed:

(10)
x ¨ + 2 d x ˙ + p 0 2 x + k 1 x - x d = 0 .

To find stability boundaries, we express unknown variable in exponential form. After separating real and imaginary term we get system of transcendental equations:

(11)
- ω 2 + p 0 2 + k 1 1 - cos ω T = 0 ,             2 d ω + k 1 sin ω T = 0 .

Solving these equations in respect to ω, T for different values of Kc we obtain classical periodically altering regions of stability and instability (Fig. 3). Values of parameters used in this paper are d= 0.0493, ey=6.17e-3, a=1.67 e-4.

The behavior of the system in vicinity of stability borders is of great interest. Bifurcation analysis of the dynamical system requires analysis of full nonlinear Eq. (8).

Fig. 3. Stability diagram of linearized system. A and B – regions where bifurcation analysis is performed

 Stability diagram of linearized system. A and B – regions where bifurcation analysis is performed

3.2. Nonlinear analysis

To investigate dynamic behavior after stability loss we perform bifurcation analysis of nonlinear differential delay equations by means of multiple scale method [13]. We introduce time, variable and time delay as:

(12)
d d t = T 0 + ε T 1 + ε 2 T 2 ,
(13)
x = ε x 1 T 0 , T 1 , T 2 + ε 2 x 2 T 0 , T 1 , T 2 + ε 3 x 3 T 0 , T 1 , T 2 ,
(14)
T = T c r + ε τ 1 + ε 2 τ 2 ,

where Tcr – critical delay corresponding to stability boundary.

The delayed variable is expanded as:

(15)
x i t - T = x i T 0 - T , T 1 - ε T , T 2 - ε 2 T = x i d + ε - T c r D 1 x i d - τ 1 D 0 x i d +
        + ε 2 - T c r D 2 x i d - τ 1 D 1 x i d - τ 2 D 0 x i d + 1 2 ( T c r 2 D 1,1 x i d + 2 T c r τ 1 D 0,1 x i d + τ 1 D 0,0 x i d ) ,

where Di=Ti, Di,j=2TiTj,  xid=xiT0-Tcr,T1,T2.

Substituting Eqs. (12)-(15) and equating the same terms of power ε, the following partial differential equations are obtained:

(16)
ε 1 :   D 0,0 x 1 + 2 d D 0 x 1 + k x 1 + k 1 x 1 - x 1 d = 0 ,
(17)
ε 2 :   D 0,0 x 2 + 2 d D 0 x 2 + k x 2 + k 1 x 2 - x 2 d =
        = - 2 D 0,1 x 1 - 2 d D 1 x 1 - k 1 T c r D 1 x 1 d + τ 1 D 0 x 1 d - a 1 x 1 2 - x 1 d 2 - a 2 x 1 x 1 - x 1 d ,
(18)
ε 3 :   D 0,0 x 3 + 2 d D 0 x 3 + k x 3 + k 1 x 3 - x 3 d =
          = - 2 D 0,1 x 2 - 2 D 0,2 x 1 - D 1,1 x 1 - 2 d D 1 x 2 - 2 d D 2 x 1 - k 1 e x p r 1 - e x p r 2 ,
(19)
e x p r 1 = T c r D 2 x 1 d + τ 1 D 1 x 1 d + τ 2 D 0 x 1 d - 1 2 T c r 2 D 1,1 x 1 d + 2 T c r τ 1 D 0,1 x 1 d + τ 1 D 0,0 x 1 d
        + + T c r D 1 x 2 d + τ 1 D 0 x 2 d ,
(20)
e x p r 2 = 2 a 1 x 1 x 2 - x 1 d x 2 d + x 1 d T c r D 1 x 1 d + τ 1 D 0 x 1 d +
        + a 2 x 1 x 2 - x 2 d + T c r D 1 x 1 d + τ 1 D 0 x 1 d + x 2 x 1 - x 1 d + a 3 x 1 x 1 2 - x 1 d 2 .

General solution of first equation is written as:

(21)
x 1 = A T 1 , T 2 cos ω T 0 + B T 1 , T 2 sin ω T 0 .

Quadratic terms do not introduce secular terms in second approximation, therefore τ1=0, and A, B are independent of time-scale T1. However, quadratic terms still contribute to third approximation equations. Second approximation has a form:

(22)
x 2 = A 2 T 1 cos 2 ω T 0 + B 2 T 1 sin 2 ω T 0 + C 2 T 1 ,

where A2, B2, C2 – coefficients, which are expressed in terms of A, B.

After elimination of secular terms in third approximation we get first order differential equations in respect to unknowns AT2, BT2.

It is convenient to introduce solution of equations in polar form:

(23)
A T 2 = C T 2 c o s φ T 2 ,   B T 2 = C T 2 s i n φ T 2 ,

where C, φ – amplitude and phase.

Eliminating φ from equations we get differential equation for amplitude C. To find limit cycle amplitude we assume, that T2C=0. As a result we will derive the algebraic equation for limit cycle amplitude C calculation. Resulting equation has a form:

(24)
f ω , T c r ,   d , p 0 , a 1 , a 2 , a 3 A 2 = g ω , T c r ,   d , p 0 , a 1 , a 2 , a 3 τ 2 .

Obtained functions f and g are quite cumbersome and are not presented here. All algebraic manipulations were performed in Matlab Symbolic toolbox.

The system behavior in the vicinity of the bifurcation point has been studied in order to analyze the nature of Andronov-Hopf bifurcation with respect to system parameters variation. Here we vary delay parameter, which corresponds to rotational velocity of gundrill.

Bifurcation diagram for two regions in the vicinity of the stability border. (Fig. 3, regions A, B) have been analyzed by means of asymptotic method. Amplitudes of limit cycle vs parameter T in Fig. 6 are presented. To verify results obtained by asymptotic method we perform numerical time integration of Eq. (8) using implicit scheme with iterations. Examples of resulting displacement time histories and its Fourier spectrum are shown in Figs. 4-5. Amplitudes of steady motion, obtained by numerical integration are compared with the results of asymptotic method. Bifurcation diagrams for two cases determined by two methods in Fig. 6 are presented.

Fig. 4. Time history for x and Fourier spectrum of steady motion. Parameters: Kc= 8000, T= 7.2433

 Time history for x and Fourier spectrum of steady motion. Parameters: Kc= 8000, T= 7.2433

a) Time history for x

 Time history for x and Fourier spectrum of steady motion. Parameters: Kc= 8000, T= 7.2433

b) Fourier spectrum of steady motion

Fig. 5. Time history for x and Fourier spectrum of steady motion. Parameters: Kc= 9000, T= 7.9261

Time history for x and Fourier spectrum of steady motion. Parameters: Kc= 9000, T= 7.9261

a) Time history for x

Time history for x and Fourier spectrum of steady motion. Parameters: Kc= 9000, T= 7.9261

b) Fourier spectrum of steady motion

Results obtained by asymptotic method are found to be in good agreement with results obtained by numerical integration.

Fig. 6. Bifurcation diagram determined by asymptotic analysis and numerical integration

 Bifurcation diagram determined by asymptotic analysis and numerical integration

a) Zone A, Kc= 8000, Tcr= 7.2258

 Bifurcation diagram determined by asymptotic analysis and numerical integration

b) Zone B, Kc= 9000, Tcr= 7.9061

4. Conclusions

In this paper nonlinear model of chatter vibrations of gundrilling tool is presented. Stability diagram of the linearized system was obtained. The relative cutting speed is the main critical parameter which effects on process stability. Post-critical behavior after stability loss was analyzed by means of method of multiple scales. The bifurcation diagrams of the dimensionless cutting speed (speed of detail rotation) influence on chatter vibration amplitude were determined. The results of numerical simulation and the nonlinear equation asymptotic solution showed good agreement. Presented model could be used for stability analysis and for prediction of amplitudes of chatter in unstable zone. The results of modeling can be exploited in industry for optimal cutting conditions determination.

References

  1. Biermann D., Bleicher F., Heisel, Klocke F., Möhring H. C., Shih A. Deep hole drilling. CIRP Annals, Vol. 67, Issue 2, 2018, p. 673-694. [Publisher]
  2. Sakuma K., Taguchi K., Katsuki A., Takeyama H. Self-guiding action of deep-hole-drilling tools. CIRP Annals, Vol. 30, Issue 1, 1981, p. 311-315. [Publisher]
  3. Ivanov I. I., Voronov S. A. Processing parameters influence on dynamics of vibratory drilling with adaptive control. MATEC Web of Conferences, Vol. 226, 2018, p. 02001. [Publisher]
  4. Kiselev I. A., Zhukov N. A., Selivanov A. N., Barysheva D. V., Voronov S. A., Gouskov A. M. Three-dimensional modeling of deep hole vibratory drilling dynamics. Procedia Engineering, Vol. 176, 2017, p. 50-55. [Publisher]
  5. Ivanov I., Pleshcheev I., Larkin A. Investigation of controlled vibratory drilling dynamics with variable velocity feedback gain. International Conference on Industrial Engineering, 2018, p. 585-594. [Publisher]
  6. Kuts V., Kiselev I., Voronov S. Vibrations reduction in milling by using singular spectral analysis. MATEC Web of Conferences, Vol. 224, 2018, p. 01102. [Publisher]
  7. Voronov S. A., Gouskov A. M., Svetlitsky V. A. Vozbuzhdenie poperechnyh avtokolebanij steblja instrumenta pri glubokom sverlenii (Excitation of lateral self-oscillations of tool shank during deep hole drilling). Raschety na prochnost' (Strength Analysis) Mashinostroenie, 1979, p. 172-182, (in Russian). [Search CrossRef]
  8. Batzer S. A., A. M., Gouskov, Voronov S. A. Modeling vibratory drilling dynamics. Transactions ASME. Journal of Vibration and Acoustics, Vol. 123, Issue 4, 2001, p. 435-443. [Publisher]
  9. Altintas Y. Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design. Cambridge University Press, 2000. [Publisher]
  10. Gouskov A. M., Voronov S. A., Butcher E., Sinha S. C. Nonlinear flexural-torsional vibrations of a gundrilling tool. International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, 2005, p. 971-980. [Publisher]
  11. Bayly P. V., Metzler S. A., Shaut A. J., Young S. G. Theory of torsional chatter in twist drills: model, stability analysis and composition test. Transactions ASME. Journal of Manufacturing Science and Engineering, Vol. 123, 2002, p. 552-561. [Publisher]
  12. Bleicher F., Reiter M., Brier J. Increase of chip removal rate in single-lip deep hole drilling at small diameters by low-frequency vibration support. CIRP Annals, Vol. 68, Issue 1, 2019, p. 93-96. [Publisher]
  13. Nayfeh A. H. Perturbation Methods. John Wiley & Sons, 2008. [Search CrossRef]