Chatter stability analysis of variable speed milling with helix angled cutters

Abstract

The number of factors considered in the establishing of chatter model for regenerative chatter would directly affect the accuracy of numerical analysis. In this paper, considering the helix angle of cutter, the establishment method of three degrees of freedom complex model is investigated synthetically. An improved semi-discretization method is adopted for the prediction of milling stability of variable speed milling with helix angle. The influences of milling cutter helix angle, radial immersion rate and variable speed modulation coefficient on chatter suppression in milling process are studied. The results show that the selection of milling cutters with different helix angles has a certain influence on the overall chatter stability, especially the effect cannot be ignored under local special machining conditions. Excessive radial immersion rate will increase the probability of chatter. Compared with constant speed milling, the chatter can be suppressed obviously by variable speed milling, and the selection of variable speed modulation coefficient will directly affect the suppression effect of chatter.

Introduction

During the milling process, the tool is in the state of high-speed rotation, the cutting force and the actual cutting thickness of the tool are always changing dynamically.^1,2 Regenerative chatter is caused by the difference of the phase between the chatter marks formed in the last cutting and the vibration displacement of this cutting, which results in the different cutting thickness of the tool.^3–6 The occurrence of milling chatter will affect the surface quality of work-piece, and even lead to tool wear. In order to avoid chatter effectively and not affecting the efficiency of actual processing, it is very important to select the normal milling parameters under the current processing conditions. ⁷ Therefore, chatter stability analysis is the key to predict milling chatter and suppress chatter effectively.^8,9

Research on stability of milling chatter has never stopped. Altintas and Budak ¹⁰ put forward a theory of stability lobes based on the average time period dynamic milling coefficient, that is, the corresponding cutting depth without chatter could be selected for different rotating speeds, which laid a theoretical foundation for the later study of chatter stability analysis. Inserger and Stépán ¹¹ put forward a semi-discretization method that could deal with time-delay terms, using a series of ordinary differential equations with known solutions to replace the original time-delay differential equations for numerical calculation. This method improved the accuracy of low-speed stability analysis and was widely used in the stability prediction of various milling systems. Ding et al.^12,13 proposed the total discretization method based on numerical integration and direct integration method. Compean et al. ¹⁴ innovatively developed an enhanced multi-level homotopic perturbation method for chatter stability prediction. Ding et al. ¹⁵ and Huang et al. ¹⁶ used a simple and efficient semi-analytical time domain method for stability analysis, but this method was based on linear acceleration approximation.

Other researchers are also focused on using different ways to improve the accuracy and convergence. Li et al. ¹⁷ proposed a zero order semi-discretization stability analysis method based on Magnus Gaussian truncation. Niu et al. ¹⁸ proposed the analysis method based on Runge Kutta, etc. Under the same discrete time interval, as the matrix exponential function of the full discrete method only depends on the rotation speed, the calculation time was shorter, but the convergence speed was lower than that of the first-order method. ¹⁹ The total discrete method can also be used in the stability classification of multiple time-delay terms, and could further improve the convergence.^20,21 Li et al. ²² put forward a new method which was different from the traditional full discrete method. It discretized all time related terms, and obtained the iterative formula by numerical iterative method, thus obtained the Floquet transfer matrix. Xie and Zhang ²³ proposed an improved semi-discretization method for milling chatter stability prediction. According to the two-dimensional cutting force model, the stable lobes of variable speed milling and constant speed milling were compared, which proved that the chatter could be suppressed by changing the spindle speed under the current assumption model. Altintas and Ko ²⁴ considered the influence of vibration in two transverse, one axial and torsional directions on the chatter stability. It was found that the model can be applied to the chatter stability analysis and prediction by changing the meshing conditions between component and work-pieces. Long et al. ²⁵ proposed a four degree of freedom chatter model based on cutter and work-piece, and compared constant speed milling with variable speed milling, and found that within a certain range of spindle speed, the change of spindle speed was beneficial. Niu et al. ²⁶ proposed a general method based on VNIM for stability study and modulation parameter optimization of variable axis speed milling process. By comparing different modulation schemes in a unified framework, it was found that sinusoidal modulation can achieve high processing efficiency, while triangular modulation was allowed for greater selection of modulation coefficients. Seguy et al. ²⁷ also studied the influence of variable speed on the chatter of high-speed milling. The influence of frequency and amplitude on the stability was analyzed by parameterization, and also confirmed that the change of spindle speed could inhibit the vibration. Zatarain et al. ²⁸ found that the helix angle is very important to the study of chatter stability, and found that when the cutting depth equaled to the multiple of helix angle pitch, the FILP lobes were closed and separated by horizontal lines. Masud Akhtar et al. ²⁹ studied the influence of helix angle on down milling, up milling and gear shaping. A multi frequency method was proposed to solve the milling force with the influence of helix angle. The results showed that the stable leaf structure with small helix angle followed the Altintas (MFS) model, and the stable leaf structure with large helix angle followed the Altintas (ZOA) model. Song et al. ³⁰ established the dynamic balance equation of milling process stability considering the helix structure, and used the finite dimension matrix to determine the stability of milling. It was found that the helix angle was important in the stability of milling process. Jin et al. ³¹ adopted a method to obtain the three degrees of freedom stable lobes of thin-walled milling, aiming at the helix angle effect of tools and the dynamic characteristics of thin-walled component. The experimental results also showed that when considering the model established by the helix angle, the stability analysis would have different results. Compared with that without considering the helix angle effect, the stability limit curve with the helix angle effect was more accurate. Wang et al. ³² used a new calculation method of milling force considering helix angle and bending, and improved the calculation method of instantaneous cutting area. It was found that tool bending and helix were crucial in stabilizing lobes. Chong et al. ³³ proposed an implicit four-order and three-step ADAMS numerical analysis method to solve the milling dynamic differential equation, which provided the optimization method of cutting parameters for circular corner milling. Chen et al. ³⁴ employed a model to predict the surface morphology of peripheral milling based on Z-map model and dynamic cutting force model of peripheral milling considering regeneration effect. It was found that when the frequency of the milling system was an integral multiple of the machining frequency, a stable cutting state coul be obtained but the surface finish was poor. A stable cutting state could be obtained on the left side of the resonance spindle speed and a high surface gloss was guaranteed. Li et al. ³⁵ established the milling model of process damping and mode coupling based on the 2-DOF model of regenerative chatter. Modal coupling can increase the area of the stable region, and process damping has a greater influence on the stability of the low-speed region.

Although there are many researches on milling chatter stability, the current research situation is limited to the establishment of chatter model under ideal conditions. ³⁶ Too many assumptions and simplifications on the model may affect the accuracy of stability analysis.^1,10,31 And the stability analysis of thin-walled milling is still not well enough. In the processing of thin-walled component, the cutting amount has been optimized compared with the traditional milling process, but the actual thickness of the part processed is too small, and the ball-end milling cutter is often used in thin-walled processing, so the regenerative chatter caused by the axial displacement in the processing process has to be considered. In the previous study of chatter stability, in order to facilitate the analysis, the model is usually simplified, and the helix angle of cutter is defaulted to zero. However, in the actual complex machining process, this assumption has not met the requirements. The introduction of helix angle into the contact angle can make the chatter stability analysis more suitable for the actual situation. some researchers have shown that the method of variable spindle speed can suppress the occurrence of chatter effectively, but these research is either based on the low degree of freedom model, or not considering enough actual factors to model. Therefore, it will directly affect the evaluation of variable speed in the actual complex machining process. In this paper, it is assumed that the tool is flexible relative to the rigid work-piece. According to the helical angle effect of cutter and the dynamic characteristics of thin-wall component to establish three degrees of freedom chatter model of thin-walled milling. The actual milling state is divided into five kinds based on the geometrical shape and machining principle of milling tool, so as to calculate the cutting force more accurately. At the same time, the chatter suppression effect of variable speed in three-dimensional model is also studied.

Three degrees of freedom milling chatter model

A variable speed milling model with sinusoidal modulation is established. The dynamic milling process is expressed by a periodic system with variable time delay. Using an improved semi-discretization method to analyzed the stability of variable speed milling considering the helix angle. In the establishment of the cutting force model, the cutting area is divided into five cases to discuss, so as to avoid the unreasonable situation that the forces on the cutter teeth are equal everywhere. Through the establishment of the above chatter model, the stability critical curve in line with the actual machining state can be obtained by digital simulation, so as to provide the parameter selection of axial cutting depth and spindle speed before machining.

Radial contact angle and helix angle

The helix cutters are widely used in the machining of thin-walled component because they can obtain higher precision surface quality. Therefore, when modeling the cutting contact angle, besides considering the number of milling cutter teeth, it is also necessary to consider the milling cutter helix angle. So the cutting radial contact angle can be expressed as:

$${\varphi }_j(t)=\frac{2\pi }{60}{\int }_0^t\mathrm{\Omega }(x)\mathrm{dx}+\frac{2\pi }{N}(j-1)+\frac{\mathrm{da}\tan \beta }{R}(l-1)$$

(1)

where ${\varphi }_j(t)$ is also the instantaneous angular immersion of tooth (j) measured clockwise from the normal (Y) axis, ⁹ as shown in Figure 1, Ω is the spindle speed(r/min), N is the number of milling cutter teeth, j is the j-th cutter tooth, a is the axial cutting thickness (da is the micro element thickness in the axial direction), β is the helix angle of cutter, R is the milling cutter radius, and l is the l-th micro unit from bottom to top in the axial direction.

Figure 1.

Diagram of 2D dynamic milling system.

The traditional machining methods usually set the spindle speed as a fixed value to machine the components. In order to study the effect of variable speed and constant speed on chatter suppression in three degrees of freedom complex model, this paper adds sine modulation to make the spindle speed change periodically in the process of machining. The spindle speed is expressed as:

$$\mathrm{\Omega }(t)=n_0+n_1\sin ({\omega }_{\mathrm{\Omega }}t)$$

(2)

where $n_0$ is the rated speed, $n_1$ is the amplitude of speed change, and ${\omega }_{\mathrm{\Omega }}$ is the frequency of speed change.

The spindle speed cannot be adjusted arbitrarily, it needs to be adjusted from the two dimensions of amplitude and frequency based on the rated speed, that is, the introduction of two modulation coefficients $\gamma$ and $\rho$ . The specific calculation formula is:

$$\{\begin{array}{c}{n}_1=\gamma n_0\hfill \\ \multicolumn{1}{c}{{\omega }_{\mathrm{\Omega }}=\rho (\pi n_0/30)}\end{array}$$

(3)

Take equations (2) and (3) into equation (1), then the contact of angle the j-th cutter tooth at time t can be more accurately expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{{\varphi }_j(t)=\frac{\pi }{30}(n_{0}t+\frac{30\gamma }{\pi \rho }(1-\cos ({\omega }_{\mathrm{\Omega }}t)))}\\ +\frac{2\pi }{N}(j-1)+\frac{\mathrm{da}\tan \beta }{R}(l-1)\hfill \end{array}$$

(4)

The position of the milling cutter head is regarded as the coordinate with zero cutting depth. When the helix angle of the milling cutter is not considered, there is no angle difference between the cutting points of the same tooth at depth 0 and depth a. As shown in Figure 2, $o\left({\varphi }_{\mathrm{st}}\right)$ and $o\left({\varphi }_{\mathrm{ex}}\right)$ are the position points corresponding to the cut-in angle and cut-out angle, respectively. So the differential units of the same tooth at different cutting positions can use the same formula to calculate the cutting force.

Figure 2.

Cutting state without helix angle.

In the actual machining, the milling cutter mostly has a helix angle. In order to comprehensively analyze the influence of helix angle under different cutting states, taking a certain cutter tooth of milling cutter as an example, the actual cutting situation is divided into five types to discuss. When calculating the cutting force of each differential unit, firstly determine which of the five cutting states is currently in operation, and then calculate the cutting force of the differential unit, rather than assuming that the cutting force of the tooth is equal everywhere. As shown in Figure 3. The first case is that the lower half of the current tooth is in the cutting state, expressed as:

$$\{\begin{array}{c}{{\varphi }_j|}_{z=0}\in ({\varphi }_{st},{\varphi }_{ex})\\ {{\varphi }_j|}_{z=a}\in (0,{\varphi }_{st})\end{array}\Rightarrow \{\begin{array}{c}{o}_1=0\\ o_2=\left({{\varphi }_j|}_{z=a}-{\varphi }_{st}\right)R/\tan \beta \end{array}$$

(5)

Figure 3.

Five cutting states considering the helix angle.

The second case is that the current tooth is completely in the cutting state, expressed as:

$$\{\begin{array}{c}{{\varphi }_j|}_{z=0}\in ({\varphi }_{\mathrm{st}},{\varphi }_{\mathrm{ex}})\hfill \\ \multicolumn{1}{c}{{{\varphi }_j|}_{z=a}\in ({\varphi }_{\mathrm{st}},{\varphi }_{\mathrm{ex}})}\end{array}\Rightarrow \{\begin{array}{c}{\mathrm{o}}_1=0\hfill \\ \multicolumn{1}{c}{{\mathrm{o}}_2=a}\end{array}$$

(6)

The third case is that the middle part of the current tooth is in the cutting state, expressed as:

$$\{\begin{array}{c}{{\varphi }_j|}_{z=0}\in ({\varphi }_{ex},\propto )\\ {{\varphi }_j|}_{z=a}\in (0,{\varphi }_{st})\end{array}\Rightarrow \{\begin{array}{c}{o}_1=\left({{\varphi }_j|}_{z=0}-{\varphi }_{ex}\right)R/\tan \beta \\ o_2=\left({{\varphi }_j|}_{z=a}-{\varphi }_{st}\right)R/\tan \beta \end{array}$$

(7)

The fourth case is that the upper half of current tooth is in the cutting state, expressed as:

$$\{\begin{array}{c}{{\varphi }_j|}_{z=0}\in ({\varphi }_{ex},\propto )\\ {{\varphi }_j|}_{z=a}\in ({\varphi }_{st},{\varphi }_{ex})\end{array}\Rightarrow \{\begin{array}{c}{o}_1=\left({{\varphi }_j|}_{z=0}-{\varphi }_{ex}\right)R/\tan \beta \\ {\mathrm{o}}_2=a\end{array}$$

(8)

The fifth case is that the current tooth is completely out of the cutting state, expressed as:

$$\{\begin{array}{c}{{\varphi }_j|}_{z=0}\in ({\varphi }_{\mathrm{ex}},\propto )\hfill \\ \multicolumn{1}{c}{{{\varphi }_j|}_{z=a}\in ({\varphi }_{\mathrm{ex}},\propto )}\end{array}\Rightarrow \mathrm{Break}\mathrm{away}\mathrm{from}\mathrm{cutting}\mathrm{area}$$

(9)

Three degrees of freedom milling system model

As observed in Figures 1 and 4, the differential equation of dynamics of the three degrees of freedom system is described as follows:

$$\{\begin{array}{c}{m}_x\stackrel{··}{x}\left(t\right)+c_x\stackrel{·}{x}\left(t\right)+k_{x}x\left(t\right)=F_x\hfill \\ \multicolumn{1}{c}{m_y\stackrel{··}{y}\left(t\right)+c_y\stackrel{·}{y}\left(t\right)+k_{y}y\left(t\right)=F_y}\\ \multicolumn{1}{c}{m_z\stackrel{··}{z}\left(t\right)+c_z\stackrel{·}{z}\left(t\right)+k_{z}x\left(t\right)=F_z}\end{array}$$

(10)

where $m,c,k$ represent the modal mass, damping, and stiffness of the milling structure, respectively. $x(t),y(t),z(t)$ represent the vibration displacement of the milling structure in the X, Y, Z direction, respectively. $F_x,F_y,F_z$ represent the cutting force of the milling structure in the X, Y, Z direction, respectively.

Figure 4.

Mechanical model of 3-DOF milling system.

The cutting force of milling cutter in three axes of X, Y, Z is expressed as follows:

$$\begin{array}{l}\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)dz}{\sin \delta }\\ \left[\begin{array}{ccc}-\cos {\phi }_j& -\sin {\phi }_j\sin \delta & -\sin {\phi }_j\cos \delta \\ \sin {\phi }_j& -\cos {\phi }_j\sin \delta & -\cos {\phi }_j\cos \delta \\ 0& \cos \delta & -\sin \delta \end{array}\right]\left[\begin{array}{c}{K}_t\\ K_r\\ K_a\end{array}\right]\\ \left[\begin{array}{c}\sin {\phi }_j\sin \delta \\ \cos {\phi }_j\sin \delta \\ -\cos \delta \end{array}\right]\left[\begin{array}{l}x(t)-x(t-\tau (t))\\ y(t)-y(t-\tau (t))\\ z(t)-z(t-\tau (t))\end{array}\right]\end{array}$$

(11)

where $o_1(t)$ and $o_2(t)$ are the actual cut-in point and cut-out point of cutting tool in the z-axis direction, $\delta$ is the axial contact angle, $K_t,K_r,K_a$ represent the tangential, radial, and axial cutting force coefficients, respectively, $\tau (t)$ is the time delay. $g({\varphi }_j)$ is the step function used to judge the cutting state of the cutter teeth, namely:

$$g({\varphi }_j)=\{\begin{array}{c}1{\varphi }_j\in ({\varphi }_{\mathrm{st}},{\varphi }_{\mathrm{ex}})\hfill \\ \multicolumn{1}{c}{0\mathrm{otherwise}}\end{array}$$

(12)

where ${\varphi }_{\mathrm{st}}$ is the starting angle of cutting. ${\varphi }_{\mathrm{ex}}$ is the end angle of cutting. According to the research results of reference, ³⁸ variable speed up milling is more conducive to avoid the occurrence of chatter. This paper adopts the up-milling, so ${\varphi }_{\mathrm{st}}=0$ and ${\varphi }_{\mathrm{ex}}=\arccos (1-2a/D)$ . $a/D$ is the radial immersion rate.

Analysis of improved semi-discretization method

Since Insperger et al.^39,40 proposed semi-discretization methods, various stability analysis methods have been put forward. In this paper, practical complex factors are taken into account to carry out chatter modeling, which not only improves the accuracy of the model, but also increases the difficulty of calculation. Therefore, an improved semi-discretization method is adopted for stability analysis, which ensures the accuracy of stability analysis and greatly optimizes the difficulty of solution.

The first-order expression is obtained by simplifying equation (10), which is expressed as follows:

$$\stackrel{·}{\nu }\left(t\right)=A_0\nu \left(t\right)+A\left(t\right)\nu \left(t\right)-A\left(t\right)\nu (t-\tau \left(t\right))$$

(13)

where

$$\mathbf{\nu }\left(\mathit{t}\right)={\left[\mathit{q}\left(\mathit{t}\right)\mathit{p}\left(\mathit{t}\right)\right]}^T$$

(14)

$$\mathit{q}\left(\mathit{t}\right)=\left\{\begin{array}{c}x\left(t\right)\hfill \\ \multicolumn{1}{c}{y\left(t\right)}\\ \multicolumn{1}{c}{z\left(t\right)}\end{array}\right\}$$

(15)

$$p\left(t\right)=\left\{\begin{array}{c}{m}_x\stackrel{·}{x}\left(t\right)+m_x\zeta {\omega }_{n}x\left(t\right)\hfill \\ \multicolumn{1}{c}{m_y\stackrel{·}{y}\left(t\right)+m_y\zeta {\omega }_{n}y\left(t\right)}\\ \multicolumn{1}{c}{m_z\stackrel{·}{z}\left(t\right)+m_z\zeta {\omega }_{n}z\left(t\right)}\end{array}\right\}$$

(16)

$$\mathit{A}\left(\mathit{t}\right)=\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{K_{\mathrm{xx}}}& K_{\mathrm{xy}}\hfill & K_{\mathrm{xz}}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{K_{\mathrm{yx}}}& K_{\mathrm{yy}}\hfill & K_{\mathrm{yz}}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{K_{\mathrm{zx}}}& K_{\mathrm{zy}}\hfill & K_{\mathrm{zz}}\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(17)

$$\{\begin{array}{c}{K}_{\mathrm{xx}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\sin {\phi }_j\sin \delta (-K_t\cos {\phi }_j-K_r\sin {\phi }_j\sin \delta -K_a\sin {\phi }_j\cos \delta )\hfill \\ K_{\mathrm{xy}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\cos {\phi }_j\sin \delta (-K_t\cos {\phi }_j-K_r\sin {\phi }_j\sin \delta -K_a\sin {\phi }_j\cos \delta )\hfill \\ K_{\mathrm{xz}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\cos \delta (K_t\cos {\phi }_j+K_r\sin {\phi }_j\sin \delta +K_a\sin {\phi }_j\cos \delta )\hfill \\ K_{\mathrm{yx}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\sin {\phi }_j\sin \delta (K_t\sin {\phi }_j-K_r\cos {\phi }_j\sin \delta -K_a\cos {\phi }_j\cos \delta )\hfill \\ K_{\mathrm{yy}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\cos {\phi }_j\sin \delta (K_t\sin {\phi }_j-K_r\cos {\phi }_j\sin \delta -K_a\cos {\phi }_j\cos \delta )\hfill \\ K_{\mathrm{yz}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\cos \delta (-K_t\sin {\phi }_j+K_r\cos {\phi }_j\sin \delta +K_a\cos {\phi }_j\cos \delta )\hfill \\ K_{\mathrm{zx}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\sin {\phi }_j\sin \delta (K_r\cos \delta -K_a\sin \delta )\hfill \\ K_{\mathrm{zy}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\cos {\phi }_j\sin \delta (K_r\cos \delta -K_a\sin \delta )\hfill \\ K_{\mathrm{zz}}=\sum _{j=1}^N\sum _{l=1}^m{\int }_{o_1(t)}^{o_2(t)}\frac{g({\phi }_j)\mathrm{dz}}{\sin \delta }\cos \delta (-K_r\cos \delta +K_a\sin \delta )\hfill \end{array}$$

(18)

Discrete time period $T$ into time intervals $\mathrm{\Delta }t$ , that is, $T=\lambda \mathrm{\Delta }t$ . Let the initial condition be $\mathbf{\nu }\left(t_i\right)={\mathbf{\nu }}_{\mathit{i}}$ . During the period $[t_i,t_{i+1}]$ the differential equation is further calculated:

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{\nu }(t)=e^{{\mathit{A}}_{\mathit{0}}(t-t_i)}{\mathbf{\nu }}_{\mathit{i}}}\\ +{\int }_{t_i}^t\{e^{{\mathit{A}}_{\mathit{0}}(t-\xi )}\mathit{A}(\xi )[\mathbf{\nu }(\xi )-\mathbf{\nu }(\xi -{\tau }_i)]\}d\xi \hfill \end{array}$$

(19)

it can also be equivalent to:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\nu }}_{i+1}=\mathbf{\nu }(t_i+\mathrm{\Delta }t)=e^{{\mathit{A}}_{\mathit{0}}(\mathrm{\Delta }t)}{\mathbf{\nu }}_i}\\ +{\int }_{t_i}^{t_i+\mathrm{\Delta }t}\{e^{{\mathit{A}}_{\mathit{0}}(t+\mathrm{\Delta }t-\xi )}\mathit{A}(\xi )[\mathbf{\nu }(\xi )-\mathbf{\nu }(\xi -{\tau }_i)]\}d\xi \hfill \end{array}$$

(20)

where ${\tau }_i$ is the average delay of discrete interval expressed as following:

$${\tau }_i=\frac{1}{\mathrm{\Delta }t}{\int }_{t_i}^{t_i+\mathrm{\Delta }t}\tau (t)\mathrm{dt}$$

(21)

Approximate $\mathit{A}(t),\mathbf{\nu }(t),\mathit{v}(t-{\tau }_i)$ respectively:

$$\{\begin{array}{c}\mathit{A}(t)={\mathit{A}}_i+\frac{{\mathit{A}}_{i+1}-{\mathit{A}}_i}{\mathrm{\Delta }t}(t-t_i)\hfill \\ \mathbf{\nu }(t)={\mathbf{\nu }}_i+\frac{{\mathbf{\nu }}_{i+1}-{\mathbf{\nu }}_i}{\mathrm{\Delta }t}(t-t_i)\hfill \\ \mathbf{\nu }(t-{\tau }_i)={\beta }_i{\mathbf{\nu }}_{i-M}+{\alpha }_i{\mathbf{\nu }}_{i-M+1}\hfill \end{array}$$

(22)

where the approximate weight ${\beta }_i,{\alpha }_i$ and delay length $M$ are expressed as follows:

$${\beta }_i=m_i+\frac{1}{2}-\frac{\lambda \rho (1-\gamma \eta )}{N}$$

(23)

$${\alpha }_i=-m_i+\frac{1}{2}+\frac{\lambda \rho (1-\gamma \eta )}{N}$$

(24)

$$M=max\left\{m_i\right\}=max\left\{\mathrm{int}\left(\frac{{\tau }_i+\mathrm{\Delta }t/2}{\mathrm{\Delta }t}\right)\right\}$$

(25)

and

$$\eta =\frac{\lambda }{2\pi }{\int }_{\raisebox{1ex}{$2\pi i$}\!\left/ \!\raisebox{-1ex}{$\lambda $}\right.}^{\raisebox{1ex}{$2\pi \left(i+1\right)$}\!\left/ \!\raisebox{-1ex}{$\lambda $}\right.}\sin \left(t\right)dti=01\dots ,\lambda -1$$

(26)

Take equation (22) into equation (20):

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\nu }}_{i+1}=(\mathit{F}{\mathit{i}}_0+{\mathit{H}}_i){\mathbf{\nu }}_i+{\mathit{P}}_i{\mathbf{\nu }}_{i+1}}\\ +{\beta }_i{\mathit{R}}_i{\mathbf{\nu }}_{i-M}+{\alpha }_i{\mathit{R}}_i{\mathbf{\nu }}_{i-M+1}\hfill \end{array}$$

(27)

where

$$\{\begin{array}{c}{\mathit{H}}_i=(\mathit{F}{\mathit{i}}_1-\frac{2}{\mathrm{\Delta }t}\mathit{F}{\mathit{i}}_2+\frac{1}{\mathrm{\Delta }{t}^2}\mathit{F}{\mathit{i}}_3){\mathit{A}}_i+(\frac{1}{\mathrm{\Delta }t}\mathit{F}{\mathit{i}}_2-\frac{1}{\mathrm{\Delta }{t}^2}\mathit{F}{\mathit{i}}_3){\mathit{A}}_{i+1}\hfill \\ {\mathit{P}}_i=(\frac{1}{\mathrm{\Delta }t}\mathit{F}{\mathit{i}}_2-\frac{1}{\mathrm{\Delta }{t}^2}\mathit{F}{\mathit{i}}_3){\mathit{A}}_i+\left(\frac{1}{\mathrm{\Delta }{t}^2}\mathit{F}{\mathit{i}}_3\right){\mathit{A}}_{i+1}\hfill \\ {\mathit{R}}_i=(\mathit{F}{\mathit{i}}_1-\frac{1}{\mathrm{\Delta }t}\mathit{F}{\mathit{i}}_2){\mathit{A}}_i+\left(\frac{1}{\mathrm{\Delta }t}\mathit{F}{\mathit{i}}_2\right){\mathit{A}}_{i+1}\hfill \end{array}$$

(28)

and

$$\{\begin{array}{c}\mathit{F}{\mathit{i}}_0=e^{{\mathit{A}}_{\mathit{0}}\mathrm{\Delta }t}\hfill \\ \mathit{F}{\mathit{i}}_1={\mathit{A}}_0^{-1}(\mathit{F}{\mathit{i}}_0-\mathit{I})\hfill \\ \mathit{F}{\mathit{i}}_2={\mathit{A}}_0^{-1}(\mathit{F}{\mathit{i}}_1-\mathrm{\Delta }t\mathit{I})\hfill \\ \mathit{F}{\mathit{i}}_3={\mathit{A}}_0^{-1}(2\mathit{F}{\mathit{i}}_2-\mathrm{\Delta }{t}^2\mathit{I})\hfill \end{array}$$

(29)

where $I$ is the unitary matrix, discrete mapping can be defined as:

$$\begin{array}{c}{\mathrm{Z}}_{i+1}={\mathit{D}}_i{\mathit{Z}}_i\\ {\mathit{Z}}_{i+1}=\mathrm{col}{\mathbf{\nu }}_i,{\mathbf{\nu }}_{i-1}···{\mathbf{\nu }}_{i-M}\end{array}$$

(30)

and

$${\mathit{D}}_i=\left(\begin{array}{cccccc}{(\mathit{I}-\mathit{Pi})}^{-1}(\mathit{F}{\mathit{i}}_0+{\mathit{H}}_i)\hfill & 0\hfill & \cdots \hfill & 0\hfill & {(\mathit{I}-\mathit{Pi})}^{-1}{\alpha }_i{\mathit{R}}_i\hfill & {(\mathit{I}-\mathit{Pi})}^{-1}{\beta }_i{\mathit{R}}_i\hfill \\ \multicolumn{1}{c}{\mathit{I}}& 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \mathit{I}\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \cdots \hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \dots \hfill & 0\hfill & \mathit{I}\hfill & 0\hfill \end{array}\right)$$

(31)

Approximate Floquet transfer matrix $D=D_{\kappa -1}{D}_{\kappa -2}\dots D_0$ .

Based on Floquet theory, if the eigenvalue of transfer matrix is less than 1, the system is in a stable state, otherwise it is in an unstable state.^41,42

Model verification

Data in reference ⁴³ are used to verify the validity of the model presented in this paper. Dynamic milling system is two degrees of freedom, milling cutter with three cutter teeth, the diameter of 16 mm and helix angle of 30. Down-milling and constant speed machining are adopted. The cutting force coefficients and the first modal parameters are given in Table 1. In order to ensure the reliability of comparison and verification, the three degree of freedom model in this paper is reduced the dimension to fit two degrees of freedom. In addition, the modulation parameters of variable speed are set to 0.01 and 0.05 to make it infinitely close to the constant speed processing conditions. Then, the stability prediction is carried out by using the adjusted model and compared with the stability lobe diagram in reference. ⁴³ It can be seen from Figure 5 that when the modulation parameters are set very small, the stability lobe predicted by the model in this paper is similar to that predicted by the reference. Therefore, the model established in this paper can be verified and is effective for the analysis and prediction of stability.

Table 1.

Cutting parameters.

Type	Parameter items
Modal parameters	Directions	Natural frequency (Hz)	Damping ratios	Modal mass (kg)
	X	898.22	0.040041	1.576
	Y	852.51	0.036768	0.852
Cutting force coefficients	Kt (Mpa)		Kr (Mpa)
	1209.355		501.095

Figure 5.

Stability lobes of the dynamic system.

Result, analysis and discussion

The system parameters used in this paper are shown in Table 2.

Table 2.

Simulation parameters.

N	dc(mm)	helix [°]	Kt(N/m 2 )	Kn(N/m 2 )	Ka(N/m 2 )	εx	εy	εz	Cutting-direction
2	12.7	10/30/45/60	6 × 108	2 × 108	1.8 × 108	0.011	0.011	0.02	Up milling

Influence of different helix angle on milling stability

Firstly, the influence of helix angle on milling is studied. In this paper, two groups of variable speed modulation coefficients are selected to study the effect of different helix angles on chatter stability under two conditions. The solid line, dashed line, dash dot line and long dash line in Figure 6 represent the stability critical curve with different milling cutter helix angles respectively, and the green filling part is the case with the largest stable area. Many studies have confirmed the necessity of considering the helix angle, so this paper will not talk about it as a research point. As observed in Figure 6, it is easy to find that no matter the modulation coefficient is $\gamma =0.1$ , $\rho =0.2$ , or $\gamma =0.1$ , $\rho =0.5$ , choosing different helix angle has little effect on the overall stability lobes, and the more obvious increase of the stability region mostly occurs in the low speed range. Especially, there are many local maxima and local minima in the stability lobes. The change of local maxima has little influence on the actual machining, because it is impossible to select the machining parameters which are near the critical value in actual machining. However, the increase of local minimum value can effectively reduce the occurrence of chatter. As observed in Figure 6(b), when the milling cutter helix angle is 10 degrees, the local minimum value of critical curve is significantly higher than other helix angles.

Figure 6.

Three degrees of freedom milling stability lobes with different helix angles: (a) $\gamma =0.1$ and $\rho =0.2$ and (b) $\gamma =0.1$ and $\rho =0.5$ .

Influence of different radial cut in rate on milling stability

Figures 7 and 8 correspond to the stability lobes with immersion rate of 0.1 and 0.25, respectively. The blue dashed line represents the stability lobes of constant speed milling, the blue solid line represents the stability lobes of variable speed milling, and the green filled area represents the stable area of variable speed milling. On the basis of three degrees of freedom milling modeling, the comparison between variable and constant spindle speed milling stability lobes shows that the blue dashed line is basically in the green filling area. And in the low speed range, the stability critical curve of variable speed milling is obviously higher than that of constant speed, but the difference between them is not significant when the spindle speed is higher. Therefore, it is confirmed that variable speed milling has a positive effect on chatter suppression even when the modeling dimension is increased and the helix angle is considered.

Figure 7.

Three degrees of freedom milling stability lobes with radial immersion rate of 0.1: (a) constant speed and variable speed stability lobes when $\gamma =0.05$ , $\rho =0.1$ , (b) constant speed and variable speed stability lobes when $\gamma =0.05$ , (c) stability lobes of constant speed and variable speed when $\gamma =0.1$ , $\rho =0.1$ , and (d) stability lobes of constant speed and variable speed when $\gamma =0.1$ , $\rho =0.5$ .

Figure 8.

Three degrees of freedom milling stability lobes with radial immersion rate of 0.25: (a) constant speed and variable speed stability lobes when $\gamma =0.05$ , $\rho =0.1$ , (b) constant speed and variable speed stability lobes when $\gamma =0.05$ , $\rho =0.5$ , (c) stability lobes of constant speed and variable speed when $\gamma =0.1$ , $\rho =0.1$ , and (d) stability lobes of constant speed and variable speed when $\gamma =0.1$ , $\rho =0.5$ .

Compared with Figures 7 and 8, regardless of constant speed or variable speed milling, excessive radial immersion rate will directly increase the possibility of chatter, so the most basic method to suppress chatter can be to reduce immersion rate to avoid chatter. In particular, the stability critical curve is also affected by the modulation coefficients when the variable spindle speed is adopted. According to Figures 7 and 8(a) and (b), it is found that the increase of stable area is not obvious when variable speed machining is used, because $\gamma =0.05$ has a smaller value, The small amplitude of the change in the spindle speed results in little difference in the stability of variable speed and constant speed milling. According to Figures 7 and 8(c) and (d), it can be observed that with the increase of $\gamma$ , the stability region of variable speed milling is significantly larger than that of constant speed milling. With the increase of $\rho$ , The stability region of variable speed milling also increases to a certain extent. But in the high-speed range, a lower critical cutting depth is obtained compared to constant speed.

Influence of different modulation coefficient on milling stability

As shown in Figure 9, when speed modulation coefficient $\gamma$ is unchanged, the speed modulation coefficient $\rho$ is taken as different values respectively. The study found that when the value of $\gamma$ is small, as the speed modulation coefficient $\rho$ increases, the stable area will show a trend of increasing first and then decreasing. And when $\rho$ is about 0.7, the stable area is the most ideal. As shown in Figure 9(c), when $\gamma$ is 0.2, the value of $\rho$ needs to be a smaller value, otherwise stable area will obviously shrink, in addition, when $\rho$ is smaller, the machining stability of the high-speed range will be significantly improved. As shown in Figure 9(d), when the value of $\gamma$ is large, the stable area will be reduced sharply. Even if modulation $\rho$ is larger, it will not be helpful. This result is not correct, and it is in line with the actual situation: the spindle speed cannot change in a large range and it must be controlled within a reasonable range.

Figure 9.

Three degrees of freedom milling stability lobes with different modulation coefficients: (a) the stability lobes corresponding to different modulation coefficients $\rho$ when $\gamma =0.05$ (b) the stability lobes corresponding to different modulation coefficients $\rho$ when $\gamma =0.1$ , (c) the stability lobes corresponding to different modulation coefficients $\rho$ when $\gamma =0.05$ , and (d) the stability lobes corresponding to unreasonable modulation coefficients $\gamma$ and $\rho$ .

Conclusions

In this paper, the three degrees of freedom milling model is established by considering the cutting state corresponding to different helix angles. The chatter stability is analyzed by the improved semi-discretization method. The suppression effects of various machining conditions and modulation coefficients on chatter are compared. Concluded as follow:

1. The change of the milling cutter helix angle affects the position of the critical chatter curve to a certain extent. The inhibition of chatter is obvious under special machining conditions.

2. The modulation coefficient and its influence on the chatter stability are great. With the increase of the variation range of spindle speed, the influence of the variation frequency of speed on stability will be deepened, and too small variation of speed will weaken the gain effect of variable speed on stability. Therefore, the appropriate value of the modulation coefficient has important influence in the chatter suppression.

3. On the basis of three degrees of freedom milling modeling, the influence of variable speed on chatter stability is studied. It is found that the change amplitude of speed should not be too large, the change frequency of speed should not be too low or too high, and the modulation parameters that can make the stable area as large as possible should be selected within a reasonable range to avoid the occurrence of chatter.