Chatter stability and surface location error predictions in milling with mode coupling and process damping

Abstract

Together with machining chatter, surface location error induced by forced vibration may also inhibit productivity and affect workpiece surface quality in milling process. Addressing these issues needs the combined consideration of stability lobes diagram and surface location error predictions. However, mode coupling and process damping are seldom taken into consideration. In this article, an extended dynamic milling model including mode coupling and process damping is first built based on classical 2-degree-of-freedom dynamic model with regeneration. Then, a second-order semi-discretization method is proposed to simultaneously predict the stability lobes diagram and surface location error by solving this extended dynamic model. The rate of convergence of the proposed method is also investigated. Finally, a series of experiments are conducted to verify the veracity of the extended dynamic model. The modal parameters including direct and cross terms are identified by impact experiments. Via experimental verification, the experimental results show a good correlation with the predicted stability lobes diagram and surface location error based on the extended dynamic model. Also, the effects of mode coupling and process damping are revealed. Mode coupling increases the whole stability region; however, process damping plays a vital role in stability improvement mainly at low spindle speeds.

Introduction

Increasing demands on manufacturing precision products require high-efficiency and high-accuracy machining.^1,2 Tool path optimization and machining parameters’ optimization are two common technologies that meet the demands. For tool path optimization, Li and colleagues^3,4 constructively proposed a rank 2 tensor–based global optimal tool path generation method, which can effectively improve the machining efficiency of the workpiece with complex free-form surfaces. Sun et al. ⁵ proposed a smooth tool path generation method for five-axis machining of triangular mesh surface with nonzero genus. For machining parameters’ optimization, the researches mainly focus on the prediction of stability lobes diagram (SLD) and surface location error (SLE).^6,7 Because both chatter and forced vibrations have a great effect on machining process, they may degrade machining performance and deteriorate workpiece surface quality. Guo et al. ⁸ studied SLD in milling with multi-delays considering helix angle effect. Shi et al. ⁹ predicted SLD for face milling processes considering the dynamic interaction between cutting tool and workpiece. Schmitz and colleagues^10,11 revealed the effect of helix angle and run-out on SLE. Generally, the combined consideration of SLD and SLE is important for choosing optimal machining parameters.

Mathematically, classic dynamic milling model based on regeneration can be formulated as delay differential equations (DDEs) with time-varying coefficients. By solving this type of equation, the SLD and SLE can be obtained. Throughout the past few decades, tremendous efforts have been made for accurately and efficiently predicting SLD and SLE. For SLD prediction, a lot of remarkable methods have been proposed such as zero-order analytical (ZOA) method, ¹² temporal finite element analysis (TFEA) method, ¹³ multi-frequency solution (MFS), ¹⁴ semi-discretization method (SDM),^15,16 full-discretization method (FDM),^17,18 and Runge–Kutta methods. ¹⁹ For SLE prediction, the major algorithms focus on numerical simulation method, ²⁰ frequency domain method, ¹⁰ harmonic balance method, ²¹ TFEA, ²² and FDM. ²³ Especially, TFEA and FDM can achieve simultaneous predictions of SLD and SLE.

However, mode coupling and process damping were not taken into account in the above works. As is well known, regenerative chatter is the most familiar form of self-excited vibration, which appears as a time-delay item in the chip thickness equation from the modeling standpoint. Likewise, mode coupling is also an indispensable chatter mechanism which results from the coupled effect among the natural frequencies of the different vibration modes.^24,25 Tlusty and Polacek ²⁶ first discussed that chatter may occur owing to non-regenerative mechanisms such as mode coupling. Gradišek et al. ²⁷ revealed that even weak mode coupling influences the predicted stability boundary greatly via numerical simulation. Zhang et al. ²⁸ found that mode coupling and regenerative effect actually co-exist during the practical milling process, and the usually neglected mode coupling has a great effect on the stability boundary. Recently, Lu et al. ²⁹ studied the effects of tool orientation on the five-axis milling process stability with flat-end cutters, and the formulation for the dynamic milling system is developed by simultaneously considering mode coupling and regenerative effect. Apart from mode coupling, process damping originating from the interaction between tool and workpiece has also been one of the hot issues in machining dynamics.^30,31 Process damping is well known as a beneficial phenomenon that may suppress chatter and increase stability, especially in low-speed milling of difficult-to-machine materials such as titanium and nickel alloy. Given its complicated mechanism, most of the previous researches ignored this beneficial phenomenon when modeling the dynamic milling process. However, this small negligence may lead to significant errors in stability prediction at low spindle speeds. Wu ³² modeled process damping forces as a function of the tool flank–workpiece indentation volume and a material dependent coefficient. Subsequently, the principal research work focused on the coefficient calibration and indentation volume calculation. Chiou and Liang ³³ approximated the indentation volume as a function of cutting speed, vibration velocity, and wear land length in the case of small amplitude vibration. Ahmadi and Ismail^34,35 analytically derived the computation procedure of indentation volume and integrated equivalent linear viscous damper of process damping into the MFS and SDM to obtain the SLD in milling. Tunç and Budak ³⁶ proposed a practical identification method for process damping coefficient through energy balance formulation. The theoretical expressions of indentation volume is established by including some important geometrical parameters such as tool hone radius, clearance angle, separation angle, and geometry of clearance face. Ahmadi and Altintas ³⁷ presented an identification method based on frequency domain decomposition, where process damping coefficient was estimated using vibrations measured at two locations of the tool during chatter free orthogonal turning tests. They calculated wear land length as the summation of three segments by considering the effects of hone radius, separation angle, and clearance angle. Recently, Wan et al. ³⁸ proposed a generalized method for exploring the plowing mechanism of both static and dynamic cuts. The process damping coefficient can be directly identified using the static plowing forces, which are separated from the measured milling forces.

In this article, an extended dynamic model of the milling process with mode coupling and process damping is presented. It is noteworthy that SDM can be applied here due to its versatility in analysis of advanced milling models. Hence, a second-order SDM is proposed to simultaneously obtain the SLD and SLE by solving the extended dynamic milling model. The remainder of this article is organized as follows. Based on classical 2-degree-of-freedom (2-DOF) dynamic milling model with regeneration, an extended dynamic model including mode coupling and process damping is first built in section “Extended dynamic milling model with mode coupling and process damping.” In this extended dynamic model, regenerative effect varies instantaneous uncut chip thickness which, in turn, generates a variable cutting force. Process damping effect is described with an equivalent linear viscous damper model in Ahmadi and Ismail, ³⁴ and the cross terms of modal parameters are nonzero when mode coupling is considered. Through further transformation, the extended dynamic model in the state-space form can be obtained for the convenience of solution. In section “Milling stability and SLE analysis,” a second-order SDM is proposed by approximating time-delay term with second-order Newton interpolation polynomial. Using the proposed method, the extended dynamic model in the state-space form can be solved to simultaneously obtain SLD and SLE via Floquet theory and fixed point theory, respectively. The rate of convergence of the proposed method is investigated through the comparison with the first-order SDM ¹⁶ and the second-order FDM. ¹⁷ Section “Experimental verification” first calibrates the cutting force coefficients by experiments. The modal parameters including direct terms and cross terms are then identified by impact experiments. Later, a series of experiments are conducted to verify the predicted SLD and SLE based on the extended dynamic model. Also, the effects of mode coupling and process damping are revealed. Section “Conclusion” summarizes the article.

Extended dynamic milling model with mode coupling and process damping

Milling system, as shown in Figure 1, can be regarded to have 2 DOFs in the X and Y directions. The tool is assumed to be flexible with respect to the rigid workpiece. For helical end mills with N teeth and equal pitch angle, the instantaneous immersion angle is actually variable along the axis of the tool. However, considering that the axial cutting depth is small compared with the diameter of the tool, helix angle is ignored according to Altintas. ¹ The instantaneous immersion angle of the jth tooth is defined as a function of the spindle speed Ω and time t

$${\phi }_j=\frac{2\pi }{60}\Omega t+\frac{2\pi }{N}(j-1)$$

(1)

Figure 1.

The schematic diagram of 2-DOF dynamic milling system with process damping.

The cutting forces F_r and F_t acting on the jth tooth are attributed to two physical mechanisms: the chip shearing and process damping effect. The shearing forces, with subscript s, are generally supposed to be proportional to the chip area

$$F_{ts}=K_{t}g({\phi }_j)h\left(t\right)a_p,F_{\mathrm{rs}}=K_{r}g({\phi }_j)h\left(t\right)a_p$$

(2)

where a_p represents the axial cutting depth. K_t and K_r are tangential and radial shearing force coefficients, separately. The unit step function, $g({\phi }_j)$ , is equal to 1 when the jth tooth is engaged in the cut and 0 otherwise. Owing to the regeneration effect shown in Figure 1, the uncut instantaneous chip thickness h(t) is determined using the following equation

$$h\left(t\right)=f_t\sin {\phi }_j+r\left(t\right)-r(t-T),T=60N\Omega$$

(3)

where $f_t\sin {\phi }_j$ denotes the mean chip thickness, r(t − T) is the vibration in the radial direction by the previous tooth, and T is the tooth passing period. r(t) represents the vibration by the current tooth which can be expressed in the light of the dynamic displacements in the X and Y directions

$$r\left(t\right)=\left[\begin{array}{cc}\hfill \sin {\phi }_j\hfill & \hfill \cos {\phi }_j\hfill \end{array}\right]\mathbf{q}\left(t\right),\mathbf{q}\left(t\right)={\left[\begin{array}{cc}\hfill x\left(t\right)\hfill & \hfill y\left(t\right)\hfill \end{array}\right]}^T$$

(4)

where q(t) is the dynamic displacement vector. As for process damping force with subscript p, it is assumed to be proportional to the indented volume between the workpiece and the flank face of the tool. Following the development in Ahmadi and Ismail, ³⁴ the equivalent viscous model is utilized as follows

$$F_{tp}=\mu F_{rp},F_{rp}=g({\phi }_j)C_{eq}\stackrel{·}{r},C_{eq}=K_{sp}{C}_d\frac{a_p}{v}$$

(5)

where $\mu$ denotes the coulomb friction coefficient, $K_{\mathrm{sp}}$ is the specific indentation force coefficient, and $v$ is the tangential velocity. The damping shape factor, $C_d$ , is computed from the geometrical parameters of the tool and the surface undulations. Through a lot of experiments, Ahmadi and Ismail ³⁹ generated a database for C_d. In view of small amplitude vibration assumption, the damping shape factor can be expressed as $C_d=0.25W^2$ accordingly, where W is the tool wear land width, as shown in Figure 1.

The total cutting forces F in the X and Y directions can be then calculated by the summation of the cutting force F_r and F_t on all teeth

$$\begin{array}{c}{F}_t=F_{ts}+F_{tp},F_r=F_{rs}+F_{rp}\hfill \\ \mathit{F}=\left\{\begin{array}{c}\hfill F_x\hfill \\ \multicolumn{1}{c}{F_y}\end{array}\right\}=\sum _{j=1}^N\left[\begin{array}{cc}\hfill -\cos {\phi }_j\hfill & \hfill -\sin {\phi }_j\hfill \\ \multicolumn{1}{c}{\sin {\phi }_j}& \hfill -\cos {\phi }_j\hfill \end{array}\right]\left\{\begin{array}{c}\hfill F_t\hfill \\ \multicolumn{1}{c}{F_r}\end{array}\right\}\hfill \end{array}$$

(6)

Rearranging the equation of total cutting forces in matrix form yields the following

$$\begin{array}{c}\mathit{F}=a_p\mathbf{H}(\left[\begin{array}{c}\hfill f_t\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]+\mathbf{q}\left(t\right)-\mathbf{q}(t-T))+C_{eq}\mathbf{G}\stackrel{·}{\mathbf{q}}\left(t\right)\hfill \\ \mathbf{H}=\left[\begin{array}{cc}\hfill -h_{\mathrm{xx}}\hfill & \hfill -h_{xy}\hfill \\ \multicolumn{1}{c}{-h_{yx}}& \hfill -h_{yy}\hfill \end{array}\right],\mathbf{G}=\left[\begin{array}{cc}\hfill -g_{xx}\hfill & \hfill -g_{xy}\hfill \\ \multicolumn{1}{c}{-g_{yx}}& \hfill -g_{yy}\hfill \end{array}\right]\hfill \end{array}$$

(7)

where

$$\begin{array}{c}\left\{\begin{array}{c}\hfill h_{xx}\hfill \\ \multicolumn{1}{c}{g_{xx}}\end{array}\right\}=\sum _{j=1}^{N}g({\phi }_j)\sin {\phi }_j(\left\{\begin{array}{c}\hfill K_t\hfill \\ \multicolumn{1}{c}{\mu }\end{array}\right\}\cos {\phi }_j+\left\{\begin{array}{c}\hfill K_r\hfill \\ \multicolumn{1}{c}{1}\end{array}\right\}\sin {\phi }_j)\hfill \\ \left\{\begin{array}{c}\hfill h_{\mathrm{xy}}\hfill \\ \multicolumn{1}{c}{g_{\mathrm{xy}}}\end{array}\right\}=\sum _{j=1}^{N}g({\phi }_j)\cos {\phi }_j(\left\{\begin{array}{c}\hfill K_t\hfill \\ \multicolumn{1}{c}{\mu }\end{array}\right\}\cos {\phi }_j+\left\{\begin{array}{c}\hfill K_r\hfill \\ \multicolumn{1}{c}{1}\end{array}\right\}\sin {\phi }_j)\hfill \\ \left\{\begin{array}{c}\hfill h_{yx}\hfill \\ \multicolumn{1}{c}{g_{yx}}\end{array}\right\}=\sum _{j=1}^{N}g({\phi }_j)\sin {\phi }_j(\left\{\begin{array}{c}\hfill K_r\hfill \\ \multicolumn{1}{c}{1}\end{array}\right\}\cos {\phi }_j-\left\{\begin{array}{c}\hfill K_t\hfill \\ \multicolumn{1}{c}{\mu }\end{array}\right\}\sin {\phi }_j)\hfill \\ \left\{\begin{array}{c}\hfill h_{yy}\hfill \\ \multicolumn{1}{c}{g_{yy}}\end{array}\right\}=\sum _{j=1}^{N}g({\phi }_j)\cos {\phi }_j(\left\{\begin{array}{c}\hfill K_r\hfill \\ \multicolumn{1}{c}{1}\end{array}\right\}\cos {\phi }_j-\left\{\begin{array}{c}\hfill K_t\hfill \\ \multicolumn{1}{c}{\mu }\end{array}\right\}\sin {\phi }_j)\hfill \end{array}$$

(8)

The governing equation of this extended 2-DOF dynamic model can be formulated as a DDE, with time-varying coefficients

$$\begin{array}{c}\mathbf{M}\stackrel{··}{\mathbf{q}}\left(t\right)+\mathbf{C}\stackrel{·}{\mathbf{q}}\left(t\right)+\mathbf{K}\mathbf{q}\left(t\right)\hfill \\ =a_p\mathbf{H}(\left[\begin{array}{c}\hfill f_t\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]+\mathbf{q}\left(t\right)-\mathbf{q}(t-T))+C_{eq}\mathbf{G}\stackrel{·}{\mathbf{q}}\left(t\right)\hfill \\ \mathbf{M}=\left[\begin{array}{cc}\hfill m_{\mathrm{xx}}\hfill & \hfill m_{\mathrm{yx}}\hfill \\ \multicolumn{1}{c}{m_{\mathrm{xy}}}& \hfill m_{\mathrm{yy}}\hfill \end{array}\right],\mathbf{K}=\left[\begin{array}{cc}\hfill k_{\mathrm{xx}}\hfill & \hfill k_{\mathrm{yx}}\hfill \\ \multicolumn{1}{c}{k_{\mathrm{xy}}}& \hfill k_{\mathrm{yy}}\hfill \end{array}\right],\hfill \\ \mathbf{C}=\left[\begin{array}{cc}\hfill c_{\mathrm{xx}}\hfill & \hfill c_{\mathrm{yx}}\hfill \\ \multicolumn{1}{c}{c_{\mathrm{xy}}}& \hfill c_{\mathrm{yy}}\hfill \end{array}\right]\hfill \end{array}$$

(9)

where M, C, and K are the modal mass matrix, modal damping matrix, and modal stiffness matrix, respectively. As noted, the cross terms are nonzero, which indicates the mode coupling is taken into account. The subscript “xy” denotes that the response in the X direction is excited by the force in the Y direction. The second-order differential equations of q(t) is derived from equation (9) in the following general form

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{··}{\mathbf{q}}\left(t\right)=(a_p{\mathbf{M}}^{-1}\mathbf{H}-{\mathbf{M}}^{-1}\mathbf{K})\mathbf{q}\left(t\right)}\\ +(C_{eq}{\mathbf{M}}^{-1}\mathbf{G}-{\mathbf{M}}^{-1}\mathbf{C})\stackrel{·}{\mathbf{q}}\left(t\right)\hfill \\ -a_p{\mathbf{M}}^{-1}\mathbf{H}\mathbf{q}(t-T)+a_p{\mathbf{M}}^{-1}\mathbf{H}\left[\begin{array}{c}\hfill f_t\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]\hfill \end{array}$$

(10)

Then, let $\mathbf{X}(t)=[\mathbf{q}(t)\stackrel{·}{\mathbf{q}}(t){]}^T=[x(t)y(t)\stackrel{·}{x}(t)\stackrel{·}{y}(t){]}^T$ . Through algebraic transformations with equation (10), equation (9) is rewritten in the state-space form

$$\stackrel{·}{\mathbf{X}}\left(t\right)=\mathbf{A}\left(t\right)\mathbf{X}\left(t\right)+\mathbf{B}\left(t\right)\mathbf{X}(t-T)+\mathbf{f}\left(t\right)$$

(11)

where

$$\begin{array}{c}\mathbf{A}\left(t\right)=\left[\begin{array}{cc}\hfill 0\hfill & \hfill \mathbf{I}\hfill \\ \multicolumn{1}{c}{a_p{\mathbf{M}}^{-1}\mathbf{H}-{\mathbf{M}}^{-1}\mathbf{K}}& \hfill C_{eq}{\mathbf{M}}^{-1}\mathbf{G}-{\mathbf{M}}^{-1}\mathbf{C}\hfill \end{array}\right]\hfill \\ \mathbf{B}\left(t\right)=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{-a_p{\mathbf{M}}^{-1}\mathbf{H}}& \hfill 0\hfill \end{array}\right]\hfill \\ \mathbf{f}\left(t\right)=\left[\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{a_p{\mathbf{f}}_0\left(t\right)}\end{array}\right],{\mathbf{f}}_0\left(t\right)={\mathbf{M}}^{-1}\left[\begin{array}{c}\hfill -h_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{-h_{\mathrm{yx}}}\end{array}\right]f_t\hfill \end{array}$$

(12)

Milling stability and SLE analysis

Second-order SDM

Without loss of generality, the extended dynamic model in the state-space form can be solved by the SDM. In order to keep high computation accuracy, a second-order SDM is proposed to simultaneously satisfy SLD and SLE predictions. Regarding initial condition ${\mathbf{X}}_0=\mathbf{X}(0)$ , the general solution of equation (11) can be expressed in the integration form

$$\mathbf{X}\left(t\right)=e^{\mathbf{A}t}{\mathbf{X}}_0+\underset{0}{\overset{t}{\int }}{e}^{\mathbf{A}(t-s)}[\mathbf{B}\left(s\right)\mathbf{X}(s-T)+\mathbf{f}\left(s\right)]\mathrm{ds}$$

(13)

First, the time period T is discretized into $m$ small time intervals with each interval length $\tau =T/m$ . Equation (13) on the time interval [kτ, (k + 1)τ] can be formulated as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{X}}_{k+1}=e^{\mathbf{A}\tau }{\mathbf{X}}_k+\underset{0}{\overset{\tau }{\int }}{e}^{\mathbf{A}s}\mathbf{f}(k\tau +\tau -s)\mathrm{ds}}\\ +\underset{0}{\overset{\tau }{\int }}{e}^{\mathbf{A}s}\mathbf{B}(k\tau +\tau -s)\mathbf{X}(k\tau +\tau -s-T)ds\hfill \end{array}$$

(14)

where f(kτ + τ − s) and B(kτ + τ − s) are the periodic coefficient terms. X(kτ + τ − s − T) is the time-delay term. Then, the periodic coefficient terms are approximated by linear interpolation as follows

$$\begin{array}{c}\mathbf{f}(k\tau +\tau -s)={\mathbf{f}}_{k+1}+\frac{({\mathbf{f}}_k-{\mathbf{f}}_{k+1})s}{\tau }\hfill \\ \mathbf{B}(k\tau +\tau -s)={\mathbf{B}}_{k+1}+\frac{({\mathbf{B}}_k-{\mathbf{B}}_{k+1})s}{\tau }\hfill \end{array}$$

(15)

For the sake of improving the approximation accuracy, time-delay term is approximated using the following second-order Newton interpolation polynomial

$$\begin{array}{c}\mathbf{X}(k\tau +\tau -s-T)=(\frac{s^2}{2{\tau }^2}-\frac{s}{2\tau }){\mathbf{X}}_{k+2-m}\hfill \\ +(1-\frac{s^2}{{\tau }^2}){\mathbf{X}}_{k+1-m}+(\frac{s^2}{2{\tau }^2}+\frac{s}{2\tau }){\mathbf{X}}_{k-m}\hfill \end{array}$$

(16)

Through further algebraic transformations with equations (15) and (16), equation (14) can be transformed as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{X}}_{k+1}={\mathit{F}}_0{\mathbf{X}}_k+{\mathit{F}}_{m-2}{\mathbf{X}}_{k+2-m}}\\ +{\mathit{F}}_{m-1}{\mathbf{X}}_{k+1-m}+{\mathit{F}}_m{\mathbf{X}}_{k-m}+{\mathbf{D}}_k\hfill \end{array}$$

(17)

where

$$\begin{array}{c}{\mathit{F}}_0={\mathbf{\Phi }}_0\hfill \\ {\mathit{F}}_{m-2}=\frac{({\mathbf{\Phi }}_4-\tau {\mathbf{\Phi }}_3){\mathbf{B}}_k+(2\tau {\mathbf{\Phi }}_3-{\mathbf{\Phi }}_4-{\tau }^2{\mathbf{\Phi }}_2){\mathbf{B}}_{k+1}}{2{\tau }^3}\hfill \\ {\mathit{F}}_{m-1}=\frac{({\tau }^2{\mathbf{\Phi }}_2-{\mathbf{\Phi }}_4){\mathbf{B}}_k+({\mathbf{\Phi }}_4-\tau {\mathbf{\Phi }}_3-{\tau }^2{\mathbf{\Phi }}_2+{\tau }^3{\mathbf{\Phi }}_1){\mathbf{B}}_{k+1}}{{\tau }^3}\hfill \\ {\mathit{F}}_m=\frac{({\mathbf{\Phi }}_4+\tau {\mathbf{\Phi }}_3){\mathbf{B}}_k+({\tau }^2{\mathbf{\Phi }}_2-{\mathbf{\Phi }}_4){\mathbf{B}}_{k+1}}{2{\tau }^3}\hfill \\ {\mathbf{D}}_k={\mathbf{\Phi }}_1{\mathbf{f}}_{k+1}+\frac{{\mathbf{\Phi }}_2({\mathbf{f}}_k-{\mathbf{f}}_{k+1})}{\tau }\hfill \end{array}$$

(18)

and

$$\begin{array}{c}{\mathbf{\Phi }}_0=e^{\mathbf{A}\tau },{\mathbf{\Phi }}_1=\underset{0}{\overset{\tau }{\int }}{e}^{\mathbf{A}s}\mathrm{ds},{\mathbf{\Phi }}_2=\underset{0}{\overset{\tau }{\int }}s{e}^{\mathbf{A}s}\mathrm{ds}\hfill \\ {\mathbf{\Phi }}_3=\underset{0}{\overset{\tau }{\int }}{s}^2{e}^{\mathbf{A}s}\mathrm{ds},{\mathbf{\Phi }}_4=\underset{0}{\overset{\tau }{\int }}{s}^3{e}^{\mathbf{A}s}\mathrm{ds}\hfill \end{array}$$

(19)

Thus, according to the expanded state vector ${\mathbf{Z}}_k=[{\mathbf{X}}_k,{\mathbf{X}}_{k-1},\dots ,{\mathbf{X}}_{k+1-m},{\mathbf{X}}_{k-m}{]}^T$ , equation (17) can be rewritten as follows

$${\mathbf{Z}}_{k+1}={\mathbf{Q}}_k{\mathbf{Z}}_k+\mathbf{E}{\mathbf{D}}_k$$

(20)

where

$$\begin{array}{c}\mathbf{E}={[\mathbf{I},0,\dots ,0]}^T\hfill \\ {\mathbf{Q}}_k=\left[\begin{array}{cccccc}\hfill {\mathit{F}}_0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill {\mathit{F}}_{m-2}\hfill & \hfill {\mathit{F}}_{m-1}\hfill & \hfill {\mathit{F}}_m\hfill \\ \multicolumn{1}{c}{\mathbf{I}}& \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill \mathbf{I}\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill \cdots \hfill & \hfill \mathbf{I}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill \mathbf{I}\hfill & \hfill 0\hfill \end{array}\right]\hfill \end{array}$$

(21)

For one periodic time interval, the SLD can be finally determined using the Floquet theorem with the transition relation constructed as follows

$${\mathbf{Z}}_m=\mathbf{Q}{\mathbf{Z}}_0+\mathbf{D}$$

(22)

where

$$\begin{array}{c}\mathbf{Q}={\mathbf{Q}}_{m-1}{\mathbf{Q}}_{m-2}\dots {\mathbf{Q}}_1{\mathbf{Q}}_0\hfill \\ \mathbf{D}=\mathbf{E}{\mathbf{D}}_{m-1}+\sum _{i=0}^{m-2}({\mathbf{Q}}_{m-1}{\mathbf{Q}}_{m-2}\dots {\mathbf{Q}}_{i+1}\mathbf{E}{\mathbf{D}}_i)\hfill \end{array}$$

(23)

For the multiperiod case, the fixed point theory is utilized and the transition relation is described as follows

$${\mathbf{Z}}_n=\mathbf{Q}{\mathbf{Z}}_{n-1}+\mathbf{D}$$

(24)

where

$${\mathbf{Z}}_n={\mathbf{Z}}_{n-1}={\mathit{Z}}^{*}$$

(25)

Then, the steady-state coefficients Z * can be obtained in the following form

$${\mathit{Z}}^{*}={(\mathbf{I}-\mathbf{Q})}^{-1}\mathbf{D}$$

(26)

Finally, the SLE, which corresponds to the entry region for up-milling and exit region for down-milling, can be determined using the steady-state coefficients Z *.

Rate of convergence analysis

As a vital evaluation indicator, rate of convergence is widely utilized in recognized algorithms. In this section, the rate of convergence of the proposed method is investigated through the comparison with first-order SDM ¹⁶ and second-order FDM. ¹⁷ Using the error estimation method, the local discretization errors of the first-order SDM and second-order FDM are both proved as o(τ ³ ). Following the same way, the local discretization error of the proposed method is determined as o(τ ³ ).

In order to achieve a more in-depth analysis of the rate of convergence, a 1-DOF milling system in the study by Insperger and Stépán ¹⁵ is utilized here as an example, where the cutting conditions and modal parameters are listed in Table 1. Computer programs of the first-order SDM, the second-order FDM, and the proposed method are all written in MATLAB R2014a, and run on a personal computer (PC) (Intel^® Core™ i7-4790K CPU, 4.00 GHz, 16 GB). Through calculation, the convergence rate comparisons of above three algorithms at large and low radial immersions are illustrated in Figures 2 and 3, respectively. Note that the label $||\mu |-|{\mu }_0||$ means the deviation between the approximate critical eigenvalue $|\mu |$ and the exact one $|{\mu }_0|$ . In terms of accuracy, $|{\mu }_0|$ , as reference value, is determined by the first-order SDM with a large discretization parameter m = 600. For large radial immersions milling in Figure 2, the radial immersion ratio is a_D = 1. The spindle speed is chosen as Ω = 5000 r/min. The axial cutting depths are chosen as a_p = 0.2 and 1 mm, respectively. It can be seen clearly that the proposed method converges faster than the first-order SDM and the second-order FDM. For low radial immersion milling in Figure 3, the radial immersion ratio is a_D = 0.1. The spindle speed is chosen as Ω = 3000 r/min. The axial cutting depths are chosen as a_p = 1 and 2 mm, respectively. Although the proposed method suffers oscillations due to the loss-of-contact effect, this method converges slightly faster than the first-order SDM and the second-order FDM. Generally, the proposed method has a satisfactory convergence rate both at large and low radial immersions.

Table 1.

Cutting conditions and modal parameters

Cutting conditions	Milling mode	Down-milling
	Number of teeth	2
	Tangential shearing force coefficients	600 MPa
	Radial shearing force coefficients	200 MPa
Modal parameters	Natural frequency	922 Hz
	Modal mass	0.03993 kg
	Damping ratio	0.011

Figure 2.

Convergence rate comparisons of the first-order SDM, the second-order FDM, and the proposed method at large radial immersions (radial immersion ratio a_D = 1).

Figure 3.

Convergence rate comparisons of the first-order SDM, the second-order FDM, and the proposed method at low radial immersions (radial immersion ratio a_D = 0.1).

Experimental verification

In order to verify the veracity of the extended dynamic milling model, the practical experiments are conducted. Prior to calculating the SLD and SLE with the proposed method, the cutting force coefficients and modal parameters must be determined. Later on, almost all the experimental efforts focus on verifying the accuracy of prediction. In light of the accuracy concerns, all the experiments are carried out on the same 5-Axis CNC Vertical Machining Center (DMU 50 DECKEL MAHO).

Cutting force coefficients’ identification

As mentioned in dynamic model, the cutting forces are attributed to two physical mechanisms: the chip shearing and process damping effect. The process damping coefficients mainly depend on the material of the workpiece. According to Tunç and Budak, ³⁶ the specific indentation force coefficient and the coulomb friction coefficient corresponding to titanium alloy (Ti-6Al-4V) are identified as K_sp = 3×10⁴ N/mm³ and µ = 0.3. As for the shearing force coefficients K_t and K_r, the mechanistic approach ⁴⁰ is put into use. Considering the effect of process damping at low spindle speed, shearing force coefficients are determined at high spindle speed. A helical solid carbide end mill with four teeth, 10 mm diameter, 30° helix angle, and 60 mm overhang length is used during this experiment. The workpiece material is a titanium alloy block and it is machined in full-immersion conditions. The experimental configuration of shearing force coefficients is as sketched in Figure 4, where the experiment systems are made up of Kistler 9257B dynamometer, Kistler 5070A charge amplifier, Kistler 5697A1 data acquisition, and DynoWare software on PC. With the spindle speed and axial cutting depth kept as constant, the feed per tooth f_t varies over a range of 0.02–0.10 mm/tooth in increments of 0.02. Averaging the measured cutting forces and implementing linear regression, the shearing forces coefficients can be calibrated as K_t = 2061.9 MPa, K_r = 636.6 MPa.

Figure 4.

Experimental configuration of shearing force coefficients.

Modal parameters’ identification

Modal impact experiment plays an increasingly essential role in the identification of mechanical structures’ modal parameters, which needs good engineering judgment and experience. Figure 5 shows the impact directions of modal impact experiment, in which the subscripts including x and y mean same as mentioned in dynamic model. Note that the impact experiment is conducted four times at each direction so as to reduce the uncertainty. Also shown in this figure is the schematic diagram of obtaining experimental frequency response function (FRF). Analog signal obtained from sensor is first filtered by anti-aliasing filter. Then, analog signal is converted to digital signal with the function of analog digital converter (ADC). Due to the unavoidable leakage of signal, window function is generally applied. The FRF can be acquired later via fast Fourier transform (FFT).

Figure 5.

(a) The impact directions of modal impact experiment and (b) the schematic diagram of obtaining experimental FRF.

The procedure of modal impact experiment is as follows. An impact hammer with a piezoelectric force transducer is used to excite the tool tip. The corresponding response is measured by an accelerometer, which is mounted at tool tip. Once the data are processed by signal conditioning instrument and data acquisition instrument, the corresponding FRF can be measured by a dynamic signal analysis system. Through further conversion from acceleration to displacement, the measured FRFs (H_xx, H_xy, H_yx, and H_yy) of the tool tip are shown in Figure 5. As can be observed, the most flexible natural frequencies of H_xx, H_xy, H_yx, and H_yy are 1547.3, 1550.2, 1563.7, and 1549.6 Hz, respectively, which are basically equal to 1550 Hz. The magnitudes of cross FRFs (H_xy and H_yx) are smaller than that of the direct FRFs (H_xx and H_yy). This phenomenon indicates that mode coupling indeed exists in the practical milling. The modal parameters are identified by peak picking approach. ² H_xx is selected for an example to illustrate briefly. Figure 6 shows the real and imaginary parts of H_xx. The dominant mode is observed at f_n = 1547.3 Hz which corresponds to the lowest point C having negative value $4.12\times {10}^{-5}$ in the imaginary part. Damping ratio ζ_xx and modal stiffness k_xx are determined by the frequencies of points A and B and the peak value of point C

$${\zeta }_{\mathrm{xx}}=\frac{f_B-f_A}{2f_n},k_{\mathrm{xx}}=\frac{-1}{2{\zeta }_{\mathrm{xx}}\mathrm{Imag}(\mathrm{point}C)}$$

(27)

Figure 6.

Measured frequency response functions.

Due to w = 2πf_n, the modal mass and modal damping can be computed by $m_{\mathrm{xx}}=k_{\mathrm{xx}}/w^2$ and c_xx = 2m_xxζ_xxw. Analogous to H_xx, the relatively accurate modal parameters can be estimated from H_xy, H_yx, and H_yy which are all listed in Table 2. The modal fitting results of H_xx, H_xy, H_yx, and H_yy are shown in Figure 7. As can be seen, the fitted curves agree well with the measured curves.

Table 2.

Modal parameters of the tool.

Modal parameters	Modal mass (kg)	Modal damping (Ns/m)	Modal stiffness (N/m)
xx	0.00618	2.499	$5.8415\times {10}^5$
xy	0.01795	7.518	$1.7032\times {10}^6$
yx	0.03920	16.715	$3.7840\times {10}^6$
yy	0.00609	2.360	$5.7731\times {10}^5$

Figure 7.

Modal fitting results of H_xx, H_yx, H_xy, and H_yy.

Stability verification

For the stability verification, this section first predicts the stability with up-milling cut for a rigid titanium alloy block. The radial cutting depth is a_e = 0.25 mm and the corresponding radial immersion ratio is a_D = 0.025. According to the identified modal parameters and cutting force coefficients, the predicted SLD based on the extended dynamic model is calculated by the proposed method with m = 80, as shown in Figure 8. Note that unstable regions are shaded while stable regions are left unshaded.

Figure 8.

Experimental results and predicted SLD based on the extended dynamic model.

For the sake of verifying the stability prediction, a series of stability experiments at constant feed rate of 0.1 mm/tooth are performed with different parameter combinations of spindle speed and axial cutting depth. The axial cutting depth increases gradually from 0.5 to 1.5 mm in increments of 0.25 mm, and the spindle speed varies from 1200 to 8000 r/min. Once again with Kistler 9257B dynamometer, the dynamic cutting force can be measured during the milling process. Based on the machined surface quality and the FFT of measured cutting forces, the experimental stability results are obtained as depicted in Figure 8. The crosses, triangles, and circles stand for chatter, marginal cut, and stable cut, separately. As can be observed, in spite of some small quantitative discrepancies, it should be emphasized that experimental results reach a satisfying correlation with the prediction. The slight discrepancies may be due to the neglect of run-out, which could be a subject for future research.

Besides, in order to make the analysis more targeted and representative, five typical parameter combinations of spindle speed and axial cutting depth, namely, A (5000 r/min, 1.5 mm), B (5000 r/min, 1.25 mm), C (5000 r/min, 1.0 mm), D (5000 r/min, 0.75 mm), and E (5000 r/min, 0.5 mm), are selected so as to illustrate the presence of chatter and stable cuts. Figure 9 shows the experimental system of stability experiment. The corresponding frequency spectrum of measured cutting forces and the machined surface quality is as shown in Figure 10. In the case of stable point E, the measured cutting forces are mainly dominated by the tooth passing frequency $(f_{\mathrm{tp}}=5000\times 4/60=333.3\mathrm{Hz})$ and harmonic frequencies (666, 1000, and 1333 Hz). Also, there are no apparent chatter marks appearing on the machined surface. In the case of unstable point C, the chatter frequency occurs in the vicinity of 1577 Hz, which is close to the natural frequency of the tool (about 1550 Hz). Apparent chatter marks appear on the machined surface. With the increase in axial cutting depth, the vibration grows and new chatter frequency occurs at 1911 Hz in the case of unstable points A and B. Note that the difference of two adjacent chatter frequencies basically equals tooth passing frequency. Finally, for the marginal point D near the stability boundary, slightly unstable status occurs.

Figure 9.

Experimental system of stability experiment.

Figure 10.

(a) Frequency spectrum of measured cutting forces and (b) the machined surface quality.

After experimental stability verification, the effects of mode coupling and process damping are further studied. Figure 11 shows the predicted SLDs based on three different dynamic models which correspond to the solid line, dashed line, and dotted line. The radial immersion ratio is kept as a_D = 0.025. Solid line means the extended dynamic model with mode coupling and process damping. Dashed line represents the dynamic model with mode coupling and without process damping. Dotted line denotes traditional dynamic model without mode coupling and process damping. Through the comparison of the dashed line and dotted line, we observe that the stability regions are increased within the whole spindle speeds range due to the mode coupling effect. The reason for this phenomenon may be that the transferred energy is scattered due to mode coupling effect, although the major energy still focuses on regeneration effect. Thus, regeneration needs more energy to yield chatter and the critical stability limit increases correspondingly. Another discovery is that the critical stability limit of the solid line SLD has an obvious trend of enhancement at low spindle speeds. In contrast, the critical stability limit of dashed line SLD is always kept as constant. Additionally, the SLDs with and without process damping are basically consistent at high spindle speeds. This phenomenon indicates that process damping plays a vital role in stability improvement mainly at low spindle speeds.

Figure 11.

The predicted SLDs with different dynamic models.

SLE verification

For SLE verification, both experimental and theoretical values need to be determined. In view of high spindle speed sensitivity, a series of cutting experiments are performed during the spindle speed range 7300–8000 r/min in increments of 100 r/min, and the axial cutting depth is chosen as a constant value of a_p = 0.75 mm. The other milling condition is identical with the one in stability verification. The schematic diagram of SLE experiments is shown in Figure 12, where the original thickness of workpiece is L = 20 mm and theoretical radial cutting depth is a_e = 0.25 mm. After cutting experiment, 10 measured points with equal space are selected along the machined surface at the position of half axial cutting depth. The coordinate values y of all measured points are measured by a coordinate measuring machine (Zeiss Prismo), which are depicted in Figure 13. Then, the experimental SLE can be calculated as follows

$$\mathrm{SLE}=\{\begin{array}{c}\hfill \max |y-(L-a_e)|\mathrm{when}y\prime <(L-a_e)\hfill \\ \multicolumn{1}{c}{-\max |y-(L-a_e)|\mathrm{when}y\prime >(L-a_e)}\end{array}$$

(28)

where y′ denotes the y satisfying max $|y-(L-a_e)|$ . Note that, the positive value of SLE yields an overcut surface for up-milling in this experimental verification, and the negative value yields an undercut surface.

Figure 12.

The schematic diagram of SLE experiments.

Figure 13.

The measured coordinate values y of all measured points.

The experimental SLE is shown in Figure 14. Also shown in this figure are the predicted SLEs based on three different dynamic models which correspond to the solid line, dash-dot line, and dashed line. One should notice that the dashed line only shows the stable area during the spindle speed range 7300–7873 r/min, according to Figure 11, while for solid line and dash–dot line, stable condition can be satisfied over the whole experimental spindle speeds. Through the comparison, the solid line, which represents the predicted SLE based on the extended dynamic model with mode coupling and process damping, correlates better with the experimental SLE. This phenomenon indicates that the proposed method for SLE prediction based on the extended dynamic model can also reach a high accuracy. The possible reason for slight deviation may be tool run-out, geometric errors in the machine axis, and thermal errors. Additionally, an interesting phenomenon can be seen that the SLE is large and changes rapidly near 1550 × 60/4/3 = 7750 r/min. This is because the system resonance is excited when the tooth passing frequency and its harmonics are near the natural frequency about 1550 Hz, which results in the high sensitivity of SLE. Moreover, for up-milling, the positive SLE yields an overcut surface, while negative SLE is opposite. In actual milling process, overcut is more undesirable than undercut. Because in the case of undercut, the surface quality of the workpiece can be improved by additional material removal operations; however, in the case of overcut, the error is irreversible. Thus, a positive SLE should be avoided when selecting suitable spindle speed and axial cutting depth in up-milling.

Figure 14.

The comparison of the predicted and experimental SLE.

Conclusion

In this article, the SLD and SLE based on the extended dynamic milling model with mode coupling and process damping are predicted simultaneously by a second-order SDM. Experimental results verify the veracity of predictions based on the extended dynamic milling model. The following conclusions are drawn from this work:

• Based on traditional 2-DOF dynamic milling process with regeneration, an extended dynamic model including mode coupling and process damping is established as DDE, with time-varying coefficients. Through transformation, the extended dynamic model in the state-space form is obtained.

• A second-order SDM is proposed to solve the extended dynamic model in the state-space form. The SLD and SLE can be obtained simultaneously via Floquet theory and fixed point theory, respectively. Through the comparison with the first-order SDM and the second-order FDM, the proposed method shows a faster rate of convergence.

• The cutting force coefficients and modal parameters are experimentally determined. The SLD and SLE based on the extended dynamic model are predicted by the proposed second-order SDM. Via experimental verification, the experimental results are in good correlation with the predictions. This phenomenon indicates that the SLD and SLE prediction based on the extended dynamic model can achieve higher prediction accuracy compared with traditional dynamic model.

• The direct and cross FRFs are measured by modal impact experiment. The magnitudes of cross FRFs are smaller than that of the direct FRFs, which demonstrates that the mode coupling indeed exists in the practical milling. Through experiments and simulation, it is shown that mode coupling seems to increase the whole stability region. As for process damping, it mainly plays a role at low spindle speeds and affects little at high spindle speeds.