# An Efficient Third-Order Full-Discretization Method for Prediction of Regenerative Chatter Stability in Milling

Authors: Chao Huang; Wen-An Yang; Xulin Cai; Weichao Liu; YouPeng You
Journal: Shock and Vibration
 (2020)

---

## Abstract

The prediction of regenerative chatter stability has long been recognized as an important issue of concern in the field of machining community because it limits metal removal rate below the machine’s capacity and hence reduces the productivity of the machine. Various full-discretization methods have been designed for predicting regenerative chatter stability. The main problem of such methods is that they can predict the regenerative chatter stability but do not efficiently determine stability lobe diagrams (SLDs). Using third-order Newton interpolation and third-order Hermite interpolation techniques, this study proposes a straightforward and effective third-order full-discretization method (called NI-HI-3rdFDM) to predict the regenerative chatter stability in milling operations. Experimental results using simulation show that the proposed NI-HI-3rdFDM can not only efficiently predict the regenerative chatter stability but also accurately identify the SLD. The comparison results also indicate that the proposed NI-HI-3rdFDM is very much more accurate than that of other existing methods for predicting the regenerative chatter stability in milling operations. A demonstrative experimental verification is provided to illustrate the usage of the proposed NI-HI-3rdFDM to regenerative chatter stability prediction. The feature of accurate computing makes the proposed NI-HI-3rdFDM more adaptable to a dynamic milling scenario, in which a computationally efficient and accurate chatter stability method is required.


---

## Body

## 1. Introduction

The control of cutting vibration and machining instability plays a significant role in high quality and high-efficiency machining of a workpiece, such as titanium alloy and high-temperature alloy thin-walled components. The vibrations existing in the machining processes consist of free vibration, forced vibration, and self-excited vibration, where the last one can be further divided into three main categories: regenerative chatter, model coupling chatter, and frictional chatter [1, 2]. With respect to milling operations, it has been proved that the main factor that influences the machining stability is regenerative chatter [3], which may limit metal removal rate below the machine’s capacity and hence reduce the productivity of the machine. Therefore, accurate and efficient regenerative chatter stability prediction is a critical task in milling operations for determining the potential regions of stable machining.The main advance in the study of machining dynamics in the past few years has come in recognizing and interpreting a variety of different delay-differential equations (DDEs) that arise in the milling processes when considering the regenerative effects [4, 5]. Stability lobe diagrams (SLDs) obtained with DDEs [6, 7] are the most widely applied tools for revealing the boundary between stable (i.e., chatter-free) and unstable (i.e., with chatter) machining and locating their possible optimal chatter-free parameters (e.g., spindle speed and depth of cut). The theoretical advances in numerical analysis methods and the use of on-line computers for solving DDEs in machining process modeling have led to increased interest in creating SLDs for stable machining parameter selection. These techniques are referred to as frequency-time methods (e.g., zero order analytical (ZOA) [8] and multifrequency (MF) [9]) and time-domain methods (e.g., temporal finite element analysis (TFEA) [10], first-order semidiscretization method (1stSDM) [11], updated semidiscretization method (USDM) [12], first-order full-discretization method (1stFDM) [13], second-order full-discretization method (2ndFDM) [14], third-order full-discretization method (3rdFDM) [15], third-order least squares approximated full-discretization method (3rdLSAFDM) [16], hyper-third-order full-discretization method (HTOFDM) [17], numerical integration method (NIM) [18], improved numerical integration method (INIM) [19], complete discretization scheme (CDS) [20], Runge–Kutta methods (RKM) [21], Runge–Kutta-based complete discretization method (RK-CDM) [22], third-order Hermite approximation method (3rdHAM) [23], third-order Hermite–Newton approximation method (3rdH-NAM) [24], high-order full-discretization method (HFDM) [25], and third-order updated full-discretization method (3rdUFDM) [26]). ZOA methods that expand time varying coefficients of dynamic milling equations into Fourier series and truncate the series to include only the average component have been applied in regenerative chatter stability prediction tasks. Although it was proved that the ZOA can perform well in regenerative chatter stability prediction tasks, not readily applicable to low radial immersion milling is a common problem for ZOAs in regenerative chatter stability prediction. The MF method extended the study of ZOA to consider the harmonics of the tooth-passing frequencies for low radial immersion machining processes. Generally, however, this method is limited to harmonically excited systems and as the harmonics in the expansion of the nonlinear force becomes increasingly high, there are growing concerns that time cost could become a threat.In recent years, intensive research efforts have been made in the use of time-domain techniques as an effective approach to regenerative chatter stability prediction. Bayly et al. [10] proposed a TFEA method to predicting milling stability. A comparative study indicated that the performance of their proposed method was superior to that of traditional frequency-domain methods with regard to computation accuracy. Insperger et al. [11] proposed the1stSDM, where the latter had higher computation efficiency and computation accuracy. Long and Balachandran [12] developed the USDM for a milling dynamics model with consideration of the loss-of-contact and time-delay effects and indicated that period-doubling bifurcations of periodic orbits can occur in low-immersion operations. Ding et al. [13] proposed the 1stFDM for milling stability prediction. A comparative study indicated that the proposed 1stFDM was superior to the 1stSDM with regard to computation efficiency. Owning to its wonderful computation accuracy, various variants of the 1stFDM [14–17] such as the 2ndFDM, 3rdFDM, 3rdLSAFDM, and HTOFDM were subsequently developed for milling stability prediction. Ding et al. [18] employed NIM to recognize milling stability boundaries and showed that their NIM outperformed the traditional SDM and TFEA method with regard to convergence rate and computation efficiency. Liang et al. [19] extended the study of the NIM to develop an INIM for creating SLD for milling operations with time-varying delay. Li et al. [20] proposed the CDS for milling stability prediction by discretizing the delay time term, state term, and periodic-coefficient matrices simultaneously. A comparative study indicated that their CDS could outperform the SDM and FDM with regard to computation efficiency. Niu et al. [21] applied the fourth-order Runge–Kutta method to handle the Duhamel term of DDEs for milling stability prediction. Their approach was superior to 1stSDM in terms of convergence rate and computation efficiency. Li et al. [22] presented the RK-CDM to predict the chatter stability for milling processes. An extensive comparison based on convergence rate and computation efficiency as the performance index showed that the proposed RK-CDM could outperform the SDM and CDS in predicting milling stability. However, nearly all the abovementioned methods in this research area have used the coordinates of each node as interpolation conditions; the curve was maintained relatively flat when interpolating with lower-order functions while it would not be maintained smooth at the joints of the piecewise curves. The result might as well not provide the near-optimum SLD. On the contrary, these methods might lead to a poor overall machining parameters selection.Recently, the Hermite interpolation polynomial is one of many polynomial interpolation algorithms that have drawn increasing attention for their attractive properties and excellent error estimates on a wide range of numerical summation problems. Hermite interpolation polynomials not only define the coordinate value of each node but also specify the derivative of the curve at each node, which has been shown to be superior to the Newton interpolation and Lagrange interpolation in keeping the stability of the curve. While the Hermite interpolation polynomial solution may be the globally stable optimum, current FDM methods tend to fall into a local unstable solution. The unsmooth connection problem is unlikely to occur with Hermite interpolation polynomial. Therefore, Hermite interpolation polynomial can provide excellent performance for milling stability prediction problems.To the best of the authors’ knowledge, only three studies in the literature have focused on the Hermite interpolation polynomial-based machining stability prediction for milling operations. Liu et al. [23] predicted the machining stability in milling operations using the 3rdHAM and compared their results with those obtained using the 1stSDM and 2ndFDM. In the efficient FDM, linear interpolation and Hermite interpolation were approximate to the state item and time-delay item, respectively. Compared with the results obtained using the 1stSDM and 2ndFDM, their FDM reduced computation time by almost 53% and 10%, respectively. Ji et al. [24] used the 3rdH-NAM to resolve the problem of machining stability prediction for milling processes. In the 3rdH-NAM, third-order Hermite interpolation and third-order Newton interpolation were used to approximate the state item and time-delay item, respectively. Simulations showed that their approach was superior to Guo’s 3rdFDM and Liu’s 3rdHAM in terms of the accuracy of identifying the milling stability boundaries. Liu et al. [25] proposed a HFDM to resolve the problem of machining stability prediction for turning processes and showed that their HFDM improved calculation efficiency by about 50% compared to the very efficient 2ndFDM. In the HFDM, Lagrange interpolation and Hermite interpolation were used to approximate the state values and delayed state values, respectively. Although the results of these studies were promising, the problem of stability prediction accuracy has still remained unsolved, and more accurate trials are needed.In the present research, based on third-order Newton interpolation polynomial and third-order Hermite interpolation polynomial techniques, this study proposes a straightforward and effective third-order full-discretization Method (called NI-HI-3rdFDM) to predict the regenerative chatter stability in milling operations. Experimental results using simulation show that the proposed NI-HI-3rdFDM can not only efficiently predict the regenerative chatter stability but also accurately identify the SLD. The comparison results also indicate that the proposed NI-HI-3rdFDM is very much more accurate than that of other existing methods for predicting the regenerative chatter stability in milling operations. A demonstrative experimental verification is provided to illustrate the usage of the proposed NI-HI-3rdFDM approach to regenerative chatter stability prediction. The feature of accurate computing makes the proposed NI-HI-3rdFDM more adaptable to a dynamic milling scenario, in which a computationally efficient and accurate method is required.The study is organized as follows. Section2 describes the mathematical model of milling dynamics adopted in this study. The proposed NI-HI-3rdFDM for solving the mathematical model of milling dynamics for regenerative chatter stability prediction is described in Section 3. Section 4 presents the performance comparisons of the proposed NI-HI-3rdFDM with those of existing methods in the literature. Experimental results are given in Section 5, and these indicate that the proposed NI-HI-3rdFDM is capable of efficiently creating superior SLDs for accurate and efficient regenerative chatter stability prediction. Finally, Section 6 is devoted to conclusions and directions for future research.


## 2. Mathematical Model of Milling Dynamics

It is generally accepted that the behavior characteristic of a milling operation is just as much the product of DDE as the morphology. Thus, by considering regenerative effects, the behavior which can be adequately explained by the classical mathematical theory of milling dynamics is described in state-space form by the equation as follows:(1)x˙t=A0xt+Btxt−xt−T,where A0 is a constant matrix and represents the time-invariant nature of the system and B(t) is one periodic-coefficient matrix and satisfies Bt+T=Bt.The general response of equation (1) can thus be obtained via direct numerical integration as follows:(2)xt=eA0t−t0xt0+∫t0teA0t−ξBξxξ−xξ−Tdξ.In order to solve equation (2), by introducing a time interval τ, the time period T is equally divided into m small time intervals, and therefore, T=mτ. On each time interval t∈kτ,k+1τk=0,1,…,m, equation (2) can be rewritten as the following equation:(3)xt=eA0t−kτxkτ+∫kτteA0t−ξBξxξ−xξ−Tdξ.Letxk denote xkτ, then xk+1 can be got from equation (3) as the following equation:(4)xk+1=eA0τxk+∫0τeA0τ−δBδxδ−xδ−Tdδ,where δ=ξ−kτ, δ∈0,τ. In earlier investigations, it has been demonstrated that the performance of several approaches including NIM, CDS, RK-CDM, FDM, and its variants are very sensitive to interpolation polynomial for approximating the state item and time-delay item. In the current study, the state term x(δ) is interpolated by third-order Newton interpolation polynomial, while the time-delay term xδ−T is approximated by third-order Hermite interpolation polynomial. Specifically, on the time interval kτ,k+1τ, the nodal values xkτ+τ, xkτ, xkτ−τ, and xkτ−2τ(denoted as xk+1, xk, xk−1, and xk−2, respectively) are used to approximate xδ, and the nodal values xkτ−T, xkτ+τ−T, and xkτ+2τ−T (denoted as xk−m, xk+1−m and xk+2−m, respectively) are used to approximate xδ−T. In addition, because high order interpolation of Bδ has no apparent effect on improving the computation efficiency and computation accuracy of the dynamic milling system [27], linear interpolation is employed to approximate Bδ with the boundaries Bkτ,Bkτ+τ, which is denoted here as Bk,Bk+1.The state termxδ is determined by third-order Newton interpolation polynomial as follows:(5)xδ=a1xk+1+b1xk+c1xk−1+d1xk−2,where(6)a1=δ3+3τδ2+2τ2δ6τ3,b1=δ2+3τδ+2τ22τ2−δ3+3τδ2+2τ2δ2τ3,c1=δ+2ττ−δ2+3τδ+2τ2τ2+δ3+3τδ2+2τ2δ2τ3,d1=1+δ+2ττ+δ2+3τδ+2τ22τ2−δ3+3τδ2+2τ2δ6τ3.Using the derivative formula, the first time derivative of thex(t) at t=tk−m and t=tk+1−m could be got from(7)x˙k−m=xk+1−m−xk−mτ,x˙k+1−m=xk+2−m−xk+1−mτ.Then, the time-delay termxδ−T is obtained by third-order Hermite interpolation polynomial at the time span tk−m,tk+1−m and expressed as(8)xδ−T=a2xk−m+b2xk+1−m+c2xk+2−m,where(9)a2=τ3−τ2δ−τδ2+δ3τ3,b2=τ2δ+2τδ2−2δ3τ3,c2=δ3−τδ2τ3.Similarly, the periodic-coefficient term is achieved by linear interpolation at the time spankτ,k+1τ, resulting in(10)Bδ=δτBk+1+τ−δτBk.Substituting equations (5), (8), and (10) into equation (4) leads to(11)Hk,1−Ixk+1+H0+Ηk,0xk+Hk,−1xk−1+Hk,−2xk−2=Hk,mxk−m+Hk,m−1xk+1−m+Hk,m−2xk+2−m,where(12)H0=Φ0=eA0τ,Hk,1=Φ5+3τΦ4+2τ2Φ36τ4Bk+1+−Φ5−2τΦ4+τ2Φ3+2τ3Φ26τ4Bk,Hk,0=Φ5+4τΦ4+5τ2Φ3+2τ3Φ22τ4Bk+1+−Φ5−3τΦ4−τ2Φ3+3τ3Φ2+2τ4Φ12τ4Bk,Hk,−1=Φ5+τΦ4−2τ2Φ32τ4Bk+1+−Φ5+3τ2Φ3−2τ3Φ22τ4Bk,Hk,−2=−Φ5+τ2Φ36τ4Bk+1+Φ5−τΦ4−τ2Φ3+τ3Φ26τ4Bk,Hk,m−2=Φ5−τΦ4τ4Bk+1+−Φ5+2τΦ4−τ2Φ3τ4Bk,Hk,m−1=−2Φ5+2τΦ4+τ2Φ3τ4Bk+1+2Φ5−4τΦ4+τ2Φ3+τ3Φ2τ4Bk,Hk,m=Φ5−τΦ4−τ2Φ3+τ3Φ2τ4Bk+1+−Φ5+2τΦ4−2τ3Φ2+τ4Φ1τ4Bk,Φ1=∫0τe−A0δdδ,Φ2=∫0τe−A0δδdδ,Φ3=∫0τe−A0δδ2dδ,Φ4=∫0τe−A0δδ3dδ,Φ5=∫0τe−A0δδ4dδ.The matricesΦ1–Φ5 could be calculated by Φ0 and A0 as(13)Φ1=A0−1Φ0−I,(14)Φ2=A0−1Φ1−τI,(15)Φ3=A0−12Φ2−τ2I,(16)Φ4=A0−13Φ3−τ3I,(17)Φ5=A0−14Φ4−τ4I.Note that it is thus highly desirable to find a one-to-one discrete mapping between the current and previous period [27]; all the current period variables x (i.e., xk+1,xk,xk−1,xk−2) should be located in the left hand side of equation (11), while all the previous period variables x (i.e., xk+2−m,xk+1−m,xk−m) in the right hand side of equation (11). If not, corresponding terms need to be converted into a required range. In equation (11), when k=1, the left hand side variables xk−1 and xk−2 of equation (11) are equal to x0 and x−1 which are also presented as xm−m and xm−1−m in the previous period. Therefore, equation (11) is rewritten as(18)H1,1−Ix2+H0+Η1,0x1=H1,mx1−m+H1,m−1x2−m+H1,m−2x3−m−H1,−1xm−m−H1,−2xm−1−m.Whenk=2, the left hand side variable xk−2 of equation (11) is equal to x0 which is also equal to xm−m in the previous period. Then, equation (11) can be rewritten as(19)H2,1−Ix3+H0+Η2,0x2+H2,−1x1=H2,mx2−m+H2,m−1x3−m+H2,m−2x4−m−H2,−2xm−m.Whenk=m, the right hand side variable xk+2−m of equation (11) is equal to x2 which is in the current period. Equation (11) can be rewritten as(20)Hm,1−Ixm+1+H0+Ηm,0xm+Hm,−1xm−1+Hm,−2xm−2−Hm,m−2x2=Hm,mxm−m+Hm,m−1xm+1−m.According to equations (11) and (18)–(20), the discrete mapping is obtained as(21)D1x1x2⋮xmxm+1=D2x1−mx2−m⋮xm−mxm+1−m,where(22)D1=I000⋯00000H0+H1,0H1,1−I00⋯00000H2,−1H0+H2,0H2,1−I0⋯00000H3,−2⋮H3,−1⋮H0+H3,0⋮H3,1−I⋮⋯⋱0⋮0⋮0⋮0⋮0⋮0000⋯Hm−1,−2Hm−1,−1H0+Hm−1,0Hm−1,1−I00−Hm,m−200⋯0Hm,−2Hm,−1H0+Hm−1,0Hm−1,1−I,D2=00000⋯00IH1,mH1,m−1H1,m−200⋯−H1,−2−H1,−100H2,mH2,m−1H2,m−20⋯0−H2,−2000H3,mH3,m−1H3,m−2⋯000⋮⋮⋮⋮⋮⋱⋮⋮⋮00000⋯Hm−1,mHm−1,m−1Hm−1,m−200000⋯0Hm,mHm,m−1.The transition matrix of dynamic systemΦ is defined over one period T as(23)Φ=D1−1D2.Finally, if the maximum modulus for the eigenvalue of the transition matrixФ is the at most one, then conclude that the dynamic milling system is in stable (i.e., chatter-free) machining in accordance with Floquet theory [28]; otherwise, conclude that the dynamic milling system is in unstable (i.e., with chatter) machining in accordance with Floquet theory [28].


## 3. The Proposed Method for Milling Stability Analysis

The motion of a single DOF milling system with single discrete time delay can be expressed as a functional differential equation as follows:(24)q¨t+2ζωnq˙t+ωn2qt=−whtmqt−qt−T,where ζ, ωn, mt, w, and q(t) denote the relative damping, natural frequency, modal mass, depth of cut, and displacement vector of the cutter, respectively. Tooth-passing period T is equal to 60/NΩ, where N is the number of teeth and Ω is the spindle speed. The directional cutting force coefficient ht is one periodic function with period T and is defined as(25)ht=∑j=1NgϕjtsinϕjtKtcosϕjt+Knsinϕjt,where Kt and Kn are the tangential and normal cutting force coefficients, and ϕjt denotes the jth tooth’s angular position and is expressed as(26)ϕjt=2πΩ60t+2πj−1N.The window functiongϕjt is defined by(27)gϕjt=1,ϕst<ϕjt<ϕex,0,otherwise,where ϕst and ϕex are the start and exit angles, respectively. As for up-milling, ϕst=0 and ϕex=arcos1−2a/D; as for down-milling, ϕst=arcos2a/D−1 and ϕex=π, where a/D is the radial depth of the cut ratio [5].Letpt=mq˙t+mζωnqt and xt=qtpt, then equation (16) can be expressed as(28)x˙t=A0xt+Btxt−xt−T,where(29)A0=−ζωnmt−1mtζωn2−mtωn2−ζωn,Bt=00−wht0.The modal parameters relating to the dominant mode of the milling system are given in Table1. For more information about periodic DDEs, one can refer to the literature [5].Table 1
Modal parameters.
Number of teeth (N)Frequencyf (Hz)Damping (Ζ)Massmt (kg)Tangential cutting force coefficientKt (N/m2)Normal cutting force coefficientKn (N/m2)29220.0110.039936 × 1082 × 108


### 3.1. Convergence Rate

In order to evaluate the convergence rate of different time-domain methods such as the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM, the local error ||u| − |u0|| is expressed as the function of the discrete number. Noting that |u| is the maximum modulus of the eigenvalue of the transition matrix and |u0| is obtained by the 1stSDM [11] at m = 200. All the other parameters for the milling process are summarized in Table 2. All of the calculation programs are written in MATLAB 9.2 and run on a personal computer Intel® Core™ i5-7400 CPU and 8 GB memory.Table 2
Parameters for the milling process.
Milling typeRadial depth of cut ratioa/DSpindle speed Ω (rpm)Depth of cutw (mm)Up-milling150000.2, 0.5, 0.7, 1Figure1 presents the convergence rates of the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM in predicting the milling stability. Figure 1 indicates that the proposed NI-HI-3rdFDM can converge sooner than the NIM, 3rdHAM, and 3rdFDM can. Specifically, the local errors ||u| − |u0|| of the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM were 0.0155, 0.0084, 0.0233, and 0.6272 × 10−3, respectively, when Ω = 5000 rpm, w = 0.5 mm, and m = 40. The local error of the proposed NI-HI-3rdFDM decreased by 90%, 92.5%, and 97.3% in comparison to that of the NIM, 3rdHAM, and 3rdFDM, respectively. Although the proposed NI-HI-3rdFDM converged a little slower than 3rdHAM in case of w = 0.5 mm and |u0| = 1.0726, for other three cases, the proposed NI-HI-3rdFDM was fairly faster than 3rdHAM. This slow convergence rate can be explained by the fact the derivability of the state item was utilized to enhance the calculation precision in 3rdHAM.Figure 1
Convergence rate for the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM. (a)w = 0.2 mm, u0=0.8192; (b) w = 0.5 mm, u0=1.0726; (c) w = 0.7 mm, u0=1.2197; (d) w = 1 mm, u0=1.4040.
(a)(b)(c)(d)


### 3.2. Stability Lobe Diagrams

In order to demonstrate the computation accuracy and computation time performances of the proposed NI-HI-3rdFDM, the SLDs obtained with the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM form = 20, 30, and 40 are summarized in Table 3. Note that the machining operation is up-milling operation, the ratio of radial cutting depth and tool diameter is a/D = 1, the spindle speed Ω, and depth of cut w is set to be 5000 rpm to 10000 rpm and 0 mm to 4 mm over a 200×100 sized grid of parameters, respectively [14]. It is also worth noting that in Table 3 the stability charts obtained with the 1stSDM for m = 200 in red color were regarded as the criterion-referenced SLDs. As can be seen, the stability charts of the proposed NI-HI-3rdFDM were much closer to the criterion-referenced lobe diagrams than those of the NIM, 3rdHAM and 3rdFDM. Table 3 indicates that the proposed NI-HI-3rdFDM can also achieve better or at least comparable performance with existing methods such as the NIM, 3rdHAM, and 3rdFDM with respect to the computation time. Specifically, in the case of m = 30, the computation time of the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM were 10 s, 92 s, 38 s, and 23 s, respectively. In this case, the computation efficiency of the proposed NI-HI-3rdFDM increased by 75% and 39.5% compared with that of the 3rdHAM and 3rdFDM, respectively, except for the NIM. Although the computation time of the proposed NI-HI-3rdFDM was slightly greater than that of the NIM, the computation accuracy obtained with the proposed NI-HI-3rdFDM (corresponding computation time = 23 s) was even equal to that of the NIM when m = 60 (corresponding computation time = 37 s). Moreover, in the case of m = 40, the stability charts of the proposed NI-HI-3rdFDM were in good agreement with the criterion-referenced lobe. A visual examination of Table 3 also revealed that the proposed NI-HI-3rdFDM not only obtained successfully accurate milling stability prediction but also saved a lot of computation time.Table 3
SLD for the single DOF milling model using the NIM, 3rdHAM, 3rd FDM, and the proposed NI-HI-3rdFDM.
m = 20Computation time 6 sComputation time 55 sComputation time 22 sComputation time 15 sm = 30Computation time 10 sComputation time 92 sComputation time 38 sComputation time 23 sm = 40Computation time 21 sComputation time 125 sComputation time 54 sComputation time 37 sIn addition, the proposed NI-HI-3rdFDM can also be extended to a multifreedom milling system. A stability prediction of a two-DOF milling system was employed to evaluate the accurate and efficient performances of NI-HI-3rdFDM in predicting milling stability. Note that the same system parameters taken from the literature [5] were used. The SLDs calculated from up-milling operations with m = 40 and a/D = 1, 0.5, 0.1, and 0.05 are shown in Table 4. As can be seen in Table 4, the SLD of the proposed NI-HI-3rdFDM was relatively identical to that of the 1stFDM and the 1stSDM, but the computation time of the proposed NI-HI-3rdFDM was much less than that of the 1stFDM and the 1stSDM under the same milling parameters. The proposed NI-HI-3rdFDM achieved better and comparable efficiency compared with the 1stFDM and the 1stSDM. Especially for a/D = 0.05, the runtime of the proposed NI-HI-3rdFDM decreased from 200 s to 160 s when compared with the 1stFDM and from 1197 s to 160 s when compared with the 1stSDM.Table 4
SLD for the two-DOF milling model using the 1stSDM, 1stFDM, and the proposed NI-HI-3rdFDM.
Radial depth of cut ratio1stSDM [11]1stFDM [13]The proposed NI-HI-3rdFDMa/D = 1Computation time 1247 sComputation time 282 sComputation time 256 sa/D = 0.5Computation time 1189 sComputation time 243 sComputation time 190 sa/D = 0.1Computation time 1204 sComputation time 202 sComputation time 163 sa/D = 0.05Computation time 1197 sComputation time 200 sComputation time 160 s


## 3.1. Convergence Rate

In order to evaluate the convergence rate of different time-domain methods such as the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM, the local error ||u| − |u0|| is expressed as the function of the discrete number. Noting that |u| is the maximum modulus of the eigenvalue of the transition matrix and |u0| is obtained by the 1stSDM [11] at m = 200. All the other parameters for the milling process are summarized in Table 2. All of the calculation programs are written in MATLAB 9.2 and run on a personal computer Intel® Core™ i5-7400 CPU and 8 GB memory.Table 2
Parameters for the milling process.
Milling typeRadial depth of cut ratioa/DSpindle speed Ω (rpm)Depth of cutw (mm)Up-milling150000.2, 0.5, 0.7, 1Figure1 presents the convergence rates of the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM in predicting the milling stability. Figure 1 indicates that the proposed NI-HI-3rdFDM can converge sooner than the NIM, 3rdHAM, and 3rdFDM can. Specifically, the local errors ||u| − |u0|| of the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM were 0.0155, 0.0084, 0.0233, and 0.6272 × 10−3, respectively, when Ω = 5000 rpm, w = 0.5 mm, and m = 40. The local error of the proposed NI-HI-3rdFDM decreased by 90%, 92.5%, and 97.3% in comparison to that of the NIM, 3rdHAM, and 3rdFDM, respectively. Although the proposed NI-HI-3rdFDM converged a little slower than 3rdHAM in case of w = 0.5 mm and |u0| = 1.0726, for other three cases, the proposed NI-HI-3rdFDM was fairly faster than 3rdHAM. This slow convergence rate can be explained by the fact the derivability of the state item was utilized to enhance the calculation precision in 3rdHAM.Figure 1
Convergence rate for the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM. (a)w = 0.2 mm, u0=0.8192; (b) w = 0.5 mm, u0=1.0726; (c) w = 0.7 mm, u0=1.2197; (d) w = 1 mm, u0=1.4040.
(a)(b)(c)(d)


## 3.2. Stability Lobe Diagrams

In order to demonstrate the computation accuracy and computation time performances of the proposed NI-HI-3rdFDM, the SLDs obtained with the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM form = 20, 30, and 40 are summarized in Table 3. Note that the machining operation is up-milling operation, the ratio of radial cutting depth and tool diameter is a/D = 1, the spindle speed Ω, and depth of cut w is set to be 5000 rpm to 10000 rpm and 0 mm to 4 mm over a 200×100 sized grid of parameters, respectively [14]. It is also worth noting that in Table 3 the stability charts obtained with the 1stSDM for m = 200 in red color were regarded as the criterion-referenced SLDs. As can be seen, the stability charts of the proposed NI-HI-3rdFDM were much closer to the criterion-referenced lobe diagrams than those of the NIM, 3rdHAM and 3rdFDM. Table 3 indicates that the proposed NI-HI-3rdFDM can also achieve better or at least comparable performance with existing methods such as the NIM, 3rdHAM, and 3rdFDM with respect to the computation time. Specifically, in the case of m = 30, the computation time of the NIM, 3rdHAM, 3rdFDM, and the proposed NI-HI-3rdFDM were 10 s, 92 s, 38 s, and 23 s, respectively. In this case, the computation efficiency of the proposed NI-HI-3rdFDM increased by 75% and 39.5% compared with that of the 3rdHAM and 3rdFDM, respectively, except for the NIM. Although the computation time of the proposed NI-HI-3rdFDM was slightly greater than that of the NIM, the computation accuracy obtained with the proposed NI-HI-3rdFDM (corresponding computation time = 23 s) was even equal to that of the NIM when m = 60 (corresponding computation time = 37 s). Moreover, in the case of m = 40, the stability charts of the proposed NI-HI-3rdFDM were in good agreement with the criterion-referenced lobe. A visual examination of Table 3 also revealed that the proposed NI-HI-3rdFDM not only obtained successfully accurate milling stability prediction but also saved a lot of computation time.Table 3
SLD for the single DOF milling model using the NIM, 3rdHAM, 3rd FDM, and the proposed NI-HI-3rdFDM.
m = 20Computation time 6 sComputation time 55 sComputation time 22 sComputation time 15 sm = 30Computation time 10 sComputation time 92 sComputation time 38 sComputation time 23 sm = 40Computation time 21 sComputation time 125 sComputation time 54 sComputation time 37 sIn addition, the proposed NI-HI-3rdFDM can also be extended to a multifreedom milling system. A stability prediction of a two-DOF milling system was employed to evaluate the accurate and efficient performances of NI-HI-3rdFDM in predicting milling stability. Note that the same system parameters taken from the literature [5] were used. The SLDs calculated from up-milling operations with m = 40 and a/D = 1, 0.5, 0.1, and 0.05 are shown in Table 4. As can be seen in Table 4, the SLD of the proposed NI-HI-3rdFDM was relatively identical to that of the 1stFDM and the 1stSDM, but the computation time of the proposed NI-HI-3rdFDM was much less than that of the 1stFDM and the 1stSDM under the same milling parameters. The proposed NI-HI-3rdFDM achieved better and comparable efficiency compared with the 1stFDM and the 1stSDM. Especially for a/D = 0.05, the runtime of the proposed NI-HI-3rdFDM decreased from 200 s to 160 s when compared with the 1stFDM and from 1197 s to 160 s when compared with the 1stSDM.Table 4
SLD for the two-DOF milling model using the 1stSDM, 1stFDM, and the proposed NI-HI-3rdFDM.
Radial depth of cut ratio1stSDM [11]1stFDM [13]The proposed NI-HI-3rdFDMa/D = 1Computation time 1247 sComputation time 282 sComputation time 256 sa/D = 0.5Computation time 1189 sComputation time 243 sComputation time 190 sa/D = 0.1Computation time 1204 sComputation time 202 sComputation time 163 sa/D = 0.05Computation time 1197 sComputation time 200 sComputation time 160 s


## 4. Comparison and Discussion

Performances of the proposed NI-HI-3rdFDM were compared with those of existing methods in the literature. In order to illustrate the fact that specific interpolations of state and time-delay items can be associated independently with different problems when relevant process knowledge is accessible, the other candidate method called 3rdH-HAM for chatter stability prediction that uses two three-order Hermite approximation algorithms to simultaneously approximate the state item and time-delay item, respectively, is also provided in the comparison between the proposed NI-HI-3rdFDM and other existing methods. Figure2 compares the convergence rate for chatter stability prediction of the proposed NI-HI-3rdFDM with that of the 3rdH-HAM, 3rdUFDM, and 3rdH-NAM. As seen in Figure 2, the proposed NI-HI-3rdFDM outperformed all the compared methods in converging local discretization differences between the approximate and the exact solution over the discretization interval, except in the case of w = 0.2 mm and |u0| = 0.8192. In this case, there existed many weak points for the proposed NI-HI-3rdFDM when compared with the 3rdH-NAM. This relatively slow convergence can be explained by the fact that the derivative of xδ in the 3rdH-NAM was utilized to approximate the state term with low axial depth of cut, which resulted in the increase of accumulative precision. It can also be concluded from Figure 2, especially from this case, that the method using third-order Newton interpolation algorithm and third-order Hermite interpolation algorithm simultaneously approximated the state item and time-delay item, respectively, (i.e., NI-HI-3rdFDM) and outperformed the method using two three-order Hermite interpolation algorithms to simultaneously approximate the state item and time-delay item, respectively (i.e., 3rdH-HAM). This means that three-order Hermite interpolation algorithm was ideal only when it was a time-delay item interpolator and failed to be ideal for the interpolation of state item.Figure 2
Comparison of convergence rates between the proposed NI-HI-3rdFDM and existing methods. (a)w = 0.2 mm, u0=0.8192; (b) w = 0.5 mm, u0=1.0726; (c) w = 0.7 mm, u0=1.2197; (d) w = 1 mm, u0=1.4040.
(a)(b)(c)(d)Table5 presents a comparison of the SLDs of the 3rdH-HAM, 3rdH-NAM, 3rdUFDM, and the proposed NI-HI-3rdFDM in creating lobe diagrams. Table 5 indicates that the proposed NI-HI-3rdFDM can create lobe diagrams more accurately than the 3rdH-HAM, 3rdH-NAM, and 3rdUFDM can. The proposed NI-HI-3rdFDM was also capable of creating lobe diagrams with smaller value of m, but without any loss in computation accuracy; the 3rdH-HAM, 3rdH-NAM, and 3rdUFDM had no such capability. As can be seen, using the third-order Hermite interpolation algorithm to approximate the time-delay term enhanced notably performances of the proposed NI-HI-3rdFDM, especially the computation accuracy. For a certain accuracy of SLD, the third-order Hermite interpolation algorithm was more effective than the third-order Newton interpolation algorithm in approximating the time-delay term. The abovementioned characterization of ideal interpolation implies that when special interpolation algorithms existed approximating the state item and time-delay item simultaneously in creating SLD for chatter stability prediction, it was more reasonable to combine the state-item interpolation and time-delay-item interpolation information into one scheme and analyzed their behavior jointly. The simultaneous analysis of state and time-delay items was very helpful in the approximation search for each discretization step. The third-order Newton interpolation algorithm was obviously not good enough for the proposed NI-HI-3rdFDM to find the optimal interpolation solution of the state item. However, when the third-order Newton interpolation algorithm approximated the state item, the proposed NI-HI-3rdFDM increased the efficiency by 84.8% (1-37/243) compared with the 3rdH-NAM and by 17.8% (1-37/45) compared with the 3rdUFDM under the similar accuracy condition. Therefore, it can be concluded that a third-order Newton interpolation of state item and a third-order Hermite interpolation of time-delay item were a good compromise between the performance and the execution time of the proposed NI-HI-3rdFDM.Table 5
Comparison of SLDs between the proposed NI-HI-3rdFDM and existing methods.
Computation time 129 s,m = 40Computation time 243 s,m = 62Computation time 45 s,m = 57Computation time 37 s,m = 40


## 5. Experimental Verification

In this section, a two-DOF milling system withm = 40 was utilized to verify how the proposed NI-HI-3rdFDM can function effectively in the chatter stability prediction. The same experimental verification was also carried out by Gradišek et al. [29]. The modal parameters and cutting experiment data applied to the proposed NI-HI-3rdFDM were the same as those used in the study of Gradišek et al. [29], thereby ensuring unbiased verification. Figure 3 shows the SLDs for up-milling operations with a/D = 1, 0.5, 0.1, and 0.05. The results of the investigation indicate that for a/D = 1 and 0.5 the chatter stability boundaries of the proposed NI-HI-3rdFDM absolutely agreed with Gradišek’s experimental observations. For a/D = 0.1 and 0.05, the proposed NI-HI-3rdFDM also possessed good stability prediction capability, but not without weak points for some specific cutting regimes. It was relatively difficult for the proposed NI-HI-3rdFDM to predict the stability boundaries as radial immersion was decreased. The abovementioned weak points of the chatter stability prediction capability could be tackled by means of additional consideration of certain nonlinear factors such as tool deflection and effects of edge forces residing in the low radial immersion milling. The variation of the radial immersion during cutting due to the tool deflection and the effects of edge forces could be attached to the milling dynamics modeling to improve the chatter stability prediction performance of the proposed NI-HI-3rdFDM with regard to a/D = 0.1 and 0.05.Figure 3
Experimental and predicted stability boundaries.
(a)(b)(c)(d)


## 6. Conclusions

Regenerative chatter can be related to different assignable causes that affect adversely a machining operation, such as limiting metal removal rate below the machine’s capacity and reducing the productivity of the machine. Therefore, the prediction of regenerative chatter stability is very helpful for optimally stable cutting parameter selection in ensuing a high-performance machining. In past years, various full-discretization methods have been designed for predicting regenerative chatter stability. The main problem of such methods is that they can predict the regenerative chatter stability but do not efficiently determine SLDs. In this study, using third-order Newton interpolation and third-order Hermite interpolation techniques, a straightforward and effective NI-HI-3rdFDM was developed for the prediction of regenerative chatter stability. The experimental results obtained using simulation indicate that the proposed NI-HI-3rdFDM can not only efficiently predict the regenerative chatter stability but also accurately identify the SLD. Empirical comparisons showed that the proposed NI-HI-3rdFDM outperformed the existing approaches in the literature in predicting the regenerative chatter stability, while also providing the capability of faster and more accurate computation that facilitates dynamic milling scenario, in which a computationally efficient and accurate chatter stability prediction method is required. This study also demonstrated that the chatter stability boundaries of the proposed NI-HI-3rdFDM agree with experimental observations.In this study, only a linear milling dynamics model was adopted; thus, the predictions and experiments agree only qualitatively (i.e., with respect to the structure of the stability boundary) for lower radial immersions. Hence, for applying the proposed NI-HI-3rdFDM in a practical chatter stability prediction scenario, other common types of nonlinear factors (e.g., tool deflection and effects of edge forces) should also be included in the milling dynamics modeling to overcome the usage limitation of the proposed method.



---

*Source: 9071451-2020-06-20.xml*