Modeling of chatter vibrations in gun drilling process

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. 1Scheme 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 ($F_{pl}$) 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. 2Loading conditions of gun drill cutting edge and drawing of mathematical model of gun drilling system

a) Loading conditions of gun drill cutting edge

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:

where $w$ – lateral displacement of beam axis, $t$ – time, $s$ – axial coordinate, $m_0$ – mass of unit length, $EI$ – flexural stiffness of the beam, $F_a$ – axial force acting at the beam end, $M_b$ – bending moment, due to cutting force eccentricity, $\delta$ – 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$:

where $K_c$ – cutting coefficient, $e_y$ – 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:

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:

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]:

where:

$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:

Term $2\pi K_c{e}_{y}a$ represents constant component of bending moment, which defines the static deflection $x_{st}$. 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:

where:

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-K_c$). Linearized Eq. (8) reads as followed:

10

$\ddot{x}+2d\dot{x}+{p_0}^{2}x+k_1\left(x-x_d\right)=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

$-{\omega }^2+{p_0}^2+k_1\left(1-\cos \left(\omega T\right)\right)=0,2d\omega +k_1\sin \left(\omega T\right)=0.$

Solving these equations in respect to $\omega$, $T$ for different values of $K_c$ we obtain classical periodically altering regions of stability and instability (Fig. 3). Values of parameters used in this paper are $d=$ 0.0493, $e_y=6.17e^{-3}$, $a=1.67e^{-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. 3Stability 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

$\frac{d}{dt}=\frac{\partial }{\partial T_0}+\epsilon \frac{\partial }{\partial T_1}+{\epsilon }^2\frac{\partial }{\partial T_2},$

13

$x=\epsilon x_1\left(T_0,T_1,T_2\right)+{\epsilon }^2{x}_2\left(T_0,T_1,T_2\right)+{\epsilon }^3{x}_3\left(T_0,T_1,T_2\right),$

14

$T=T_{cr}+\epsilon {\tau }_1+{\epsilon }^2{\tau }_2,$

where $T_{cr}$ – critical delay corresponding to stability boundary.

The delayed variable is expanded as:

15

$x_i\left(t-T\right)=x_i\left(T_0-T,T_1-\epsilon T,T_2-{\epsilon }^{2}T\right)=x_{id}+\epsilon \left(-T_{cr}{D}_1{x}_{id}-{\tau }_1{D}_0{x}_{id}\right)+\dots$
$+{\epsilon }^2\left(-T_{cr}{D}_2{x}_{id}-{\tau }_1{D}_1{x}_{id}-{\tau }_2{D}_0{x}_{id}+\frac{1}{2}({T_{cr}}^2{D}_{\mathrm{1,1}}{x}_{id}+2T_{cr}{\tau }_1{D}_{\mathrm{0,1}}{x}_{id}+{\tau }_1{D}_{\mathrm{0,0}}{x}_{id})\right),$

where $D_i=\frac{\partial }{\partial T_i}$, $D_{i,j}=\frac{{\partial }^2}{\partial T_i\partial T_j}$, $x_{id}=x_i\left(T_0-T_{cr},T_1,T_2\right)$.

Substituting Eqs. (12)-(15) and equating the same terms of power $\mathrm{\epsilon }$, the following partial differential equations are obtained:

16

${\epsilon }^1:D_{\mathrm{0,0}}{x}_1+2d{D}_0{x}_1+k{x}_1+k_1\left(x_1-x_{1d}\right)=0,$

17

${\epsilon }^2:D_{\mathrm{0,0}}{x}_2+2d{D}_0{x}_2+k{x}_2+k_1\left(x_2-x_{2d}\right)=\dots$
$=-2D_{\mathrm{0,1}}{x}_1-2d{D}_1{x}_1-k_1\left(T_{cr}{D}_1{x}_{1d}+{\tau }_1{D}_0{x}_{1d}\right)-a_1\left({x_1}^2-{x_{1d}}^2\right)-a_2{x}_1\left(x_1-x_{1d}\right),$

18

${\epsilon }^3:D_{\mathrm{0,0}}{x}_3+2d{D}_0{x}_3+k{x}_3+k_1\left(x_3-x_{3d}\right)=\dots$
$=-2D_{\mathrm{0,1}}{x}_2-2D_{\mathrm{0,2}}{x}_1-D_{\mathrm{1,1}}{x}_1-2d{D}_1{x}_2-2d{D}_2{x}_1-k_1{expr}_1-{expr}_2,$

19

${expr}_1=T_{cr}{D}_2{x}_{1d}+{\tau }_1{D}_1{x}_{1d}+{\tau }_2{D}_0{x}_{1d}-\frac{1}{2}\left({T_{cr}}^2{D}_{\mathrm{1,1}}{x}_{1d}+2T_{cr}{\tau }_1{D}_{\mathrm{0,1}}{x}_{1d}+{\tau }_1{D}_{\mathrm{0,0}}{x}_{1d}\right)$
$+\dots +T_{cr}{D}_1{x}_{2d}+{\tau }_1{D}_0{x}_{2d},$

20

${expr}_2=2a_1\left(x_1{x}_2-x_{1d}{x}_{2d}+x_{1d}\left(T_{cr}{D}_1{x}_{1d}+{\tau }_1{D}_0{x}_{1d}\right)\right)+\dots$
$+a_2\left(x_1\left(x_2-x_{2d}+T_{cr}{D}_1{x}_{1d}+{\tau }_1{D}_0{x}_{1d}\right)+x_2\left(x_1-x_{1d}\right)\right)+a_3{x}_1\left({x_1}^2-{x_{1d}}^2\right).$

General solution of first equation is written as:

21

$x_1=A\left(T_1,T_2\right)\cos \left(\omega T_0\right)+B\left(T_1,T_2\right)\sin \left(\omega T_0\right).$

Quadratic terms do not introduce secular terms in second approximation, therefore ${\tau }_1=0$, and $A$, $B$ are independent of time-scale $T_1$. However, quadratic terms still contribute to third approximation equations. Second approximation has a form:

22

$x_2=A_2\left(T_1\right)\mathit{cos}\left(2\omega T_0\right)+B_2\left(T_1\right)\mathit{sin}\left(2\omega T_0\right)+C_2\left(T_1\right),$

where $A_2$, $B_2$, $C_2$ – 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 $A\left(T_2\right)$, $B\left(T_2\right)$.

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

23

$A\left(T_2\right)=C\left(T_2\right)\mathrm{c}\mathrm{o}\mathrm{s}\phi \left(T_2\right),B\left(T_2\right)=C\left(T_2\right)\mathrm{s}\mathrm{i}\mathrm{n}\phi \left(T_2\right),$

where $C$, $\phi$ – amplitude and phase.

Eliminating $\phi$ from equations we get differential equation for amplitude $C$. To find limit cycle amplitude we assume, that $\frac{\partial }{\partial T_2}C=0$. As a result we will derive the algebraic equation for limit cycle amplitude $C$ calculation. Resulting equation has a form:

24

$f\left(\omega ,T_{cr},d,p_0,a_1,a_2,a_3\right)A^2=g\left(\omega ,T_{cr},d,p_0,a_1,a_2,a_3\right){\tau }_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. 4Time history for x and Fourier spectrum of steady motion. Parameters: Kc= 8000, T= 7.2433

a) Time history for $x$

b) Fourier spectrum of steady motion

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

a) Time history for $x$

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. 6Bifurcation diagram determined by asymptotic analysis and numerical integration

a) Zone A, $K_c=$ 8000, $T_{cr}=$ 7.2258

b) Zone B, $K_c=$ 9000, $T_{cr}=$ 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.