Dynamic modeling and stability prediction in milling process of thin-walled workpiece with multiple structural modes

Abstract

Regenerative chatter can easily occur in the milling process of thin-walled workpiece due to the inherently low stiffness. This article aims to predict the stability of thin-walled workpiece in the milling process with a complete dynamic model. First, multiple structural modes of thin-walled workpiece are taken into consideration, and a complete dynamic model of thin-walled workpiece milling system is developed. Then, a numerical integration method is used to achieve the stability lobe diagrams of the milling system and identify the chatter frequency. Besides, the major structural mode, which is responsible for the occurrence of thin-walled workpiece chatter in the milling process, is predicted. A series of milling tests concerning a general cantilever plate are conducted, and the test results agree well with the predicted results, which shows the effectiveness of the proposed method. Finally, the effects of milling tool and structural modes on milling stability are discussed separately, which could provide theoretical basis for the dynamic modeling of thin-walled workpiece in milling process.

Introduction

Thin-walled workpieces are extensively used in aeronautic industry. Unfortunately, chatter occurs easily in the milling process of thin-walled workpiece due to its inherently low stiffness, which is the major limitation to product quality and productivity. Besides, milling chatter could cause some other negative effects, such as excessive noise, disproportionate tool wear and so on. Therefore, the chatter problem in the milling process of thin-walled workpiece is a topic of enormous interest in manufacturing.

Traditionally, milling tool is regarded as flexible with respect to the workpiece. ¹ Most researchers try to construct the stability lobe diagram (SLD) and then select optional machining parameters to ensure chatter-free operations. The regeneration phenomenon, which is the major cause of the milling chatter, is modeled mathematically by delay differential equation (DDE) with infinite dimensional state space, but its stability condition cannot be given analytically. Hence, several methods have been proposed to solve the problem. Smith and Tlusty ² used the peak-to-peak (PTP) force diagrams in time domain to predict the milling stability limit. Altintas and Budak ³ developed a zeroth order approximation (ZOA) method for chatter prediction in frequency domain. Since the time varying dynamic cutting force coefficient matrix can be approximated by its Fourier series component, the ZOA method allows a rapid calculation of SLD compared with time domain method. This method has been extended to variable-pitch cutters, ⁴ honed tools ⁵ and so on.

Other research scholars tried to predict the stability of the milling process in state space. For example, Insperger and Stepan^6,7 and Ding et al. ⁸ proposed the semi-discretization method (SDM) and the full-discretization method (FMD) to conduct stability analysis of the linear time-delay system, respectively. Zhang et al. ⁹ proposed the numerical integration method to obtain the SLD regarding the governing equations of the milling dynamic process as initial value problem for differential equation. To improve the convergence rate and computational accuracy, Niu et al. ¹⁰ presented a generalized form of the Runge–Kutta method (GRKM) based on a Volterra integral equation of the second kind. Li et al. ¹¹ built an extended dynamic milling model including mode coupling and process damping, and proposed a second-order SDM to simultaneously predict the stability lobes diagram and surface location error by solving this extended dynamic model. Yue et al. ¹² used FMD to predict the stability of the milling process by taking the characteristics of the contact between the tool and the workpiece into consideration. Qin et al. ²⁵ developed a Chebyshev-wavelet-based method for accurate and fast milling stability prediction.

In the milling process of thin-walled workpiece, the stiffness of thin-walled workpiece is much lower than the milling tool, which means that the milling tool–workpiece system changes from the traditional flexible–rigid system to the rigid–flexible system. Hence, the variations of thin-walled workpiece dynamic characteristics, mainly caused by the tool position, as well as multiple modes of the structure should be taken into consideration to achieve the SLD. Bravo et al. ¹³ presented a method to obtain the three-dimensional SLD when both milling tool and workpiece have similar dynamic behaviors, where the third dimension is the tool position along the workpiece. Seguy et al. ¹⁴ analyzed the relationship between chatter instability and surface roughness in the milling operations of thin walls. Mane et al. ¹⁵ analyzed the dynamic interaction of milling tool and thin-walled workpiece by a finite element approach. Wu et al. ¹⁶ investigated the stiffness-strengthening effect of longitudinal torsional ultrasonic milling of thin-plate structures by analyzing the forced vibrations of the thin-plate structure under conventional milling and ultrasonic milling loading. Sun and Jiang ¹⁷ established a novel dynamic model for thin-wall milling systems by considering force-induced deformation effect.

Most studies in the field have focused on milling chatter prediction, but only a few research efforts have tried to investigate the dynamic behavior of thin-walled workpiece, such as chatter frequency prediction, influence of the structural mode interaction on milling chatter and identification of the responsible structural mode for the occurrence of milling chatter. In this article, a clear picture about the chatter frequency arising in the milling process of thin-walled workpiece will be given by establishing a complete dynamic model of thin-walled workpiece. Besides, the effect of the structural mode interaction on the milling chatter will be studied to identify the responsible structural mode for the occurrence of milling chatter, which could provide a theoretical basis for the active control of thin-walled workpiece milling chatter.

The structure of this article is organized as follows. In section “Dynamic milling model,” a general cantilever plate model is employed, and the dynamic model of the milling process of thin-walled workpiece is developed. A method for the prediction of the SLD and the identification of the chatter frequency is presented in section “Stability analysis and chatter frequency.” In section “Experimental validation and results,” the experimental procedures are presented and the dynamic behavior of thin-walled workpiece is analyzed. The effects of milling tool and structural modes on milling stability are discussed separately in section “Discussions.” The conclusions of the article are drawn in the last section.

Dynamic milling model

In this section, a general cantilever plate model is used to investigate the dynamic behavior of thin-walled workpiece, as shown in Figure 1. Traditionally, the end mill can be considered as a stack of differential disk elements. The linear milling force model is employed, and the milling force acting on the jth tooth corresponding to a differential element can be expressed as ¹⁸

$$\left[\begin{array}{c}d{F}_t\left({\phi }_j\right)\hfill \\ \multicolumn{1}{c}{d{F}_r\left({\phi }_j\right)}\end{array}\right]=g\left({\phi }_j\right)\{\left[\begin{array}{c}{K}_{\mathrm{tc}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{rc}}}\end{array}\right]h\left({\phi }_j\right)+\left[\begin{array}{c}{K}_{\mathrm{te}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{re}}}\end{array}\right]\}\mathrm{dz}$$

(1)

where $d{F}_t({\phi }_j)$ and $d{F}_r({\phi }_j)$ are the differential forces in the tangential and radial directions; $K_{\mathrm{tc}}$ and $K_{\mathrm{rc}}$ are the cutting force coefficients, $K_{\mathrm{te}}$ and $K_{\mathrm{re}}$ are the edge force coefficients; $h({\phi }_j)$ is the instantaneous uncut chip thickness and the window function $g({\phi }_j)$ is used to explain state whether the jth tooth is cutting or not

$$g\left({\phi }_j\right)=\{\begin{array}{cc}1\hfill & {\phi }_{\mathrm{st}}<{\phi }_j\left(t\right)<{\phi }_{\mathrm{ex}}\hfill \\ \multicolumn{1}{c}{0}& \mathrm{others}\hfill \end{array}$$

(2)

${\phi }_{\mathrm{st}}$ and ${\phi }_{\mathrm{ex}}$ are the teeth enter and exit angle.

Figure 1.

Dynamic milling model for thin-walled workpiece.

The angular position of the jth tooth ${\phi }_j$ can be defined by

$${\phi }_j=\frac{2\pi \Omega }{60}t-\frac{\mathrm{ztan}\beta }{R}-(j-1)\frac{2\pi }{N}j=1,2,\dots ,N$$

(3)

where $\Omega$ is the spindle speed, $\beta$ denotes the cutter helix angle, R is the radius of the tool, z is the local axial direction depth and N is the number of the teeth.

Therefore, the regenerative effect is considered and the instantaneous uncut chip thickness $h({\phi }_j)$ can be written as

$$\begin{array}{c}h\left({\phi }_j\right)=f_t\sin \left({\phi }_j\right)+\left[\begin{array}{cc}\sin \left({\phi }_j\right)\hfill & \cos \left({\phi }_j\right)\hfill \end{array}\right]\hfill \\ \left[\begin{array}{c}[x_t(t-T)-x_t\left(t\right)]-[x_w(t-T)-x_w\left(t\right)]\hfill \\ \multicolumn{1}{c}{[y_t(t-T)-y_t\left(t\right)]-[y_w(t-T)-y_w\left(t\right)]}\end{array}\right]\hfill \end{array}$$

(4)

where $f_t$ denotes the feed rate. T is the time delay, namely the tooth passing period, $T=2\pi /N\Omega$ . $x_t$ and $y_t$ are the vibration displacements of the milling tool in the x and y directions and $x_w$ and $y_w$ are the vibration displacements of the thin-walled workpiece in the x and y directions.

According to Figure 1, the total milling force acting on the milling tool can be given as ¹⁹

$$\left[\begin{array}{c}{F}_x\left(t\right)\hfill \\ \multicolumn{1}{c}{F_y\left(t\right)}\end{array}\right]=\sum _{j=1}^N\underset{0}{\overset{a_p}{\int }}g\left({\phi }_j\right){\mathbf{T}}_j\{\left[\begin{array}{c}{K}_{\mathrm{tc}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{rc}}}\end{array}\right]h\left({\phi }_j\right)+\left[\begin{array}{c}{K}_{\mathrm{te}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{re}}}\end{array}\right]\}\mathrm{dz}$$

(5)

where $a_p$ is the axial depth of cut, ${\mathbf{T}}_j=\left[\begin{array}{cc}\cos \left({\phi }_j\right)\hfill & \sin \left({\phi }_j\right)\hfill \\ \multicolumn{1}{c}{-\sin \left({\phi }_j\right)}& \cos \left({\phi }_j\right)\hfill \end{array}\right].$

In this article, the first $n_1$ order structural modes of the milling tool and the first $n_2$ order structural modes of the thin-walled workpiece are employed, and the response of the milling tool and the thin-walled workpiece can be obtained based on the mode superposition method.

The modal shape of the milling tool is normalized at the tip of the milling tool, and the dynamic response of the milling tool can be determined by

$$\begin{array}{c}{m}_{\mathrm{pi}}{\stackrel{··}{p}}_i\left(t\right)+c_{\mathrm{pi}}{\stackrel{·}{p}}_i\left(t\right)+k_{\mathrm{pi}}{p}_i\left(t\right)={\varphi }_{\mathrm{pi}}{F}_x\left(t\right)i=1,2,\dots ,n_1\hfill \\ m_{\mathrm{qi}}{\stackrel{··}{q}}_i\left(t\right)+c_{\mathrm{qi}}{\stackrel{·}{q}}_i\left(t\right)+k_{\mathrm{qi}}{q}_i\left(t\right)={\varphi }_{\mathrm{qi}}{F}_y\left(t\right)i=1,2,\dots ,n_1\hfill \\ x_t\left(t\right)=\sum _{i=1}^{n_1}{\varphi }_{\mathrm{pi}}{p}_i\left(t\right)\hfill \\ y_t\left(t\right)=\sum _{i=1}^{n_1}{\varphi }_{\mathrm{qi}}{q}_i\left(t\right)\hfill \end{array}$$

(6)

where $p_i$ and $q_i$ are the ith order modal coordinates of the milling tool in the x and y directions. $m_{\mathrm{pi}}$ , $c_{\mathrm{pi}}$ , $k_{\mathrm{pi}}$ and ${\varphi }_{\mathrm{pi}}$ are the ith order modal mass, damping, stiffness and shape of the milling tool in the x direction. $m_{\mathrm{qi}}$ , $c_{\mathrm{qi}}$ , $k_{\mathrm{qi}}$ and ${\varphi }_{\mathrm{qi}}$ are the ith order modal mass, damping, stiffness and shape of the milling tool in the y direction.

In terms of the milling tool, with the excitation of the milling force at $x_0$ , the dynamic response of the thin-walled workpiece can be determined by

$$\begin{array}{c}{m}_{\mathrm{ri}}{\stackrel{··}{r}}_i\left(t\right)+c_{\mathrm{ri}}{\stackrel{·}{r}}_i\left(t\right)+k_{\mathrm{ri}}{r}_i\left(t\right)=-{\varphi }_{\mathrm{ri}}\left(x_0\right)\hfill \\ \multicolumn{1}{c}{F_x\left(t\right)i=1,2,\dots ,n_2}\\ m_{\mathrm{si}}{\stackrel{··}{s}}_i\left(t\right)+c_{\mathrm{si}}{\stackrel{·}{s}}_i\left(t\right)+k_{\mathrm{si}}{s}_i\left(t\right)=-{\varphi }_{\mathrm{si}}\left(x_0\right)\hfill \\ \multicolumn{1}{c}{F_y\left(t\right)i=1,2,\dots ,n_2}\\ x_w\left(t\right)=\sum _{i=1}^{n_2}{\varphi }_{\mathrm{ri}}{r}_i\left(t\right)\hfill \\ y_w\left(t\right)=\sum _{i=1}^{n_2}{\varphi }_{\mathrm{si}}{s}_i\left(t\right)\hfill \end{array}$$

(7)

where $r_i$ and $s_i$ are the ith order modal coordinates of the thin-walled workpiece at $x_0$ in the x and y directions. $m_{\mathrm{ri}}$ , $c_{\mathrm{ri}}$ , $k_{\mathrm{ri}}$ and ${\varphi }_{\mathrm{ri}}$ are the ith order modal mass, damping, stiffness and shape of the thin-walled workpiece in the x direction. $m_{\mathrm{si}}$ , $c_{\mathrm{si}}$ , $k_{\mathrm{si}}$ and ${\varphi }_{\mathrm{si}}$ are the ith order modal mass, damping, stiffness and shape of the thin-walled workpiece in the y direction.

The dynamic parameters of the milling tool and the thin-walled workpiece mentioned above can be determined by the modal test, as shown in section “Frequency response function acquisition.”

By substituting Equations (4) and (5) into Equations (6) and (7), the dynamic control equation of the milling system can be expressed as follows

$$\mathbf{M}\stackrel{··}{\mathbf{P}}\left(t\right)+\mathbf{C}\stackrel{·}{\mathbf{P}}\left(t\right)+\mathbf{K}\mathbf{P}\left(t\right)=\mathbf{\Phi }\left\{{\mathbf{K}}_c\left(t\right){\mathbf{\Phi }}^T\left[\mathbf{P}\left(t-T\right)-\mathbf{P}\left(t\right)\right]+{\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2\right\}$$

(8)

where

$$\begin{array}{c}\mathbf{M}=\mathrm{diag}(m_{p1},\dots ,m_{p{n}_1},m_{q1},\dots ,m_{q{n}_1},m_{r1},\dots ,m_{r{n}_2},m_{s1},\dots ,m_{s{n}_2}),\hfill \\ \mathbf{C}=\mathrm{diag}(c_{p1},\dots ,c_{p{n}_1},c_{q1},\dots ,c_{q{n}_1},c_{r1},\dots ,c_{r{n}_2},c_{s1},\dots ,c_{s{n}_2}),\hfill \\ \mathbf{K}=\mathrm{diag}(k_{p1},\dots ,k_{p{n}_1},k_{q1},\dots ,k_{q{n}_1},k_{r1},\dots ,k_{r{n}_2},k_{s1},\dots ,k_{s{n}_2}),\hfill \\ \mathbf{P}={[p_1,\dots ,p_{n1},q_1,\cdots ,q_{n1},r_1,\dots ,r_{n2},s_1,\cdots ,s_{n2}]}^T,\hfill \\ \multicolumn{1}{c}{\mathbf{\Phi }={\left[\begin{array}{c}{\mathbf{\phi }}_{\mathbf{p}}\hfill \\ {\mathbf{\phi }}_{\mathbf{q}}\hfill \\ \multicolumn{1}{c}{-{\mathbf{\phi }}_{\mathbf{r}}}\\ -{\mathbf{\phi }}_{\mathbf{s}}\hfill \end{array}\right]}^{\mathbf{T}},}\\ {\mathbf{\phi }}_{\mathbf{p}}={[{\varphi }_{p1},\dots ,{\varphi }_{p{n}_1}]}^T,{\varphi }_q={[{\varphi }_{q1},\dots ,{\varphi }_{q{n}_1}]}^T,\hfill \\ \multicolumn{1}{c}{{\mathbf{\phi }}_{\mathbf{r}}={[{\varphi }_{r1}\left(x_0\right),\dots ,{\varphi }_{r{n}_2}\left(x_0\right)]}^T,}\\ {\mathbf{\phi }}_{\mathbf{s}}={[{\varphi }_{s1}\left(x_0\right),\dots ,{\varphi }_{s{n}_2}\left(x_0\right)]}^T,\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_c\left(\mathbf{t}\right)=\sum _{j=1}^N\underset{0}{\overset{a_p}{\int }}g\left({\phi }_j\right){\mathbf{T}}_j\left[\begin{array}{c}{K}_{\mathrm{tc}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{rc}}}\end{array}\right]\cos \left({\phi }_j\right)\mathrm{dz},}\\ {\mathrm{\Delta }}_1=\sum _{j=1}^N\underset{0}{\overset{a_p}{\int }}g\left({\phi }_j\right){\mathbf{T}}_j\left[\begin{array}{c}{K}_{\mathrm{tc}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{rc}}}\end{array}\right]f_t\sin \left({\phi }_j\right)\mathrm{dz},\hfill \\ {\mathrm{\Delta }}_2=\sum _{j=1}^N\underset{0}{\overset{a_p}{\int }}g\left({\phi }_j\right){\mathbf{T}}_j\left[\begin{array}{c}{K}_{\mathrm{te}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{re}}}\end{array}\right]\mathrm{dz}.\hfill \end{array}$$

Since ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$ are the harmonic parts, which do not affect the stability of the milling process, Equation (8) can be rewritten as

$$\mathbf{M}\stackrel{··}{\mathbf{P}}\left(t\right)+\mathbf{C}\stackrel{·}{\mathbf{P}}\left(t\right)+\mathbf{K}\mathbf{P}\left(t\right)={\tilde{\mathbf{K}}}_c\left(t\right)\left[\mathbf{P}\left(t-T\right)-\mathbf{P}\left(t\right)\right]$$

(9)

where ${\tilde{\mathbf{K}}}_c\left(t\right)=\mathbf{\Phi }{\mathbf{K}}_c\left(t\right){\mathbf{\Phi }}^T$ .

Considering that the stiffness of the thin-walled workpiece in the feed direction is high enough, the vibration of the workpiece in the feed direction can be neglected. Besides, it is generally considered that the bending stiffness of the thin-walled workpiece is much smaller than that of the milling tool. In other words, the milling tool can be regarded as rigid, and the vibration of the milling tool can be neglected (the effect of the milling tool on the milling stability will be discussed in section “Discussions”). Hence, the instantaneous uncut chip thickness $h({\phi }_j)$ can be represented as

$$h\left({\phi }_j\right)=f_t\sin \left({\phi }_j\right)-[y_w(t-T)-y_w\left(t\right)]\cos \left({\phi }_j\right)$$

(10)

The dynamic control Equation (9) can be simplified as

$${\mathbf{M}}_w{\stackrel{··}{\mathbf{P}}}_w\left(t\right)+{\mathbf{C}}_w{\stackrel{·}{\mathbf{P}}}_w\left(t\right)+{\mathbf{K}}_w{\mathbf{P}}_w\left(t\right)={\tilde{\mathbf{K}}}_w\left(t\right)[{\mathbf{P}}_w(t-T)-{\mathbf{P}}_w\left(t\right)]$$

(11)

where ${\mathbf{M}}_w=\mathrm{diag}(m_{s1},m_{s2},\dots ,m_{s{n}_2})$ , ${\mathbf{C}}_w=\mathrm{diag}(c_{s1},c_{s2},\dots ,c_{s{n}_2})$ , ${\mathbf{K}}_w=\mathrm{diag}(k_{s1},k_{s2},\dots ,k_{s{n}_2})$ , ${\mathbf{P}}_w=[s_1,s_2,\cdots ,s_{n2}{]}^T$ , ${\tilde{\mathbf{K}}}_w(t)={\tilde{\mathrm{K}}}_c(t){\mathbf{\phi }}_{\mathbf{s}}{\mathbf{\phi }}_{\mathbf{s}}^{\mathrm{T}}$ , ${\mathbf{\phi }}_{\mathbf{s}}=[{\varphi }_{s1}(x_0),{\varphi }_{s2}(x_0),\dots ,{\varphi }_{s{n}_2}(x_0){]}^T$ , ${\tilde{\mathrm{K}}}_c(t)=\sum _{j=1}^N{\int }_0^{a_p}g({\phi }_j)[\begin{array}{cc}-\sin ({\phi }_j)\hfill & \cos ({\phi }_j)\hfill \end{array}]\left[\begin{array}{c}{K}_{\mathrm{tc}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{rc}}}\end{array}\right]\cos ({\phi }_j)\mathrm{dz}$ .

Stability analysis and chatter frequency

In this section, the numerical integration method ⁹ will be extended to the milling process of thin-walled workpiece by taking multiple structural modes into consideration.

With the transformation $\mathbf{Y}\left(t\right)=[\begin{array}{cc}{\mathbf{P}}_w\left(t\right)\hfill & {\stackrel{·}{\mathbf{P}}}_w\left(t\right)\hfill \end{array}{]}^T$ , Equation (11) can be represented as

$$\stackrel{·}{\mathbf{Y}}(t)=\mathbf{AY}(t)+\mathbf{B}(t)[\mathbf{Y}(t-T)-\mathbf{Y}\left(t\right)]$$

(12)

where

$\mathbf{A}=\left[\begin{array}{cc}0\hfill & \mathbf{I}\hfill \\ \multicolumn{1}{c}{-{{\mathbf{M}}_w}^{-1}{\mathbf{K}}_w}& -{{\mathbf{M}}_w}^{-1}{\mathbf{C}}_w\hfill \end{array}\right]$ , $\mathbf{B}\left(t\right)=\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{{{\mathbf{M}}_w}^{-1}{\tilde{\mathbf{K}}}_w\left(t\right)}& 0\hfill \end{array}\right]$ .

The time delay parameter T is equally divided into m small intervals, and the corresponding discretized time points are

$$t_k=t_1+(k-1)\tau k=1,2,\cdots ,m+1$$

(13)

where $\tau =T/m$ .

The solution of Equation (12) over each time interval $[t_k,t_{k+1}]$ can be obtained in the integration form

$$\mathbf{Y}\left(t_{k+1}\right)=\mathbf{Y}\left(t_k\right)+\underset{t_k}{\overset{t_{k+1}}{\int }}\mathbf{\Psi }\left(s\right)\mathrm{ds}$$

(14)

where $\mathbf{\Psi }(t)=\mathbf{\Theta }(t)\mathbf{Y}(t)+\mathbf{B}(t)\mathbf{Y}(t-T)$ , $\mathbf{\Theta }(t)=\mathbf{A}-\mathbf{B}(t)$ .

At time point $t_2$ , the fourth Lagrange interpolation polynomial is employed to approximate $\mathbf{\Psi }(t)$ , and then the value $\mathbf{Y}(t_2)$ can be obtained through numerical integration method. At the other points $t_k(k=3,\cdots ,m+1)$ , $\mathbf{Y}(t_k)$ can be calculated by the Simpson numerical integral method. Moreover, by considering the following relationship $\mathbf{Y}(t_1)=\mathbf{Y}(t_{m+1}-T)$ , a discrete map can be defined as

$$\mathbf{\Gamma }\left[\begin{array}{c}\mathbf{Y}\left(t_1\right)\hfill \\ \multicolumn{1}{c}{\mathbf{Y}\left(t_2\right)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\mathbf{Y}\left(t_{m+1}\right)}\end{array}\right]=\mathbf{\Lambda }\left[\begin{array}{c}\mathbf{Y}(t_1-T)\hfill \\ \multicolumn{1}{c}{\mathbf{Y}(t_2-T)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\mathbf{Y}(t_{m+1}-T)}\end{array}\right]$$

(15)

where

$$\mathbf{\Gamma }=\left[\begin{array}{ccccccc}\mathbf{I}\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{-(\mathbf{I}+\frac{251\tau }{720}{\mathbf{\Theta }}_1)}& (\mathbf{I}-\frac{323\tau }{360}{\mathbf{\Theta }}_2)\hfill & \frac{22\tau }{60}{\mathbf{\Theta }}_3\hfill & -\frac{53\tau }{360}{\mathbf{\Theta }}_4\hfill & \frac{19\tau }{720}{\mathbf{\Theta }}_5\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{-(\mathbf{I}+\frac{\tau }{3}{\mathbf{\Theta }}_1)}& -\frac{4\tau }{3}{\mathbf{\Theta }}_2\hfill & (\mathbf{I}-\frac{\tau }{3}{\mathbf{\Theta }}_3)\hfill & \hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& .\hfill & .\hfill & .\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & .\hfill & .\hfill & .\hfill & .\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & .\hfill & .\hfill & .\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & -(\mathbf{I}+\frac{\tau }{3}{\mathbf{\Theta }}_{m-1})\hfill & -\frac{4\tau }{3}{\mathbf{\Theta }}_m\hfill & (\mathbf{I}-\frac{\tau }{3}{\mathbf{\Theta }}_{m+1})\hfill \end{array}\right],$$

$$\mathbf{\Lambda }=\left[\begin{array}{ccccccc}\mathbf{I}\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{251\tau }{720}{\mathbf{B}}_1}& -\frac{323\tau }{360}{\mathbf{B}}_2\hfill & \frac{22\tau }{60}{\mathbf{B}}_3\hfill & -\frac{53\tau }{360}{\mathbf{B}}_4\hfill & \frac{19\tau }{720}{\mathbf{B}}_5\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{\frac{\tau }{3}{\mathbf{B}}_1}& \frac{4\tau }{3}{\mathbf{B}}_2\hfill & \frac{\tau }{3}{\mathbf{B}}_3\hfill & \hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& .\hfill & .\hfill & .\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & .\hfill & .\hfill & .\hfill & .\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & .\hfill & .\hfill & .\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \frac{\tau }{3}{\mathbf{B}}_{m-1}\hfill & \frac{4\tau }{3}{\mathbf{B}}_m\hfill & \frac{\tau }{3}{\mathbf{B}}_{m+1}\hfill \end{array}\right].$$

The transition matrix over one tooth passing period can be expressed as

$${\mathbf{F}}_{\mathrm{IM}}={\mathbf{\Gamma }}^{-1}\mathbf{\Lambda }$$

(16)

Finally, the stability of the milling system can be determined based on the Floquet theory. The relevant characteristic multiplier $\mu$ is the eigenvalue of ${\mathbf{F}}_{\mathrm{IM}}$ which is the largest in modulus. If $\left|\mu \right|<1$ , the milling system is stable, and the milling operation is chatter free; otherwise, it is unstable (Figure 2).

Figure 2.

Unstable behavior in complex plane: (a) Hopf bifurcation, (b) period-doubling bifurcation and (c) period one bifurcation.

The arising chatter frequency for an unstable milling process can be determined for the case $\left|\mu \right|=1$ . If the relevant characteristic multiplier $\mu$ is a complex pair, which corresponds to the secondary Hopf bifurcation, then the chatter frequencies are^20,21

$$f_{\mathrm{hopf}}=\{\pm \frac{\omega }{2\pi }+l\frac{\Omega }{60}\}\mathrm{Hz}l=0,\pm 1,\pm 2\cdots$$

(17)

where $\omega =\mathrm{Im}(\ln (\mu ))/T$ .

If the relevant characteristic multiplier is $\mu =1$ and $\mathrm{\omega }=0$ , this case corresponds to the period one bifurcation, which is topologically equivalent to the saddle-node bifurcation of autonomous system. The chatter frequencies are

$$f_{p1}=\{0+l\frac{2\pi }{\tau }\}\mathrm{rad}/\mathrm{s}=\{0+l\frac{\Omega }{60}\}\mathrm{Hz}l=0,\pm 1,\pm 2\cdots$$

(18)

If the relevant characteristic multiplier is $\mu =-1$ and $\mathrm{\omega }=\pi /T$ , this case corresponds to the period-doubling bifurcation or flip bifurcation. The chatter frequencies are

$$f_{p2}=\{\frac{\pi }{T}+l\frac{2\pi }{\tau }\}\mathrm{rad}/\mathrm{s}=\{\frac{\Omega }{120}+l\frac{\Omega }{60}\}\mathrm{Hz}l=0,\pm 1,\pm 2\cdots$$

(19)

Experimental validation and results

In this section, a cantilever plate is employed and the impact test is conducted to determine its modal parameters. Then the SLD of the thin-walled workpiece and the chatter frequency are predicted based on the method introduced in section “Stability analysis and chatter frequency.” Finally, a series of milling tests are carried out to verify the predicted results.

A cantilever plate, which is made of 7050 Aluminum, is used in this article, and the dimension of the plate is $120\mathrm{mm}\times 60\mathrm{mm}\times 4\mathrm{mm}$ . Besides, a three-fluted flat end mill cutting tool, which is made of carbide, is used to conduct the milling test. The dimension parameters of the cutting tool are shown in Table 1.

Table 1.

Dimension parameters of cutting tool.

Helix angle (°)	Diameter (mm)	Length of tool (mm)	Length of tool blade (mm)	Tool overhang (mm)
$35$	$16$	$100$	$45$	$55$

Milling force coefficients acquisition

To predict the milling SLD and the chatter frequency, milling force coefficients should be calibrated in advance.

During the milling test of calibrating the milling force coefficients, milling parameters are selected to achieve stable milling process. The instantaneous milling forces are measured by the dynamometer Kistler 9255C and recorded by the software HRsoft. Based on the method proposed by Budak et al., ²² the milling force coefficients are calibrated through the test results, as listed in Table 2.

Table 2.

Milling force coefficients.

$K_{\mathrm{tc}}(\mathrm{N}/\mathrm{m}{\mathrm{m}}^2)$	$K_{\mathrm{rc}}(\mathrm{N}/\mathrm{m}{\mathrm{m}}^2)$	$K_{\mathrm{te}}(\mathrm{N}/\mathrm{mm})$	$K_{\mathrm{re}}(\mathrm{N}/\mathrm{mm})$
762.5	254.3	17.7	15.8

Frequency response function acquisition

As mentioned before, the diameter of the tool is $16\mathrm{mm}$ and the tool overhang is $55\mathrm{mm}$ . Hence the stiffness of the tool is expected to be much higher than that of the flexible workpiece. In other words, the workpiece is more flexible than the tool and the vibration of tool can be neglected. Consequently, only the modal parameters of the thin-walled workpiece are used to predict the SLD and the chatter frequency.

The variations of the thin-walled workpiece dynamic characteristics caused by the tool position are taken into consideration. Since the thin-walled workpiece is an axisymmetric structure, five measurement points (P1–P5) on the first half of the thin-walled workpiece in feed direction are selected, as shown in Figure 3(a). The direct FRFs of the measurement points are obtained by the impact tests.

Figure 3.

(a) Measurement points on the workpiece and (b) experimental setup for impact test.

During the impact tests, the impulse hammer Kistler 9724A2000 is used as the excitation source, and the response of the plate is acquired by the eddy current sensor. The response signals are measured and processed by the LMS SCADAS XS data acquisition system and the software LMS Test. Lab 14A, as shown in Figure 3(b). The impact test of each measurement point is conducted 10 times to eliminate the disturbance of noise. The sampling bandwidth, namely the cutoff frequency, is set as $3000\mathrm{Hz}$ .

The impact test results are illustrated in Figure 4. It is noted that point P5 is situated on the node of the second modal shape of the cantilever plate. Hence, the second-order mode cannot be seen in the FRF of point P5. Similarly, point P3 is located at the node of the third modal shape of the cantilever plate. Hence, the third-order mode cannot be detected in the FRF of point P3. The modal parameters are identified by normalizing the mode shape of the driving point based on the peak-picking method,^22–24 and the results are listed in Table 3.

Figure 4.

Measured FRFs of the thin-walled workpiece.

Table 3.

Modal parameters of the thin-walled workpiece.

Point	Mode	Frequency (Hz)	Damping ration (%)	Stiffness $(\mathrm{N}/\mathrm{m})$
P1	First	850.3	1.02	7.56e5
	Second	1363.1	0.15	6.5e5
	Third	2624	0.14	2.87e6
P2	First	850.3	1.1	7.38e5
	Second	1363.1	0.19	9.44e5
	Third	2624	0.1	1.23e7
P3	First	850.3	1.33	7.2e5
	Second	1363.1	0.17	2.06e6
P4	First	850.3	1.14	6.9e5
	Second	1363.1	0.19	5.21e6
	Third	2624	0.15	7.34e6
P5	First	850.3	1.18	6.41e5
	Third	2624	0.11	3.94e6

Milling stability identification and chatter frequency prediction

After obtaining the milling force coefficients and the dynamic characteristics of the thin-walled workpiece, the numerical integration method, which has been introduced in section “Stability analysis and chatter frequency,” is employed to calculate the SLD and the chatter frequency of each measurement point.

The stability charts are determined over a 240×200 sized grid of the spindle speed and depth of cut. The discretized number is set as 40, and the radial depth of cut is 1 mm, down milling. The calculated results are shown in Figures 5–9. ▲ means the milling process is stable when the tool passes the point, while • means chatter occurs.

Figure 5.

SLD and chatter frequency of point P1.

Figure 6.

SLD and chatter frequency of point P2.

Figure 7.

SLD and chatter frequency of point P3.

Figure 8.

SLD and chatter frequency of point P4.

Figure 9.

SLD and chatter frequency of point P5.

Clearly, there exist remarkable variations among the SLD of the five measurement points. With the increase of the structural stiffness from point P1 to point P5, the critical depth of cut is also enlarging. Besides, there does not exist any flip bifurcation.

Experimental verification

To verify the predicted SLD and the chatter frequency, a series of milling tests are conducted. Four typical parameter points (point A1, point A2, point A3 and point A4) on the predicted SLD are chosen, and the test results are discussed as below. In addition, the structural mode, which is responsible for the occurrence of the milling chatter, is also analyzed. During the milling tests, the vibration displacement signal of point P6 in Figure 3 is used to analyze the dynamic behavior of thin-walled workpiece.

The measured vibration displacement signal with the milling parameters of point A1 $(6300\mathrm{rpm},0.8\mathrm{mm})$ can be found in Figure 10. Theoretically, the milling process under the milling parameters of point A1 is unstable, which can be seen from Figures 5 to 9. The first-order structural mode of the thin-walled workpiece should be responsible for the occurrence of the milling chatter, and the predicted main chatter frequency is $860.9\mathrm{Hz}.$

Figure 10.

Test results with the milling parameters of point A1.

It is obvious that the whole milling process is unstable, which coincides with the predicted result. Three sub-signals $s1$ , $s2$ and $s3$ are selected to further verify the state of the milling process, and their frequency spectrums are illustrated in Figure 10. The measured main chatter frequency is $896\mathrm{Hz}$ , which is close to the first natural frequency of the thin-walled workpiece. Hence, it can be said with certainty that the instability of the milling process is caused by the first-order structural mode. Obviously, the test result is basically identical to the predicted result. The difference between the predicted chatter frequency and the test result can be attributed to cutter runout and so on. There are also some other chatter frequencies, such as $1211$ and $1526\mathrm{Hz}$ and these chatter frequencies are modulated by the harmonics of the rotation frequency $f=105\mathrm{Hz}$ .

The quality of the machined surface is shown in Figure 11. Evidently the surface finish is poor, and chatter marks are clearly visible.

Figure 11.

Surface finish of the workpiece with the milling parameters of point A1.

In terms of the milling process with the parameters of point A2 $(7200\mathrm{rpm},0.5\mathrm{mm})$ , it is predicted that chatter occurs when the tool located at points P1–P4 and the chatter frequency is $1365\mathrm{Hz}$ , which means that the chatter is caused by the second-order structural mode. But the milling process is stable at point P5.

The test result, as shown in Figure 12, confirms that the milling process is unstable from P1 to P4, and then gradually turns to stable from P4 to P5 as the tool advances along the tool path. But, milling chatter appears again at about $5.08\mathrm{s}$ and lasts to the end of the milling process. The frequency spectrums of the sub-signal $\mathrm{s}4$ and $\mathrm{s}5$ in Figure 12 show that the main chatter frequency is $1392\mathrm{Hz}$ . That is to say, the main chatter frequency is excited by the second-order structural mode of the thin-walled workpiece, which is in accord with the predicted result.

Figure 12.

Test results with the milling parameters of point A2.

This result can be also seen from the machined surface in Figure 13. Clear chatter marks can be observed around the two ends of the thin-walled workpiece, whereas the middle section of the workpiece shows a smooth surface.

Figure 13.

Surface finish of the workpiece with the milling parameters of point A2.

As for the point A3 $(4200\mathrm{rpm},0.7\mathrm{mm})$ , the predicted results demonstrate that the instability of the milling process from points P1 to P3 is caused by the second-order structural mode, and the chatter frequency is $1377\mathrm{Hz}$ . But the third-order structural mode is responsible for the instability of points P4 and P5, and the chatter frequency is $2638\mathrm{Hz}$ .

The test result can be found in Figure 14, which shows that the whole milling process is unstable. In simple terms, the milling process can be divided into three stages. Specifically speaking, at the beginning of the milling process, the main chatter frequency is $1418\mathrm{Hz}$ , which means that the milling chatter is excited by the second-order structural mode of the thin-walled workpiece.

Figure 14.

Test results with the milling parameters of point A3.

Then the chatter frequencies excited by the second-order structural mode fade away after the tool passes by point P3, and the chatter frequency (2675 Hz) excited by the third-order structural mode appears and dominates the milling process. In addition, the amplitude of the displacement signal decreases due to the reason that the stiffness of the workpiece in the middle section is much higher than those of the two ends. Obviously, the test result is basically identical to the predicted result. The surface finish is shown in Figure 15, and strong vibration marks can be observed, especially around the two ends.

Figure 15.

Surface finish of the workpiece with the milling parameters of point A3.

In terms of point A4 $(6900\mathrm{rpm},0.9\mathrm{mm})$ , the test result can be seen from Figures 16 and 17. It is evident that the whole milling process is also unstable. Besides, the milling chatter at points P1–P3 is caused by the second-order structural mode, and the chatter frequency is $1368\mathrm{Hz}$ . The third-order structural mode is responsible for the instability of points P4 and P5, and the chatter frequency is $2679\mathrm{Hz}$ . Hence the unstable pattern of point A4 is similar to point A3, and the test result also agrees well with the predicted result.

Figure 16.

Test results with the milling parameters of point A4.

Figure 17.

Surface finish of the workpiece with the milling parameters of point A4.

Discussions

Effect of milling tool on milling stability

To simplify the milling process, the milling tool is regarded as rigid and its vibration is neglected in section “Dynamic milling model.” In this section, the assumption will be evaluated to provide theoretical basis for simplifying the milling process.

When the milling tool is considered in the modeling process, the instantaneous uncut chip thickness $h({\phi }_j)$ can be rewritten as

$$\begin{array}{c}h\left({\phi }_j\right)=f_t\sin \left({\phi }_j\right)+\left[\begin{array}{cc}\sin \left({\phi }_j\right)\hfill & \cos \left({\phi }_j\right)\hfill \end{array}\right]\hfill \\ \left[\begin{array}{c}[x_t(t-T)-x_t\left(t\right)]\hfill \\ \multicolumn{1}{c}{[y_t(t-T)-y_t\left(t\right)]-[y_w(t-T)-y_w\left(t\right)]}\end{array}\right]\hfill \end{array}$$

(20)

The transfer matrix $\mathbf{\Phi }$ can be rewritten as

$$\mathbf{\Phi }={\left[\begin{array}{cc}{\mathbf{\phi }}_{\mathbf{p}}\hfill & \hfill \\ {\mathbf{\phi }}_{\mathbf{q}}\hfill \\ -{\mathbf{\phi }}_{\mathbf{s}}\hfill \end{array}\right]}^{\mathbf{T}}$$

(21)

As mentioned in section “Milling force coefficients acquisition,” the diameter of the milling tool used in this article is $16\mathrm{mm}$ , and the tool overhang is $55\mathrm{mm}$ . The modal impact test is employed again to obtain the model parameters of the milling tool. The impulse hammer BK 8206 and acceleration sensor PCB 352C23 (with mass $0.2\mathrm{g}$ ) are used for modal testing, and the experiment diagram is shown in Figure 18.

Figure 18.

Modal test of milling tool.

The impact test results are illustrated in Figure 19, and then the modal parameters of the milling tool are calibrated, as shown in Table 4. Obviously, the stiffness of the milling tool is much higher than that of the thin-walled workpiece.

Figure 19.

Measured FRFs of the milling tool.

Table 4.

Modal parameters of milling tool.

Direction	Frequency (Hz)	Damping ration (%)	Stiffness (N/m)
x	2114	2.8	5.45e6
y	2109	2.28	5.77e6

To evaluate the effect of the milling tool on the stability of the milling system, the numerical integration method, which have been introduced in section “Stability analysis and chatter frequency,” is employed again to calculate the SLD for its high accuracy and efficiency.

Taking the measurement points P1 and P3 for example, the calculated results are shown in Figure 20(a) and (b). Obviously, there is no difference between the SLD of the measurement points P1 and P3 with and without the consideration of the milling tool. Therefore, the assumption which neglects the vibration of the milling tool is reasonable in this article. Moreover, according to the theory of vibration, when the stiffness of the milling tool is 10 times larger than that of the workpiece, the effect of the milling tool on milling stability can be neglected.

Figure 20.

Milling stability lobes: (a) P1 and (b) P3.

Effect of structural mode on milling stability

As mentioned above, multiple structural modes of the thin-walled workpiece should be taken into consideration to study the stability of the milling system. In section “Frequency response function acquisition,” the direct FRFs of the measurement points are obtained by the impact tests, and the cutoff frequency of the impact tests is set as 3000 Hz. The test results show that there are three modes in the range of the cutoff frequency. Hence, only the first three modes are used to calculate the SLD of the thin-walled workpiece. The fourth and the higher order modes are neglected. In this section, the effect of structural mode on milling stability is investigated to provide theoretical basis for the selection of the cutoff frequency.

The fourth-order and the higher order modal parameters of the thin-walled workpiece can be also obtained by the impact test with a lager cutoff frequency. For example, the fourth-order modal frequency of the measurement point P1 is $5150\mathrm{Hz}$ , and the fourth-order modal stiffness and damping ratio of the measurement point P1 are $3.265{\mathrm{e}}^7\mathrm{N}/\mathrm{m}$ and $0.001$ , respectively.

Then the SLD of the thin-walled workpiece in the measurement point P1 can be calculated with different structural modes, as shown in Figure 21. It is noticed that there is a significant difference among the SLD with the first mode, the first two modes and the first three modes. But the SLD with the first three modes and the first four modes are identical. Hence the fourth-order and the higher order structural modes have no influence on the SLD of the thin-walled workpiece in this article.

Figure 21.

Milling stability lobes with different modes.

Therefore, the cutoff frequency should be carefully chosen to reduce the truncation error of the mode superposition method. Generally, the cutoff frequency should be at least 6 times of the exciting frequency.

Conclusion

The dynamic behavior of the thin-walled workpiece with multiple structural modes in the milling process is studied in this article. A complete cantilever plate model is proposed, and the numerical integration method is employed and extended to thin-walled workpiece milling process. In the meantime, the variations of the thin-walled workpiece dynamic characteristics caused by the tool position are considered. The stability limits are achieved, and a clear picture about the chatter frequency arising in thin-walled workpiece milling process is given. Then, the responsible structural mode for the instability of the milling process is identified. Moreover, the mechanism of the interaction of multiple structural modes when the thin-walled milling chatter occurs is revealed. The predicted results are verified by a series of milling tests, which could provide a theoretical basis for the active control of the thin-walled workpiece milling chatter. Finally, the effects of milling tool and structural modes on milling stability are discussed separately, which could provide theoretical basis for the neglect of milling tool and the selection of the cutoff frequency.