Stability Analysis of Milling Process with Multiple Delays

Chatter is a self-excited vibration induced in the machining process which reduces the machining efficiency and surface quality, accelerates the tool wear and shortens the tool durability. There have been a lot of research on the mechanism of chatter [1,2,3,4,5,6], among which regenerative effect was usually considered as a main inducement. Since the regenerative chatter is caused by the phase differences between wavy surfaces left by the adjacent teeth, it can be suppressed by adjusting the pitch angle between the adjacent teeth. The use of variable pitch cutters to improve the stability of the milling process is built on this idea. Different from the uniform pitch cutter, when the variable pitch cutter is used, the dynamics model of cutting vibration changes from DDEs with single delay to DDEs with multiple delays.

The stability analysis of the milling process with multiple delays has been studied through different methods. Altintas et al. [7] adopted the frequency domain method to analyze the milling stability of the variable pitch cutter and presented a method to select the optimal pitch angles. Budak [8,9] proposed an analytical design method for nonconstant pitch milling cutters to improve the stability of the milling process and verified it through cutting experiments. The results showed that chatter stability can be improved significantly even at slow cutting speeds by properly designing the pitch angles. The modified semi-discretization method (SDM) was used by Sims et al. [10] to investigate the stability of variable pitch and variable helix milling tools. Considering the tool runout and the pitch angle, Wan et al. [11] proposed a unified stability prediction method for a milling process with multiple delays based on the SDM. Zhang et al. proposed an improved full-discretization method (FDM) [12] and the variable-step numerical integration method [13] for milling chatter stability prediction with multiple delays. Olgac and Sipahi [14] studied the stability boundary in a milling process with variable-pitch cutters based on the cluster treatment of characteristic roots method and presented an optimization procedure for the geometry design of variable pitch cutters. The enhanced multistage homotopy perturbation method was presented by Compean et al. [15] to analyze the milling stability influenced by the variable cutting edge pitch, the variable helix, and the variable rake angle. The spectral element approach was first introduced by Khasawneh et al. to study the stability of delay systems with single delay [16]. Later, it was generalized to time-periodic delay equations with multiple delays [17]. Lehotzky et al. [18] generalized the spectral element method for time-periodic DDEs with multiple delays and distributed delay and given an explicit formula for the construction of the matrix approximation of the monodromy operator. An improved SDM was presented by Jin et al. [19,20] to predict the stability of milling process with multiple delays. On the basis of the multifrequency solution, Andreas et al. [21] presented a new method for the identification of the chatter stability lobes in milling with non-uniform pitch and variable helix tools. Tao et al. [22] proposed an extended Adams–Moulton-based method for the stability prediction of milling processes with multiple delays.

After the stability analysis is completed, the stability limits under different processing parameters are usually displayed in one diagram, which is known as the SLD. The SLD is a powerful tool, which can not only help to select the optimal process parameters to realize the chatter-free machining, but also assist the geometric design of the cutting tools. A variety of methods have been proposed to improve the accuracy or calculation efficiency of the SLD. The improved FDM was used in [12] to reduce the calculation time. However, similar to the SDM, it was also based on the equal-step discretization technique. In the case of low radical immersion cutting, it is necessary to reduce the time step to maintain the accuracy. Besides, it may not be able to guarantee that all the discretization number of the forced vibration intervals are integers in the case of multiple delays. The variable-step method proposed in [13] can discretize the forced vibration intervals accurately, and thus the calculation result has high accuracy, however, the effect of the helix angle is not considered in the derivation process. Ozoegwu et al. [23] proposed a speed-dependent discretization method to improve the calculation accuracy at lower spindle speed and save the calculation time. Li et al. [24] applied an iterative dichotomy search algorithm to determine the stability limit of the milling process. The dichotomy search algorithm has a high efficiency, but it is not suitable for the situation where the unstable islands exit [25].

It has always been desirable to improve the accuracy and calculation efficiency of stability analysis. In this work, an AVSNIM, considering the helix angle, and a novel spindle speed-dependent discretization algorithm, are developed and combined to analyze the stability of the milling process with multiple time delays. They are expected to obtain high accuracy and shorten the calculation time. The rest of this paper is organized as follows: the dynamics model of the milling process with multiple delays is outlined in Section 2. In Section 3, an AVSNIM that takes into account the helix angle is developed and verified by a benchmark milling model. In Section 4, a novel spindle speed-dependent discretization algorithm is proposed and combined with the AVSNIM to further reduce the calculation time of the SLD. Finally, conclusions are derived in Section 5.

Milling vibration considering the regeneration effect is usually modeled by DDEs which can be written as where $M$, $C$ and $K$ are the mass, damping and stiffness matrices respectively, $y\left(t\right)$ is the vibration displacements vector and $f\left(t\right)$ is the time-varying cutting force vector. When using the variable pitch cutter (as described in Figure 1), the dynamics model of the milling process becomes DDEs with multiple delays. The cutting force $f\left(t\right)$ in Equation (1) can therefore be expressed as where $N$ is the number of teeth, $T_j$ is the time delay of the j-th tooth, $y\left(t\right)$ and $y\left(t-T_j\right)$ are the vibration displacements at the current moment and the previous tooth-passing period of the j-th tooth separately. $H_j\left(t\right)$ is the time-varying cutting force coefficient matrix of the j-th tooth, which satisfies $H_j\left(t\right)=H_j\left(t+T\right)$, where $T$ is the spindle period. $H_j\left(t\right)$ can be written as where where $K_t$ and $K_n$ are the tangential and radial cutting force coefficients, respectively. $z_{u,}{}_j\left(t\right)$ and $z_l{}_{,j}\left(t\right)$ are the upper and lower bounds of the cutting edge participating in cutting on the j-th tooth as shown in Figure 1c, respectively. ${\phi }_j\left(t,z\right)$ is the angular position of the differential element of cutting edge with a height $z$ on the j-th tooth, which can be expressed as where ${\phi }_j\left(0,0\right)$ is the angular position of the lowest differential element on the j-th tooth at the initial time, $\Omega$ is the spindle speed, ${\varphi }_{lag}\left(z\right)$ is the lag angle of the differential element at height $z$ as shown in Figure 1d, and ${\varphi }_{lag}\left(z\right)=\left(2\tan \beta /D\right)\cdot z$, where $D$ is the diameter of the cutter, $\beta$ is the helix angle.

Equation (1) can be converted to the state-space form as where $A$ is a constant matrix, representing the time-invariant characteristics of the milling system, while $B_j$ is a periodic matrix satisfying $B_j\left(t\right)=B_j\left(t+T\right)$, which describes the periodic characteristics of the system. $x$ is known as the state vector where $I$ is an identity matrix with the same size as $M$.

Obviously, Equation (6) describes a linear periodic system. According to the Floquet theory, the stability of a linear periodic system depends on the spectral radius of the Floquet transition matrix: if the spectral radius is less than one, then the system is asymptotically stable; if it is larger than one, the system is unstable; if the spectral radius equals to one, then the system is in a critical stable state.

The variable-step technique can discrete the milling period accurately. Considering the effect of helix angle, an AVSNIM for the stability analysis of milling process with multiple delays is proposed in this section.

The AVSNIM can provide accurate stability limits, the accuracy of which depends on the pre-set discretization parameter $m_p$. A large discretization parameter usually brings a small time step and high calculation accuracy, and at the same time, increases the calculation time. In most of the research work, the SLD is obtained under constant discretization parameter. A constant discretization parameter will lead to different time steps and calculation accuracies under different spindle speeds. Generally, the accuracy of the stability limits in the region with low spindle speed is not as high as that in the region with a high spindle speed in the SLD. This implies that if the discretization parameter can be properly selected at each different spindle speed (usually, a large discretization parameter is used in the low-speed region and a small discretization parameter is used in the high-speed region), it is possible to improve the calculation accuracy of stability limit at a low speed and save the calculation time. Inspired by this, a spindle speed-dependent discretization algorithm is proposed to obtain the appropriate discretization parameters, and used in conjunction with the AVSNIM. The process of the algorithm is listed as follows:

In the above algorithm, parameters $\Delta \Omega$, $\epsilon$ and $\Delta m$ are variable input parameters. The value of $\Delta \Omega$ determines the density of spindle speed sampling points. The smaller the value of $\Delta \Omega$, the denser the sampling nodes, and the higher the interpolation accuracy of discretization parameters. The value of $\epsilon$ sets a relative error. The smaller the relative error is, the closer the stability limit is to the reference value. $\Delta m$ is the increment step of the discretization parameter, the value of which determines the speed of the discretization parameter approaching the reference discretization parameter. If $\Delta m$ is too small, the number of iterations will increase. On the contrary, if $\Delta m$ is too large, a discretization parameter larger than the required value may be obtained, thereby increasing the final calculation load. The values of these three parameters can be selected according to the actual engineering requirements.

Numerical simulation is performed to verify the above algorithm; two milling models are studied in this section.

The first milling model is the same as used in Section 3.2. The simulation parameters are: down milling, $\epsilon =1\times e^{-3}$, $\Delta m=2$, $m_p=200$, the initial $m_d=40$, the spindle speed $\Omega$ ranges from 2500 to 12,500 rpm and the axial depth of cut $w$ ranges from 0 to 15 mm.

Firstly, the time step under constant discretization parameter and the discretization parameters based on the spindle-speed-dependent algorithm are simulated. The results are shown in Figure 8. Figure 8a,b describes the time step under the equal-step method and the variable-step method, respectively. It can be found that the time step in the equal-step algorithm is only related to the spindle speed, while in the variable-step method, due to the effect of the helix angle, the step size depends on both the spindle speed and the axial depth of cut. Figure 8c shows the discretization parameters obtained through the spindle-speed-dependent algorithm. It can be seen that under the given relative error $\epsilon$, the discretization parameters obtained by the spindle-speed-dependent algorithm are smaller than the reference discretization parameter, which means that the time required to obtain the SLD will be effectively reduced. At the same time, we can see that the value of the discretization parameter is higher in the low-spindle-speed region than that in the high-spindle-speed region in general, which proves the previous conclusion that, in the low-speed region, in order to maintain the accuracy of the numerical calculation, larger discretization parameters are required. However, the discretization parameters do not monotonically decrease with the increase in spindle speed, but there are some fluctuations. There may be two incentives for these fluctuations. One is that they may be caused by numerical calculation, because the bisection method used in this article may randomly converge to the upper or lower boundary of the relative error limit, which usually causes some numerical fluctuations. The second is that they may be caused by the different sensitivity of the eigenvalues of Floquet transition matrix to the processing parameters (mainly the spindle speed here).

Subsequently, the SLD is provided to verify the effectiveness of the proposed spindle-speed-dependent discretization algorithm. Figure 9 shows the SLD obtained through the algorithm proposed above; the SLD is constructed over a 200 × 150-sized grid; the radial immersion in Figure 9a–c is selected as 1.0, 0.6 and 0.1, respectively. The reference stability limits demoed by the red line are obtained by the AVSNIM with a large reference discretization parameter. From Figure 9a–c, we can see that the stability limits obtained by the proposed algorithm are very close to the reference stability limits. This shows that the proposed algorithm is reliable in terms of calculation accuracy. Meanwhile, since the radial immersion in Figure 9 ranges from 1.0 to 0.1, this indicates that the proposed algorithm is applicable to all kinds of radial immersion.

Finally, the calculation time required for the AVSNIM with spindle speed-dependent discretization parameters is listed in Table 2 to compare with the AVSNIM with constant discretization parameter. As outlined in Table 2, the AVSNIM with spindle-speed-dependent discretization parameters can significantly reduce the calculation time.

Another milling model to verify the proposed algorithm is from the literature [31]. In this model, a 12.75 mm diameter cutter with one flute is adopted, the helix angle of which is 30°. The detailed values of the modal parameters and cutting force coefficients are as follows: the modal masses in x and y directions are $m_x=m_y=0.23\mathrm{kg}$; the modal damping coefficients in x and y directions are $c_x=c_y=15.73\mathrm{Ns}/\mathrm{m}$; the modal stiffness in x and y directions are $k_x=k_y=8.4\times {10}^5\mathrm{N}/\mathrm{m}$; the cutting force coefficients are $K_t=5.36\times {10}^8{\mathrm{N}/\mathrm{m}}^2$ and $K_n=1.87\times {10}^8{\mathrm{N}/\mathrm{m}}^2$. The simulation parameters are: down milling, $\epsilon =1\times e^{-3}$, $\Delta m=2$, $m_p=60$, the initial $m_d=10$, the spindle speed $\Omega$ ranges from 5000 to 30,000 rpm and the axial depth of cut $w$ ranges from 0 to 25 mm. The simulation result is shown in Figure 10. The SLD in Figure 10 is constructed over a 500 × 250-sized grid; the radial immersion is selected as 0.02. The reference stability limits demoed by the red line are obtained by the AVSNIM with a constant discretization parameter $m=100$.

As shown in Figure 10, the stability boundary is very complex, in which there are some islands. The stability limits obtained by the proposed algorithm are consistent with the reference stability limits and are the same as those obtained in the literature [31]. This shows that the algorithm proposed in this work is still effective when the stability boundary is very complex. It is worth noting that there is only a single time delay in this milling model, which is the same as the milling process with a uniform pitch milling cutter. This proves that the proposed algorithm is also suitable for the milling process with a uniform pitch milling cutter.

Based on the above analysis, we can conclude that the proposed algorithm can not only maintain the calculation accuracy but also effectively reduce the calculation time, and hence can meet the needs in the engineering field well.

This work focuses on the stability analysis of milling process with multiple delays. An AVSNIM considering the helix angle is developed. According to the geometric characteristics of the cutter (pitch helix angle, etc.) and machining parameters (spindle speed and axial cutting depth), one spindle period is divided into several segments; each segment contains a free vibration time interval and a forced vibration time interval. Then, the discretization of one spindle period is completed adaptively according to the pre-set discretization parameter. Numerical integration algorithm and Lagrange interpolation method are used to obtain the vibration state vector on each time node, and the Floquet transition matrix between two adjacent spindle periods is constructed accordingly. Finally, the stability of milling process can be analyzed on the basis of the Floquet theory. Numerical verification is conducted with the two-DOF benchmark milling model, the results of which validate that the proposed method has superior accuracy. Subsequently, a novel spindle-speed-dependent discretization algorithm is proposed to further improve the calculation efficiency of the SLD. Numerical simulation shows that this method can effectively shorten the calculation time of the SLD. The advantages of the methods proposed in this work are as follows: