Determination of the stability lobes with multi-delays considering cutter’s helix angle effect for machining process

Abstract

As a common phenomenon, chatter is one of the most important factors that inhibit the improvement of productivity and deteriorate the machined surface quality in milling process. In this article, the mathematical model of the dynamic machining process is first constructed with multi-delays, in which the effect of the cutter’s helix angle on the chatter is considered. And a new integral interpolation method is proposed to predict the stability lobes. Based on this method, the mathematical model which is divided into two parts is calculated respectively with an interpolation method and an integration method. Using the Floquet theory, the stability limits of the dynamic milling system for an arbitrary point can thus be precisely predicted. Subsequently, the convergence rate of this integral interpolation method is analyzed. In order to verify the mathematical model, the stability lobes with uniform and variable cutter pitch angle are computed and the results show a good agreement with published experimental data; also one experiment is conducted to validate the proposed model and the method. Simultaneously, the simulation of a stability lobe with non-zero helix angle is also performed and the results validate the fact that the proposed method is capable of predicting stability limits with high accuracy. Finally, the influences of the cutter and process parameters such as helix angle, pitch angle, down-milling or up-milling and radial cutting depth on the milling chatter are analyzed.

Introduction

In milling process, the resulting surface’s quality is particularly important for some complex key parts and especially for those parts that need to maintain a long work life in poor conditions. However, the chatter often leads to some negative influences on the machining process, such as poor surface quality, low material removal rate ¹ and damage of the structure ² . Therefore, for real nonlinear system, the structure vibration/chatter should be modeled. ³ And in machining process, two extreme cases always happen in the schedule of milling parameters, the first is that the choice of process parameters is over-conservative and the other is to determine those process parameters only depending on the experience of workers. For these reasons, the issues of cutting stability, including the generation mechanism of the chatter and high-accuracy prediction of the stability lobes, are attracting increasingly attentions in order to obtain a good combination of process parameters.^{4 –9}

During the milling process, the dominant reason for the chatter is the dynamic interactions between the cutting tool and part.^{10 –14} For example, the material left by the current cutting edge may be cut by the 1st, 2nd, 3rd,…,nth cutter edge previous to the current one. If the material is cut by the 1st previous cutter edge, the single regenerative effect of chatter vibration occurs; otherwise, the chatter with multi-delays is induced. ¹⁵ Actually, the chatter with multi-regenerative effect is more common due to the reasons of always existing cutter runout and the usage of the cutter with non-uniform pitch angle. Many researchers have made substantial efforts on the aspect of milling stability prediction. The control equation of the chatter with single regenerative effect is first presented by Tlusty and Polacek ¹⁶ and Tobias. ¹⁷ H Moradi et al. ¹⁸ analyze the cutting forced vibration of the milling process using this theory. B Balachandran and colleagues^{19 –21} take the dynamic cutting process as a nonlinear system and the stability limits are analyzed. Y Altintas and colleagues^{22 –24} and E Budak and colleagues^{25 –28} formulate the dynamic system as a series of delay-differential equations (DDEs), which has been become a classical model used to compute stability lobes by researchers. They transform this model from time domain to frequency domain utilizing Fourier transformation method. Then, using the zero-order Fourier expansion term, the expressions in frequency domain can be simplified. Finally, a explicate equation is obtained to calculate the relationship between axial limit milling depth and the spindle’s rotational speed under the condition of a radial cutting depth, and then stability lobes are constructed from this equation.

It is found that the process of the transformation from time domain to frequency domain has the problem of precision loss for predicting stability lobes; so the calculation method using Floquet’s theory in time domain is presented. ²⁹ Utilizing this theory, the time domain discrete methods with finite element ideology are introduced. Insperger and colleagues^{30 –33} have made significant efforts on this aspect, and the semi-discretization methods (SDM) with different order are presented by them to predict the stability lobes. The method in the literatures ³¹ has been proved by Catania and Mancinelli ³⁴ and has already aroused more and more attentions. In order to improve the prediction accuracy and computation efficiency, some promising methods are also proposed. Ding et al.^{35 –37} present a new full-discretization method (FDM). Zhang et al. ³⁸ introduce a variable-step interpolation method for milling chatter stability prediction with multiple delays. He also studies stability analysis in milling of thin-walled workpieces with emphasis on the structural effect. ³⁹ Luo et al. ⁴⁰ research material removal process optimization for milling of flexible workpiece considering machining stability. Zaghbani et al. ⁴¹ use the stability technology on a robotic machining system. In our previous works,^42,43 a high order full-discretization method is proposed to predict stability lobes with high accuracy. Huang et al. ⁴⁴ present a liner approximation of acceleration method for stability lobes prediction.

The chatter induced by multi-regenerative effect is also studied and analyzed. The literatures^38,42,43 have introduced the cutter runout’s effect on the milling systems. And researchers find that the cutter runout not only affects the actual pitch angles but also affects the actual cutting depths. That means when cutter runout exists, the actual pitch angle is not equal to the ideal cutter pitch angle and the actual cutting depth is not equal to the ideal cutting depth. As a result, this will induce the multi-regenerative effect when chatter happens in the milling process. Insperger et al. ⁴⁵ find that the cutter runout has little effect on the stability limits, but has great effect on the chatter frequencies. This is proved using our method in previous works. ⁴³ Using some special cutters with non-uniform distribution of cutter edges, the chatter with multi-regenerative effect can be depressed more easily in machining process compared with the process machined using the cutter with uniform distribution of cutter edges. Altintas and Budak ²² find this phenomenon and then give an classical example which has been proved utilizing experiments. In this example, the cutter with the variable pitch angle of 70°-110°-70°-110° is taken to mill a block, and a fact has been found that sometimes using the cutter with variable pitch angle is helpful for the stability of the milling process. And this example and view are adopted in many literatures^{26,35 –38,42,43} including the method using the time domain discretization method. In Altintas and Budak’s ²² literature, they find that the helix angle of the cutter has no effect on the stability lobes of the milling process. This principle is correct in special conditions such as small radial cutting depth. ³⁴ The researchers in Altintas and colleagues,^46,47 Dombovari and Stepan ⁴⁸ and Eftekhari et al. ² also find this problem and they consider the helix angles’ effect on the stability lobes. Wan et al. ⁴⁹ modify the mathematical model which can consider the effect of the cutter helix angle on the stability lobes, and some important conclusions are given.

In this work, a mathematical model considering multi-regenerative effect combining the cutter’s helix angle is deduced as the governing equation of the chatter in machining process. A new method to compute the stability lobes using this mathematical model is presented. Then, this governing model is divided into two parts for the prediction of the stability limits, the first part is simplified using interpolation method and the other one is simplified utilizing integration way. For next step, the convergence rate of this method is analyzed. And four examples are conduced to verify the validity of the method. The first one considers uniform pitch angle of the cutter. The second example takes variable pitch angle into consideration. The helix angle of the cutter is combined in the third example to verify the validity of the mathematical model. These examples’ data have been proved using the experiments in the literatures. One experiment is conducted in the last example to verify the model and the method. The results show that the method proposed in this article can predict the stability lobes with high accuracy. Finally, some parameters’ effects on the stability lobes are analyzed.

Governing equation

The classical mathematical model of the chatter in the milling process ignores the cutter helix angle’s effect. That means the helix angle of a cutter is taken as zero degree. And this situation cannot happen in real machining process. Therefore, the traditional model for the chatter should be made some proper modification. Figure 1 shows the chatter induced by multi-regenerative effect for a simple cutter with non-zero helix angle.

Figure 1.

(a) Schematic diagram of dynamic milling system and (b) the simple end cutter with developed diagram of the cutter edge about the cutter axis.

Figure 1(a) stands for the dynamic behavior of the milling process for a cutter with zero degree of cutter helix angle. However, for the simple cutter in Figure 1(b), the schematic diagram in Figure 1(a) cannot stand for the dynamic milling system again. Therefore, the cutter should be divided into a series of cutting elements along the cutter axis, and then every cutting element can be taken as a small cutter with zero cutter helix angle. In other words, the dynamic behavior of every cutting element can be represented using Figure 1(a). Therefore, for lth cutting element, the mathematical model of dynamic milling system can be expressed using the equation in our previous work ⁴³

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{M}}_l\frac{d^2\mathit{q}(t)}{\mathrm{dt}}+{\mathbf{C}}_l\frac{d\mathit{q}(t)}{\mathrm{dt}}+{\mathbf{K}}_l\mathit{q}(t)}\\ \multicolumn{1}{c}{=a_l\sum _{j=1}^N\left({\mathbf{K}}_{c,l,j}\left(t\right)[\mathit{q}(t-{\tau }_j)-\mathit{q}\left(t\right)]\right)}\end{array}$$

(1)

In this equation, M _l , C _l and K _l stand for the modal mass matrix, modal damping matrix and modal stiffness matrix for lth cutting element, respectively; a _l is the thickness of the lth cutting element; N is the number of the edges contained by the cutter; ${\mathbf{K}}_{c,l,j}$ is the cutting force parameters for jth cutter edge at lth cutting element; ${\tau }_j$ is the period of jth cutter edge and $q(t)$ and $\mathit{q}(t-{\tau }_j)$ are the dynamic displacement vectors of the current cutting element and the jth previous cutting element. All the cutting elements’ dynamic behavior can be expressed using equation (1). Then, after all the equations for all cutting elements’ chatter mathematical model are added, the governing equation for the cutter can be obtained and expressed utilizing

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{M}\frac{d^2\mathit{q}(t)}{\mathrm{dt}}+\mathbf{C}\frac{d\mathit{q}(t)}{\mathrm{dt}}+\mathbf{K}\mathit{q}(t)}\\ \multicolumn{1}{c}{=\sum _{l=1}^M\left(a_l\sum _{j=1}^N\left({\mathbf{K}}_{c,l,j}\left(t\right)[\mathit{q}(t-{\tau }_j)-\mathit{q}\left(t\right)]\right)\right)}\end{array}$$

(2)

where M is the amount of the cutting elements included by the cutter; M, C and K stand for the modal mass matrix, modal damping matrix and modal stiffness matrix for the cutter with non-zero helix angle, respectively, and these parameters can be expressed by the following equations

$$\{\begin{array}{ccc}\hfill \mathbf{M}=\sum _{l=1}^M{\mathbf{M}}_l,\hfill & \hfill \mathbf{C}=\sum _{l=1}^M{\mathbf{C}}_l,\hfill & \hfill \mathbf{K}=\sum _{l=1}^M{\mathbf{K}}_l\hfill \\ \multicolumn{1}{c}{\mathbf{M}=\left[\begin{array}{cc}\hfill m_x\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill m_y\hfill \end{array}\right],\mathbf{C}=\left[\begin{array}{cc}\hfill c_x\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill c_y\hfill \end{array}\right],\mathbf{K}=\left[\begin{array}{cc}\hfill k_x\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill k_y\hfill \end{array}\right]}\\ \multicolumn{1}{c}{\mathit{q}={\left[x\left(t\right)y\left(t\right)\right]}^T,{\mathbf{K}}_{c,l,j}\left(t\right)=\frac{1}{2}\left[\begin{array}{cc}\hfill a_{l,j,\mathrm{xx}}\hfill & \hfill a_{l,j,\mathrm{xy}}\hfill \\ \multicolumn{1}{c}{a_{l,j,\mathrm{yx}}}& \hfill a_{l,j,\mathrm{yy}}\hfill \end{array}\right]}\end{array}$$

(3)

$$\{\begin{array}{c}\hfill a_{l,j,\mathrm{xx}}=g\left({\phi }_{l,j}\right)[K_r\sin \left(2{\phi }_{l,j}\right)-K_t(\cos \left(2{\phi }_{l,j}\right)-1)]\hfill \\ \multicolumn{1}{c}{a_{l,j,\mathrm{xy}}=g\left({\phi }_{l,j}\right)[K_t\sin \left(2{\phi }_{l,j}\right)+K_r(\cos \left(2{\phi }_{l,j}\right)+1)]}\\ \multicolumn{1}{c}{a_{l,j,\mathrm{yx}}=g\left({\phi }_{l,j}\right)[K_t\sin \left(2{\phi }_{l,j}\right)-K_r(\cos \left(2{\phi }_{l,j}\right)-1)]}\\ \multicolumn{1}{c}{a_{l,j,\mathrm{yy}}=g\left({\phi }_{l,j}\right)[K_r\sin \left(2{\phi }_{l,j}\right)+(\cos \left(2{\phi }_{l,j}\right)-1)K_t]}\end{array}$$

(4)

K_t and K_r are the static average cutting force coefficients along the tangent and radial direction of the cutter, respectively; (m_x, m_y), (c_x, c_y) and (k_x, k_y) are the elements along x and y direction for the matrixes M, C and K; $a_{l,j,\mathrm{xx}}$ , $a_{l,j,\mathrm{xy}}$ , $a_{l,j,\mathrm{yx}}$ and $a_{l,j,\mathrm{yy}}$ are the four elements of the ${\mathbf{K}}_{c,l,j}(t)$ and $x(t)$ and $y(t)$ are the elements of the vector $\mathit{q}(t)$ . $g({\phi }_{l,j})$ is a unit step function (0-1 function) which determines whether the cutter edge is removing material or not. And this function can be expressed as

$$\{\begin{array}{cc}\hfill g\left({\phi }_{l,j}\right)=1\hfill & \hfill \mathrm{when}{\phi }_{\mathrm{st}}\le {\phi }_{l,j}\le {\phi }_{\mathrm{ex}}\hfill \\ \multicolumn{1}{c}{g\left({\phi }_{l,j}\right)=0}& \hfill \mathrm{when}{\phi }_{\mathrm{st}}>{\phi }_{l,j}\mathrm{or}{\phi }_{l,j}>{\phi }_{\mathrm{ex}}\hfill \end{array}$$

(5)

where ${\phi }_{\mathrm{st}}$ and ${\phi }_{\mathrm{ex}}$ are the entering and exiting cutting angle, respectively; ${\phi }_{l,j}$ is the rotation angle for jth cutter edge at lth cutting element, and it should transform into [0, 2π) using

$$\begin{array}{l}while{\phi }^{i,j}>2\pi ,{\phi }^{i,j}={\phi }^{i,j}-2\pi \\ andwhile{\phi }^{i,j}<0,{\phi }^{i,j}={\phi }^{i,j}+2\pi \end{array}$$

When equation (2) is obtained, the simple cutter’s dynamic behavior with non-zero helix angle can be studied. Therefore, the computation of equation (2) becomes a key for studying the stability lobes of a common cutter with non-zero helix angle.

Integral interpolation method

In order to predict the stability lobes, we introduce a new way using integral interpolation method; this method provides a new ideology combining the advantages of the integration and interpolation theories. First, we should simplify equation (2) using the following equation

$$\begin{array}{c}\multicolumn{1}{c}{d\mathit{x}\left(t\right)/\mathrm{dt}={\mathbf{A}}_0\mathit{x}\left(t\right)}\\ \multicolumn{1}{c}{+\sum _{l=1}^M\left(\sum _{j=1}^N[{\mathbf{A}}_{l,j}\left(t\right)\mathit{x}\left(t\right)+{\mathbf{B}}_{l,j}\left(t\right)\mathit{x}(t-{\tau }_j)]\right)}\end{array}$$

(6)

where the constant matrix A₀ is invariant; A _l,j (t) and B _l,j (t) are two periodic functions with the period τ_j for jth cutter edge and x (t) is the state variable of the system. A₀, x (t), A _l,j (t) and B _l,j (t) can be can be calculated using

$${\mathbf{A}}_0=\left[\begin{array}{cc}\hfill -\frac{{\mathbf{M}}^{-1}\mathbf{C}}{2}\hfill & \hfill {\mathbf{M}}^{-1}\hfill \\ \multicolumn{1}{c}{\frac{\mathbf{C}{\mathbf{M}}^{-1}\mathbf{C}}{4}-\mathbf{K}}& \hfill -\frac{{\mathbf{M}}^{-1}\mathbf{C}}{2}\hfill \end{array}\right]$$

(7)

$$\mathit{x}\left(t\right)={\left[\begin{array}{cc}\hfill \mathit{q}\left(t\right)\hfill & \hfill \mathbf{M}\frac{d\mathit{q}\left(t\right)}{\mathrm{dt}}+\frac{\mathbf{C}\mathit{q}\left(t\right)}{2}\hfill \end{array}\right]}^T$$

(8)

$${\mathbf{B}}_{l,j}\left(t\right)=-{\mathbf{A}}_{l,j}\left(t\right)=a_l\left[\begin{array}{cc}\hfill \mathbf{0}\hfill & \hfill \mathbf{0}\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_{c,l,j}\left(t\right)}& \hfill \mathbf{0}\hfill \end{array}\right]$$

(9)

The spindle rotational speed is Ω. Then, the rotational period T of the spindle can be obtained. Each cutter edge’s period can also be computed using

$$\{\begin{array}{c}\hfill \mathrm{T}=60/\Omega =\sum _{j=1}^N{\tau }_j\hfill \\ \multicolumn{1}{c}{{\tau }_j={\varphi }_{j-1,j}\mathrm{T}/2\pi }\end{array}$$

(10)

${\varphi }_{j-1,j}$ is the pitch angle between jth and (j − 1)th cutter edge; if j = 1, ${\varphi }_{j-1,j}$ is the pitch angle between 1st and Nth cutter edge. Then, the important step of computing stability lobes is to divide the spindle rotational period T into m time intervals, and jth cutter edge’s period τ_j into m_j time intervals, they meet the following conditions

$$m_j=\mathrm{round}({\tau }_j/h),h=T/m,m=\sum _{j=1}^N{m}_j$$

(11)

where h is time span for the time intervals round() stands for the rounding operation, such as round(4.2) = 4 and round(3.6) = 4.

For every time interval in the spindle period, the solution of equation (6) can be expressed using special formula. Here, taking the kth time interval in the jth cutter edge’s period as an example, this formula can be evaluated using

$$\{\begin{array}{c}\hfill {\mathit{x}}_{m(j,k+1)}=e^{{\mathbf{A}}_0(t_{m(j,k+1)}-t_{m(j,k)})}{\mathit{x}}_{m(j,k)}\hfill \\ \multicolumn{1}{c}{+\sum _{l=1,j=1}^{M,N}\left(\underset{t_{m(j,k)}}{\overset{t_{m(j,k+1)}}{\int }}{e}^{{\mathbf{A}}_0(t_{m(j,k+1)}-\epsilon )}\left[\begin{array}{c}\hfill {\mathbf{A}}_{l,j}\left(\epsilon \right)\mathit{x}\left(\epsilon \right)\hfill \\ \multicolumn{1}{c}{+{\mathbf{B}}_{l,j}\left(\epsilon \right)\mathit{x}(\epsilon -{\tau }_j)}\end{array}\right]d\epsilon \right)}\\ \multicolumn{1}{c}{m(j,k+1)=\sum _{p=1}^{j-1}{m}_p+k+1}\end{array}$$

(12)

With ϵ ∈ [t_m(j,k), t_{m(j,k + 1)}], ${\mathit{x}}_{m(j,k+1)}$ stands for the value of $\mathit{x}(t)$ at the time $t_{m(j,k+1)}$ . Then, a new ideology is introduced. The basis of this ideology is that the integration part in equation (12) is divided into two parts; the first part and the second part can be expressed using the following equations

$$\{\begin{array}{c}\hfill I\mathrm{part}:\underset{t_{m(j,k)}}{\overset{t_{m(j,k+1)}}{\int }}{e}^{{\mathbf{A}}_0(t_{m(j,k+1)}-\epsilon )}{\mathbf{A}}_{l,j}\left(\epsilon \right)\mathit{x}\left(\epsilon \right)d\epsilon \hfill \\ \multicolumn{1}{c}{\mathrm{II}\mathrm{part}:\underset{t_{m(j,k)}}{\overset{t_{m(j,k+1)}}{\int }}{e}^{{\mathbf{A}}_0(t_{m(j,k+1)}-\epsilon )}{\mathbf{B}}_{l,j}\left(\epsilon \right)\mathit{x}(\epsilon -{\tau }_j)d\epsilon }\end{array}$$

(13)

The first part can be simplified utilizing the interpolation method and the second part can be simplified using integration theory. It must be point out that part II can be taken as part I, and then part I is changed into part II.

And here, for the first part, we use the third-order Lagrange or Newton interpolation method, if using the Lagrange method. Here, ${\mathbf{A}}_{l,j}(\epsilon )$ and $\mathit{x}(\epsilon )$ can be simplified using the following equation

$$\{\begin{array}{c}{A}_{l,j}\left(\epsilon \right)\approx \left(\raisebox{1ex}{$\left(h-\zeta \right)$}\!\left/ \!\raisebox{-1ex}{$h$}\right.\right)A_{l,m\left(j,k\right)}+\left(\raisebox{1ex}{$\zeta $}\!\left/ \!\raisebox{-1ex}{$h$}\right.\right)A_{l,m\left(j,k+1\right)}\\ x\left(\epsilon \right)\approx a{x}_{m\left(j,k-2\right)}+b{x}_{m\left(j,k-1\right)}+c{x}_{m\left(j,k\right)}+d{x}_{m\left(j,k+1\right)}\end{array}$$

(14)

With ζ ∈ [0, h], ${\mathbf{A}}_{l,m(j,k)}$ and ${\mathbf{A}}_{l,m(j,k+1)}$ are the matrixes’ values at the t_m(j,k) and t_{m(j,k + 1)} for ${\mathbf{A}}_{l,j}(t)$ and a, b, c, d are the Lagrange coefficients which can be represented using the following equation

$$\{\begin{array}{c}\hfill a=-\frac{{\zeta }^3-\zeta h^2}{6h^3}\hfill \\ \multicolumn{1}{c}{b=\frac{{\zeta }^3-2h^2\zeta +h{\zeta }^2}{2h^3}}\\ \multicolumn{1}{c}{c=-\frac{{\zeta }^3-h^2\zeta +2h{\zeta }^2-2h^3}{2h^3}}\\ \multicolumn{1}{c}{d=\frac{{\zeta }^3+3h{\zeta }^2+2h^2\zeta }{6h^3}}\end{array}$$

(15)

And ${\mathbf{A}}_{l,j}(\epsilon )$ and $\mathbf{x}(\epsilon )$ can also be simplified using the Newton interpolation method which expressed in our previous works.^42,43 Then, using the direct integration method, the first part can be simplified for the subsequent computations.

Thus, the first part in equation (13) can be calculated using the following analytical expressions

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{part}I={\mathbf{F}}_{l,1}{\mathit{x}}_{m(j,k-2)}+{\mathbf{F}}_{l,2}{\mathit{x}}_{m(j,k-1)}}\\ \multicolumn{1}{c}{+{\mathbf{F}}_{l,3}{x}_{m(j,k)}+{\mathbf{F}}_{l,4}{x}_{m(j,k+1)}}\end{array}$$

(16)

${\mathbf{F}}_{l,1}$ , ${\mathbf{F}}_{l,2}$ , ${\mathbf{F}}_{l,3}$ and ${\mathbf{F}}_{l,4}$ are the coefficient matrix for the $x_{m(j,k-2)}$ , ${\mathit{x}}_{m(j,k-1)}$ , ${\mathit{x}}_{m(j,k)}$ and ${\mathit{x}}_{m(j,k+1)}$ . The parameters in equation (16) can be calculated using the following equation

$$\{\begin{array}{c}\hfill {\mathbf{F}}_{l,1}=\frac{{\mathbf{f}}_0}{6h^4}\left[\begin{array}{c}\hfill (h^3{\mathbf{f}}_1-h^2{\mathbf{f}}_2-h{\mathbf{f}}_3+{\mathbf{f}}_4){\mathbf{A}}_{l,m(j,k)}\hfill \\ \multicolumn{1}{c}{+(h^2{\mathbf{f}}_2-{\mathbf{f}}_4){\mathbf{A}}_{l,m(j,k+1)}}\end{array}\right]\hfill \\ \multicolumn{1}{c}{{\mathbf{F}}_{l,2}=\frac{{\mathbf{f}}_0}{2h^4}\left[\begin{array}{c}\hfill (3h^2{\mathbf{f}}_2-2h^3{\mathbf{f}}_1-{\mathbf{f}}_4){\mathbf{A}}_{l,m(j,k)}\hfill \\ \multicolumn{1}{c}{+({\mathbf{f}}_4+h{\mathbf{f}}_3-2h^2{\mathbf{f}}_2){\mathbf{A}}_{l,m(j,k+1)}}\end{array}\right]}\\ \multicolumn{1}{c}{{\mathbf{F}}_{l,3}=\frac{{\mathbf{f}}_0}{2h^4}\left[\begin{array}{c}\hfill ({\mathbf{f}}_4+h{\mathbf{f}}_3-3h^2{\mathbf{f}}_2-h^3{\mathbf{f}}_1+2h^4\mathbf{f}){\mathbf{A}}_{l,m(j,k)}\hfill \\ \multicolumn{1}{c}{-({\mathbf{f}}_4+2h{\mathbf{f}}_3-h^2{\mathbf{f}}_2-2h^3{\mathbf{f}}_1){\mathbf{A}}_{l,m(j,k+1)}}\end{array}\right]}\\ \multicolumn{1}{c}{{\mathbf{F}}_{l,4}=\frac{{\mathbf{f}}_0}{6h^4}\left[\begin{array}{c}\hfill (2h^3{\mathbf{f}}_1+h^2{\mathbf{f}}_2-2h{\mathbf{f}}_3-{\mathbf{f}}_4){\mathbf{A}}_{l,m(j,k)}\hfill \\ \multicolumn{1}{c}{+({\mathbf{f}}_4+3h{\mathbf{f}}_3+2h^2{\mathbf{f}}_2){\mathbf{A}}_{l,m(j,k+1)}}\end{array}\right]}\end{array}$$

(17)

where f, f₀, f₁, f₂, f₃ and f₄ are intermediate matrixes and can be expressed using the following equations

$$\{\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbf{f}=\underset{0}{\overset{h}{\int }}{e}^{-{\mathbf{A}}_0\zeta }d\zeta ,{\mathbf{f}}_0=e^{{\mathbf{A}}_{0}h}\hfill \\ \multicolumn{1}{c}{{\mathbf{f}}_1=\underset{0}{\overset{h}{\int }}{e}^{-{\mathbf{A}}_0\zeta }\zeta d\zeta ,{\mathbf{f}}_2={\int }_0^h{e}^{-{\mathbf{A}}_0\zeta }{\zeta }^{2}d\zeta }\end{array}\hfill \\ \multicolumn{1}{c}{{\mathbf{f}}_3=\underset{0}{\overset{h}{\int }}{e}^{-{\mathbf{A}}_0\zeta }{\zeta }^{3}d\zeta ,{\mathbf{f}}_4={\int }_0^h{e}^{-{\mathbf{A}}_0\zeta }{\zeta }^{4}d\zeta }\end{array}$$

(18)

And f, f₁, f₂, f₃ and f₄ can be calculated analytically utilizing f₀

$$\{\begin{array}{c}\hfill \mathbf{f}={{\mathbf{A}}_0}^{-1}(\mathbf{I}-{{\mathbf{f}}_0}^{-1}),{\mathbf{f}}_1={{\mathbf{A}}_0}^{-1}(\mathbf{f}-h{{\mathbf{f}}_0}^{-1})\hfill \\ \multicolumn{1}{c}{{\mathbf{f}}_2={{\mathbf{A}}_0}^{-1}(2{\mathbf{f}}_1-h^2{{\mathbf{f}}_0}^{-1}),{\mathbf{f}}_3={{\mathbf{A}}_0}^{-1}(3{\mathbf{f}}_3-h^3{{\mathbf{f}}_0}^{-1})}\\ \multicolumn{1}{c}{{\mathbf{f}}_4={{\mathbf{A}}_0}^{-1}(4{\mathbf{f}}_3-h^4{{\mathbf{f}}_0}^{-1})}\end{array}$$

(19)

For the second part, the numerical integration method is used in this article, two-point integration method is adopted to simplify this part. And the simplification equation can evaluated utilizing the formula with

$$\begin{array}{c}\multicolumn{1}{c}{\underset{t_{m(j,k)}}{\overset{t_{m(j,k+1)}}{\int }}{e}^{{\mathbf{A}}_0(t_{m(j,k+1)}-\epsilon )}{\mathbf{B}}_{l,j}\left(\epsilon \right)\mathit{x}(\epsilon -{\tau }_j)d\epsilon =}\\ \multicolumn{1}{c}{\frac{h}{2}(e^{{\mathbf{A}}_{0}h}{\mathbf{B}}_{l,m(j,k)}\mathit{x}(t_{m(j,k)}-{\tau }_j)}\\ \multicolumn{1}{c}{+{\mathbf{B}}_{l,m(j,k+1)}\mathit{x}(t_{m(j,k+1)}-{\tau }_j))}\end{array}$$

(20)

Then, for the reason of the periodic function of the $\mathit{x}(t-{\tau }_j)$ , the periodic condition in the following equation must be met

$$\{\begin{array}{c}\hfill \mathit{x}(t_{m(j,k)}-{\tau }_j)={\mathit{x}}_{m(j,k-m_j)}\hfill \\ \multicolumn{1}{c}{\mathit{x}(t_{m(j,k+1)}-{\tau }_j)={\mathit{x}}_{m(j,k-m_j+1)}}\end{array}$$

(21)

Therefore, equation (20) can be expressed in another way, which is

$$\begin{array}{c}\multicolumn{1}{c}{\underset{t_{m(j,k)}}{\overset{t_{m(j,k+1)}}{\int }}{e}^{{\mathbf{A}}_0(t_{m(j,k+1)}-\epsilon )}{\mathbf{B}}_{l,j}\left(\epsilon \right)\mathit{x}(\epsilon -{\tau }_j)d\epsilon \approx }\\ \multicolumn{1}{c}{\frac{h}{2}({\mathbf{f}}_0{\mathbf{B}}_{l,m(j,k)}{\mathit{x}}_{m(j,k-m_j)}+{\mathbf{B}}_{l,m(j,k+1)}{\mathit{x}}_{m(j,k-m_j+1)})}\end{array}$$

(22)

Finally, equation (12) can be evaluated using following equation after simplification of part I and part II with formulas from equations (14) to (22)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{x}}_{m(j,k+1)}={\mathbf{f}}_0{\mathit{x}}_{m(j,k)}+\sum _{l=1,j=1}^{M,N}{\mathbf{F}}_{l,1}{\mathit{x}}_{m(j,k-2)}}\\ \multicolumn{1}{c}{+\sum _{l=1,j=1}^{M,N}{\mathbf{F}}_{l,2}{\mathit{x}}_{m(j,k-1)}}\\ \multicolumn{1}{c}{+\sum _{l=1,j=1}^{M,N}{\mathbf{F}}_{l,3}{\mathit{x}}_{m(j,k)}+\sum _{l=1,j=1}^{M,N}{\mathbf{F}}_{l,4}{\mathit{x}}_{m(j,k+1)}}\\ \multicolumn{1}{c}{+\sum _{l=1,j=1}^{M,N}\frac{h}{2}{\mathbf{f}}_0{\mathbf{B}}_{l,m(j,k)}{\mathit{x}}_{m(j,k-m)}}\\ \multicolumn{1}{c}{+\sum _{l=1,j=1}^{M,N}\frac{h}{2}{\mathbf{B}}_{l,m(j,k+1)}{\mathit{x}}_{m(j,k-m+1)}}\end{array}$$

(23)

Then, equation (23) can be expressed after merging the same items with the formula

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{x}}_{m(j,k+1)}={(\mathbf{I}-{\mathbf{F}}_0)}^{-1}(({\mathbf{F}}_1+{\mathbf{f}}_0){\mathit{x}}_{m(j,k)}+{\mathbf{F}}_2{\mathit{x}}_{m(j,k-1)})}\\ \multicolumn{1}{c}{+{(\mathbf{I}-{\mathbf{F}}_0)}^{-1}({\mathbf{F}}_3{\mathit{x}}_{m(j,k-2)}+\sum _{j=1}^N\left({\mathbf{F}}_{m(j,k-m_j)}{\mathit{x}}_{m(j,k-m_j)}\right))}\\ \multicolumn{1}{c}{+{(\mathbf{I}-{\mathbf{F}}_0)}^{-1}\sum _{j=1}^N\left({\mathbf{F}}_{m(j,k-m_j+1)}{\mathit{x}}_{m(j,k-m_j+1)}\right)}\end{array}$$

(24)

With

$$\{\begin{array}{c}\hfill {\mathbf{F}}_0=\sum _{l=1,j=1}^{M,N}{\mathbf{F}}_{l,4},{\mathbf{F}}_1=\sum _{l=1,j=1}^{M,N}{\mathbf{F}}_{l,3}\hfill \\ \multicolumn{1}{c}{{\mathbf{F}}_2=\sum _{l=1,j=1}^{M,N}{\mathbf{F}}_{l,2},{\mathbf{F}}_3=\sum _{l=1,j=1}^{M,N}{\mathbf{F}}_{l,1}}\\ \multicolumn{1}{c}{{\mathbf{F}}_{m(j,k-m_j)}=\sum _{l=1}^M\frac{h}{2}{\mathbf{f}}_0{\mathbf{B}}_{l,m(j,k)}}\\ \multicolumn{1}{c}{{\mathbf{F}}_{m(j,k-m_j+1)}=\sum _{l=1}^M\frac{h}{2}{\mathbf{B}}_{l,m(j,k+1)}}\end{array}$$

(25)

Then, equation (24) named as iterated function can be expressed in matrix form

$${\mathit{Y}}_{m(j,k+1)}={\mathbf{D}}_{m(j,k)}{\mathit{Y}}_{m(j,k)}$$

(26)

In this equation, the column vector $Y_{m(j,k+1)}$ can be expressed by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{Y}}_{m(j,k+1)}=[{\mathit{x}}_{m(j,k+1)}{\mathit{x}}_{m(j,k)}{\mathit{x}}_{m(j,k-1)}}\\ \multicolumn{1}{c}{\cdots {\mathit{x}}_{m(j,k-m_1+1)}\cdots }\\ \multicolumn{1}{c}{{\mathit{x}}_{m(j,k-m_N+1)}\cdots {{\mathit{x}}_{m(j,k-m+1)}]}^T}\end{array}$$

(27)

The matrix ${\mathbf{D}}_{m(j,k)}$ can be formulated as the following expressions

$$\begin{array}{l}{D}_{m\left(j,k\right)}=\\ \left[\begin{array}{l}{U}_k{U}_{k-1}{U}_{k-2}0\cdots 0U_{k-m_j+1}{U}_{k-m_j}0\cdots 0U_{k-m_N}{U}_{k-m_N+1}\cdots \\ I00\cdots \cdots \cdots \cdots \\ 0I0\cdots \cdots \cdots \cdots \\ ⋮⋮⋮⋮⋮⋮\cdots \\ 000\cdots I0\cdots \end{array}\right]\end{array}$$

(28)

The parameters in above equation can be calculated utilizing following formulas

$$\{\begin{array}{c}\hfill {\mathbf{U}}_k={(\mathbf{I}-{\mathbf{F}}_0)}^{-1}{\mathbf{F}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{U}}_{k-1}={(\mathbf{I}-{\mathbf{F}}_0)}^{-1}{\mathbf{F}}_2}\\ \multicolumn{1}{c}{{\mathbf{U}}_{k-2}={(\mathbf{I}-{\mathbf{F}}_0)}^{-1}{\mathbf{F}}_3}\\ \multicolumn{1}{c}{{\mathbf{U}}_{k-m_j+1}={(\mathbf{I}-{\mathbf{F}}_0)}^{-1}{\mathbf{F}}_{m(j,k-m_j)}}\\ \multicolumn{1}{c}{{\mathbf{U}}_{k-m_j}={(\mathbf{I}-{\mathbf{F}}_0)}^{-1}{\mathbf{F}}_{m(j,k-m_j+1)}}\end{array}$$

(29)

After iterated function (28) is known, the column vector Y _m can be represented using Y ₀ with the equation

$${\mathit{Y}}_m={\mathbf{D}}_{m-1}\cdots {\mathbf{D}}_2{\mathbf{D}}_1{\mathbf{D}}_0{\mathit{Y}}_0$$

(30)

A new matrix $\mathbf{\Phi }$ can be defined as transition function matrix over the period T by coupling each approximate solution through the discrete map D_i i = 1, 2,…, m − 1 to get

$$\mathbf{\Phi }={\mathbf{D}}_{m-1}\cdots {\mathbf{D}}_2{\mathbf{D}}_1{\mathbf{D}}_0$$

(31)

Thus, using the above transition function and Floquet theory, the stability lobes with multiple delays can be calculated.

Convergence property of the integral interpolation method

As an important factor which influences the computation time and efficiency, the convergence of the critical eigenvalues should be studied. Here, the 1-degree-of-freedom (DOF) down-milling process is chosen and the parameters of the dynamic machining system are

$$\{\begin{array}{c}\hfill m_x=0.03993\mathrm{kg}\hfill \\ \multicolumn{1}{c}{\omega =922\mathrm{Hz}}\\ \multicolumn{1}{c}{c_x=0.011}\\ \multicolumn{1}{c}{{\mathrm{K}}_t=6\mathrm{e}8\mathrm{N}/{\mathrm{m}}^2}\\ \multicolumn{1}{c}{{\mathrm{K}}_r=2e8\mathrm{N}/{\mathrm{m}}^2}\end{array}$$

(32)

The cutter has two cutter edges. Figure 2 shows the relationship between the modulus of the critical eigenvalues |u| and the number m of the time intervals. Four groups of the machining parameters in Figure 2 are chosen as follows: A: {0.5, 4 mm, 9130 r/min}, B: {0.5, 5 mm, 4980 r/min}, C: {0.5, 6 mm, 2965 r/min}, D: {0.5, 7 mm, 2380 r/min}. The first parameter for these groups is the ratio between radial cutting depth and the cutter’s diameter, the second one is the axial cutting depth and the last one is the rotational speed of the cutter about the spindle.

Figure 2.

Convergence property of the method presented in this article.

For these examples, the convergence of the critical eigenvalues |u| for presented method in this article can be observed in Figure 2. The critical eigenvalue |u| is taken as the function of the time intervals number m. The plots show that how the critical eigenvalues |u| converge to a fix value. The difference ratios between the critical eigenvalues |u| when m = 40 and m = 60 are 0.004848 for A, 0.02011 for B, 0.01407 for C and 0.074796 for D.

Model and method verification

Stability lobes prediction for uniform pitch cutter with zero helix angle

The validity of the method presented in this article is verified using the experiment data in the literature. ²⁷ The 2-DOF dynamic milling system is chosen and the properties of the system machined by the cutter which contains three cutter edges are illustrated in Table 1, half-immersion (up-milling) case are chosen.

Table 1.

Parameters of the dynamic system.

Directions	Natural frequency (Hz)	Damping ratios	Modal stiffness (kN/m)
X	603	0.039	5600
Y	666	0.035	5700

The cutting force coefficients are K_t = 600 MPa, K_r = 42 MPa. The experiment result expressed in the literature ²⁷ and the predicting stability lobe utilizing the proposed method are illustrated in Figure 3. This stability limit lobe predicted by the presented method agrees well with that presented in the literature.

Figure 3.

Comparison of the method in this article and that in the literature. ²⁸

And, four points containing A1 (6000 r/min, 2 mm), B1 (12000 r/min, 2 mm), C1(8000 r/min, 2 mm) and D1(18000 r/min, 2 mm) are simulated to compute the vibration displacement in x-axis direction using the dde23 function in the software of MATLAB. It can be seen that points A1 and B1 lie under the stability lobe curve, the vibration displacement in Figure 4 can prove this. And for points C1 and D1, the vibration displacements in Figure 4 can prove that they should lie above the curve of the stability lobe, which is verified in Figure 3. Therefore, the method in this article can predict the stability lobes with high accuracy.

Figure 4.

Simulation of vibration displacement in x-axis direction.

Stability lobes prediction for variable pitch cutter with zero helix angle

The dynamic system with variable pitch angle cutter can suppress the self-excited vibration and improve the limit of the dynamic system. The dynamic system in Altintas et al.’s ²³ paper is used to verify the validity of the proposed method in this article. The parameters of the dynamic system are illustrated in Table 2.

Table 2.

Parameters of the dynamic system.

Direction		Natural frequency (Hz)	Modal mass (kg)	Damping ratios
X	1	441.64	11.125	0.028722
	2	563.6	1.4986	0.055801
	3	778.6	13.063	0.058996
Y	1	516.21	1.199	0.025004

The ratio between radial cutting depth and cutter diameter is 0.5, the pitch angle of cutter is 70°-110°-70°-110°, the number of cutter edge is N = 4, the radius of the cutter is R = 9.525 mm, and the cutting force coefficients are K_t = 679 MPa and K_r = 249.193 MPa. Then the stability lobes are predicted using the Altintas’ method and the proposed method. The computational results are shown in Figure 5.

Figure 5.

Comparison of the stability lobes predicted by literatures’ method and the proposed method.

Four points {A2 (6500 r/min, 3 mm), B2 (8200 r/min, 4 mm), C2(8200 r/min, 4.5 mm) and D2(8000 r/min, 5.5 mm)} are chosen to analyze the predicting accuracy for both methods. It can be found in Figure 5 that the meanings of these points are different from each other; the first point A2 (6500 r/min, 3 mm) represents the unstable state for stability lobes predicted by Altintas but stable state for that predict utilizing the method in this article; the second point B2 lies under these two stability lobe curves; the third point C2 is stable for our stability lobe but unstable for Altintas et al.’s ²³ stability lobe and the fourth point D2 lies above two stability lobe curves. Then, these points are simulated to compute the vibration displacement using the method mentioned above.

The vibration displacements are simulated and illustrated in Figure 6. It can be seen that the displacements for points A2 and D2 are divergent, so the states of these two points are unstable, and on the contrary, the states of other two points are stable. Therefore, the predicting accuracy of the method in this article is higher.

Figure 6.

Simulation of vibration displacement in x-axis direction for the points in Figure 5.

Uniform pitch cutter with non-zero helix angle

In order to verify the validity of the presented modal, a system with the non-zero helix angle in Wan et al.’s ⁴⁹ literature is chosen. The end mill with three cutting edge, helix angle of 30° and the diameter of 16 mm are selected to mill the aluminum alloy 7050. The parameters of the dynamic system are shown in Table 3.

Table 3.

Parameters of the dynamic system.

Directions	Natural frequency (Hz)	Damping ratios	Modal mass (kg)
X	898.22	0.040041	1.576
Y	852.51	0.036768	0.852

The cutting force coefficients along tangent and radial direction are as follows: K_t = 1209.355 MPa, K_r = 501.095 MPa. Down-milling process is conducted, and the radial cutting depth is 5 mm. The axial cutting depth is divided into 60 small elements along the cutting axis. Then, the stability lobe using the model is predicted by the method in this article, and this stability lobe is compared with that in Wan et al.’s ⁴⁹ paper. The results are verified in his paper by experiments and shown in Figure 7, it can be seen that the stability lobe predicted using the model in this article is consist with stability lobes predicted by Wan. Therefore, the model in this article is valid for the cutter with the non-zero helix angle.

Figure 7.

Stability lobes of the dynamic system with the cutter: helix angle of 30°, three cutting edge and the diameter of 16 mm.

Experimental verification

An experiment is conducted using thin-wall part which is taken as 1-DOF dynamic system, for the reason that the stiffness in thickness direction is more smaller than that in the other directions. The parameters of the dynamic system can be seen in Table 4.

Table 4.

Parameters of the dynamic system.

Directions	Natural frequency (Hz)	Damping ratios	Modal stiffness (N/m)
Y	950.3	0.017002	4.4e6

The diameter of the cutter is 12 mm, the helix angle of the cutter is 45°, the cutting force coefficients are K_t = 1625e6 MPa and K_r = 533e6 MPa, and the ratio between radial cutting depth and cutter diameter is 0.02. Down-milling process is adopted.

The results are shown in Figure 8. In this experiment, the parameter combinations of spindle speed and axial cutting depth are given as (2500 r/min, 6 mm), (2500 r/min, 12 mm), (3250 r/min, 6 mm), (3250 r/min, 12 mm), (4000 r/min, 6 mm), (4000 r/min, 12 mm), (4750 r/min, 6 mm) and (4750 r/min, 12 mm). The experimental results prove that the process parameters under the stability limit curve, namely (4750 r/min, 6 mm) and (4750 r/min, 12 mm), are stable, and the parameter combinations over the stability limit curve, namely (2500 r/min, 6 mm), (2500 r/min, 12 mm), (3250 r/min, 12 m), (4000 r/min, 6 mm) and (4000 r/min,12 m), are also unstable in real machining. The point (3250 r/min, 6 mm) is stable, whereas this point is above the stability lobes. The reason for this may be this point is too close to the stability curve and these points that close to the curve are influenced by the environment more easily. These show that the proposed method is of valid and effective.

Figure 8.

Experimental verification of the mode.

Parameters’ effect on the stability lobes

Utilizing the modal in this article, the parameters such as up or down-milling, helix angle and radial cutting depth’s effect on the stability lobes are studied. Here, the classical example using the variable pitch angle {70°-110°-70°-110°} cutter in section “Stability lobes prediction for variable pitch cutter with zero helix angle” is taken as the example to analyze the effect of these parameters on the stability lobes.

Up-milling and down-milling

In this section, the effect of up- or down-milling’s effect on the stability lobes is studied using some parameters. And three group parameters concluding the cutter’s helix angle β and the ratio between the radial cutting depth and the cutter diameter aD are chosen, these parameters are (a) {β = 45°, aD = 0.4}, (b) {β = 30°, aD = 0.6} and (c): {β = 0°, aD = 0.5}. Then, three graphs are plotted using these three group parameters in Figure 9(a)–(c).

Figure 9.

Up- or down-milling’s effect on the stability lobes: (a) helix angle is 45°, aD = 0.4; (b) helix angle is 30°, aD = 0.6 and (c) helix angle is 0°, aD = 0.5.

Some phenomenon can be found in Figure 9; they are as follows:

1. Different milling direction (up- or down-milling) has great influence on the stability lobes. For all the lobes in (a), (b) or (c) of Figure 9, the stability lobes are different from each other.

2. Up-milling process is more helpful for the dynamic system’s stability. For all figures in (a), (b) and (c) of Figure 9, the stability lobes of up-milling process are higher than that of the down-milling process.

The cutting force for the cutter edge of up-milling process changes from zero to maximum during the machining process, whereas the cutting force for the cutter edge of down-milling process changes from maximum to zero. Therefore, up-milling process has smaller impact on the machining system compared with down-milling process for the above two reasons. Therefore, up-milling process has smaller effect on the machining system compared with down-milling process. This is probably the reason why up-milling process is more helpful for keeping the system stable.

Helix angle

Three group parameters are chosen as: aD = 0.1 with β = 0°, 30°, 45° and 60°; aD = 0.5 with β = 0°, 30°, 45° and 60° and aD = 0.9 with β = 0°, 30°, 45° and 60°. Then, the stability lobes with these parameters are predicted using the model and the method in this article. And these stability lobes’ curves are illustrated in Figure 10.

Figure 10.

Influence of the helix angle on the stability lobes: (a) aD = 0.1, β = 0°, 30°, 45° and 60°; (b) aD = 0.5, β = 0°, 30°, 45° and 60° and (c) aD = 0.9, β = 0°, 30°, 45° and 60°.

The first phenomenon in Figure 10 is that the helix angle has effect on the stability lobes, and the possible reason for this maybe the helix angle induces the change of the instantaneous uncut chip thickness for ith cutter edge, whereas when helix angle is 0, the instantaneous uncut chip thickness is fixed for ith cutter edge.

The second phenomenon can be found that when helix angle is larger than 0° and smaller than 30°, helix angle has obvious influence on the stability of the dynamic system. However, when helix angle is larger than 30°, the helix angle’s effect on the stability lobe is small. The third phenomenon is that the larger of the aD, the smaller effect of the helix angle on the stability lobes.

Radial cutting depth

In this section, different cutting depths are chosen to analyze the effect of the axial cutting depth on the stability lobes. Here, the parameters are chosen as follows: aD = 0.5, β = 30°; aD = 0.6, β = 30°; aD = 0.7, β = 30°; aD = 0.8, β = 30°; aD = 0.9, β = 30° and aD = 1.0, β = 30°. Then, the stability lobes are predicted using the presented method and illustrated using three-dimensional picture in Figure 11.

Figure 11.

Effect of the radial cutting depth on the stability lobes.

It can be found that as the radial cutting depth increases, the limit of axial cutting depth decreases at the same spindle rotational speed. As the radial cutting depth increases, the cutting forces increase, and the cutting force coefficients change simultaneously. This may be the reason for the decrease of the axial cutting depth limit as the increase of radial cutting depth.

In order to investigate the radial cutting depth’s influences on the axial cutting depth limit, the relationship curves between the radial cutting depths given as form of aD from 0.5 to 0.9 and the axial cutting depth limit at different spindle speeds are studied. These curves are illustrated in Figure 12. From this figure, a truth can be found that the relationship between axial cutting depth limit and aD cannot be expressed by a line. The reason for this may be that when radial cutting depth increases, the cutting forces and the cutting force coefficients change nonlinearly.

Figure 12.

Relationship between radial cutting depth and the axial cutting depth.

Conclusion

A mathematical model of the dynamic system combining the effect of the cutter helix angle is deduced in this article, and then, a new integral interpolation method using this model to predict the stability lobe is introduced. Subsequently, the properties of this method are analyzed. And three examples are used to verify the validity. Utilizing the method, some key parameters’ effects on the stability lobes are analyzed. Finally, the following conclusions can be obtained:

1. The model of the dynamic system in this article can combine the effect of the cutter’s helix angle. The cutter is divided into a series of cutting elements along the cuter axis. Each cutting element in machining process is taken as a dynamic system; then combining all the elements’ control equation, the model of the cutter’s dynamic system can be obtained.

2. The prediction method of stability lobes in this article can be used in the dynamic system with the multi-regenerative chatter. This article takes the multi-regenerative chatter as a focus to study, the mathematical model which can be used for the multi-regenerative delays is built.

3. The influence of the helix angle on the stability of the dynamic system is studied. A phenomenon can be found in the example that when the helix angle is larger than 0 and smaller than the 30°, the effect of the helix angle on the stability of the dynamic system is obvious. However, when the helix angle is larger than 30°, the helix angle has little influence on the stability of the dynamic system.

4. The milling direction of the machining process has great effect on the stability lobes. Especially, the up-milling process is more helpful for keeping the system stable.

5. The three-dimensional stability lobes are computed for special cutting condition. It can be found that as the radial cutting depth increases, the limit of axial cutting depth decreases at the same spindle rotational speed.