An updated Simpson-Based Method for Chatter Stability Prediction in Milling

: An updated Simpson - based method (USBM) is presented for milling stability analysis. Firstly, the delay differential equation (DDE) is employed to describe the milling process mathematically. Then, the tooth passing period is divided into two subintervals, i.e., the free and forced vibration intervals. Only the forced vibration interval is divided into many equal small - time intervals. Subsequently, the DDE in the state space is solved based on direct integration. By combining the two - step Simpson method and the semi - discretization method, the state transition matrix of the milling system is constructed. The comparison of convergence rate is conducted to validate the accuracy of the proposed method. The results show that the proposed method converges faster than the benchmark methods. The stability lobe diagrams for the one degree of freedom (one - DOF) and two degrees of freedom (two - DOF) milling systems are also obtained by different methods for further evaluation. Meanwhile, the computation time analysis is also carried out. It is revealed that the proposed USBM has advantages in both accuracy and efficiency. Besides, the proposed method can accurately and efficiently predict the stability of milling with both large and low immersion conditions.


Introduction
Milling operation is an effective processing technology, which is suitable for the machining of workpieces with various complex shapes. Regenerative chatter, which harms surface quality and productivity, is easy to initiate in the milling process. Moreover, chatter increases the cutting forces and may chip the tool and produce a poor, wavy surface finish [1]. Therefore, avoiding chatter is crucial to ensure the high surface finish with high material removal rate. The stability lobe diagram, which divides the parameter plane of spindle speed and axial depth of cut into stable and unstable regions, can be used to select chatter-free parameters. The process parameter combinations in the stable region can be used to avoid chatter and achieve high productivity. Therefore, the demand for developing accurate and efficient milling stability prediction methods is increasing.
The delay differential equation (DDE) is usually used to describe the milling process mathematically. By solving the DDE, the milling stability prediction methods can be derived. Up to now, extensive stability prediction methods in milling have been proposed. Altintas and Budak [2] presented a well-known and efficient zeroth-order approximation (ZOA) method. In this method, the DDE is transformed into the frequency domain by Laplace transform. Science the axial depth of cut is expressed analytically, the computational efficiency of the ZOA method is high. Considering that only a single harmonic of the tooth passing period is considered during the approximation of directional coefficients, the prediction accuracy of the ZOA method is not high at low immersion milling. After that, Merdol and Altintas [3] extended the ZOA method and proposed a multi-frequency method. Since high-order harmonics of the tooth passing period are used to approximate the directional coefficients, the multi-frequency method can accurately predict the stability at low immersion milling. In addition, the ZOA method is also extended for the analysis of stability in milling with non-uniform cutters [4][5][6].
The ZOA method and multi-frequency method are frequency-domain methods. Many time-domain methods have also been reported. Bayly et al. [7] presented a temporal finite element analysis method, and Butcher et al. [8] proposed a Chebyshev collocation method. Insperger et al. presented the semi-discretization method (SDM) [9], the updated SDM [10], and the first-order SDM (1st SDM) [11]. The SDMs proposed by Insperger et al. are time-consuming methods. To improve the computational efficiency, Ding et al. [12] developed a full-discretization method (FDM), which is proved to be more efficient than the SDM. Then, the second-order FDM (2nd FDM) [13], the third-order FDM (3rd FDM) [14], and the hyper-third-order FDMs [15] were proposed one after another. The comparative study shows that the prediction accuracy reaches the highest when the order of the full-discretization method is four. In the framework of FDM, extensive methods [16][17][18][19][20][21][22][23][24][25][26][27][28] have been proposed to explore the method with both high prediction accuracy and high computational efficiency. In recent years, Jiang et al. [29] developed a second-order SDM (2nd SDM), in which the delayed term is approximated by second-order Newton interpolation polynomial. In this method, the precise time-integration (PTI) algorithm is used to improve computational efficiency. Liu et al. [30] suggested an improved semi-discretization method based on the predictor-corrector scheme (PCSDM) for milling stability analysis. Besides, many other time-domain methods have also been suggested to predict the stability of milling. Ding et al. [31] reported the numerical integration method (NIM) based on the integral equation and numerical integration formulas. Dong et al. [32] developed an updated numerical integration method (UNIM) based on the NIM. Li et al. [33] presented a complete discretization scheme (CDS). In the CSD, the differential term of the DDE is discretized using Euler's method. Inspired by the CSD, Xie et al. [34] proposed an improved complete discretization method. Niu et al. [35] recommended two Runge-Kutta-based methods, namely the classical fourth-order Runge-Kutta method (CRKM) and the generalized form of the Runge-Kutta method (GRKM). The GRKM is proved to converge faster than the CRKM. Li et al. [36] suggested a Runge-Kutta-based complete discretization method (RKCDM) by combining the CRKM and the CDS. Dai et al. [37] presented a milling stability prediction method using the precise integration method. By approximating the differential term of the DDE, several methods were proposed for milling stability analysis, such as the differential quadrature method [38], the numerical differentiation method [39], and the Chebyshev-wavelet-based method [40]. Based on the numerical solution of ordinary differential equations, Qin et al. presented the Adams-Moulton-based method (AMM) [41] and the Adams-Simpson-based method (ASM) [42]. The ASM is proved to be more accurate than the AMM.
The Simpson-based method (SBM) suggested by Zhang et al. [43] achieves satisfactory accuracy and efficiency. In this work, we attempt to propose an updated Simpson-based method by combining the two-step Simpson method and the SDM to improve accuracy and efficiency further. The framework of this paper is as follows: Section 2 illustrates the Mathematical model of milling. Section 3 shows the derivation process of the proposed method. Section 4 gives the comparison and discussion. In this section, the proposed method is validated from the aspect of convergence rate and computational time. The main conclusions are given in Section 5.

Mathematical model of milling
The milling process considering regenerative effect can be described by the following equation where M, C, and K denote the modal mass, damping, and stuffiness matrices, ap represents the axial depth of cut, q(t) represents the displacement vector, and Kc(t) is the directional cutting force coefficients matrix. T is the time delay which equals the tooth passing period, i.e., T=60/(NΩ), where N is the number of cutter teeth, and Ω is the spindle speed in rpm. and For a one-DOF milling system only considering a single dominant structural mode in the X direction, the matrices A0 and B(t) can be given as For a two-DOF milling system considering the single dominant structural mode in two orthogonal directions (X and Y), the matrices A0 and B(t) can be given as where mx, ζx, and ωnx are the modal mass, damping ratio, and angular natural frequency in the X direction, respectively, and my, ζy, and ωny are the modal mass, damping ratio, and angular natural frequency in the Y direction, respectively. The directional cutting force coefficients hxx(t), hxy(t), hyx(t), and hyy(t) are given as follows: where Ktc and Krc are the tangential and the normal cutting force coefficients, respectively. g[φj(t)] is a window function which determines whether the tooth is in or out of the cut, where φj(t) is the angular position of the tooth j, and it is given as follows:

The proposed method
In this study, the tooth passing period is divided into two subintervals, i.e., the free and forced vibration intervals. After that, the forced vibration interval Tfo is divided into n equal small-time intervals with the step length of h. Then, the time point ti can be represented as follows: where t0 is the initial time point, and Tfr is the free vibration interval. Science the milling system experiences free vibration from the time point tn+1-T to the time point t1, the U(t1) and U(tn+1-T) denoted as U1 and Un+1-T have the following relation Equation (2) is solved on the interval [ti-1, ti+1], we can get where Ui+1 and Ui are the abbreviations of the U(ti+1) and U(ti), respectively.
Then, based on two-step Simpson method, the following result can be obtained where Ui-1-T, Ui-T, and Ui+1-T are the abbreviations of the U(ti-1-T), U(ti-T) and U(ti+1-T), respectively. Equation (19) is inserted into Eq. (18), we can get where When i=1, the trapezoidal rule is used to solve the Eq. (2) on the interval [ti, ti+1], we can get According to Eqs. (16), (20), and (25), the following mapping relation can be obtained as 1 1 The state transition matrix Ψ of the milling system can be obtained as Ψ=(D1) -1 D2 (31) Based on Floquet theory, the milling stability can be predicted according to the spectral radius of the state transition matrix [10].
To improve the computational efficiency of the proposed method, the precise integration algorithm [44] is used for calculating the matrix exponent i h e A , and the precision exponent is chosen as 2, which is the same as that adopted in Ref. [29].

Convergence rate
The convergence rate can be taken as the criterion to evaluate the accuracy of the stability prediction methods. Normally, the method with higher convergence rate is relatively more accurate. In Refs. , the local discretization error calculated by ||μ(n)|-|μ0|| is employed for the analysis of convergence rate, where |μ(n)| and |μ0| are the approximated and exact spectral radius of the state transition matrix. The convergence rate reflects how fast the local discretization error approaches zero. When the discrete number n is large enough, the local discretization error calculated by the prediction methods will approach zero. In this work, |μ0| is determined using 1st SDM with n=600. To evaluate the convergence rate of the proposed method, the 1st SDM, 2nd SDM, SBM, and PCSDM are taken as the benchmark. The spindle speed Ω is set as 5000 rpm, and the axial depth of cuts ap is set as 0.   Fig. 1(a) and (b), when ap=0.2 and 0.5 mm, the local discretization error calculated by the USBM approaches zero when n equals 40, while that calculated by the benchmark methods tends to be zero when n is much larger than 40. In Fig. 1(c), when ap=1 and n=30, the local discretization error obtained by the USBM almost equals zero, while that obtained by the benchmark methods is larger than 0.02. 30

One-DOF milling
In section 4.1, the accuracy of the proposed method is verified from the aspect of convergence rate. In this section, the accuracy and efficiency of the proposed method are further validated from the aspect of the stability lobe diagram (SLD). The proposed method is compared with the 1st SDM, SBM, 2nd SDM, and PCSDM in terms of the SLD. Firstly, the one-DOF milling system is used for validation. The As shown in Fig. 2, when n takes the same value (30 and 40), the SLDs obtained by the USBM and PCSDM are closer to the reference than those obtained by the 1st SDM, SBM, and 2nd SDM. Regarding the computational time, the SBM and USBM consume less time than the 1st SDM, 2nd SDM, and PCSDM to generate SLDs. Although the USBM takes slightly more time than the SBM, the result predicted by the USBM is more accurate than that predicted by the SBM. The time increment between the USBM and SBM is small. Through contrastive analysis, it is revealed that the proposed USBM has advantages in both prediction accuracy and computational efficiency.
The stability of the one-DOF milling system under the low immersion condition (a/D=0.05) is also predicted. The SLD is also plotted on a 200×200 sized equidistance grid with Ω[5000, 10000] rpm and ap[0, 0.01] m. The comparison of the SLDs calculated by the 1st SDM, SBM, 2nd SDM, PCSDM, and USBM with n=20 and 30 for the one-DOF milling system is presented in Fig. 3. As shown in Fig. 3, PCSDM has advantages in prediction accuracy but not in computation efficiency. The existing SBM is an efficient method, but whose prediction accuracy needs to be improved. Unlike the benchmark methods, the proposed USBM can achieve both high accuracy and high efficiency.

Two-DOF milling
To further evaluate the effectiveness of the USBM, the two-DOF milling system is also considered in this study. The cutting force coefficients are the same as those used in section 4.1. The modal parameters in the X and Y directions are assumed to be equal and listed as follows: ωnx=ωny=5793 rad/s, ζx=ζy=0.011, and mx=my=0.03993 kg. The SLD is plotted on a 200×200 sized equidistance grid with Ω[5000, 25000] rpm and ap[0, 0.001] m. The SLD determined by the SBM with n=300 is regarded as the reference. The SLDs obtained by different methods with n=20 and 30 for the two-DOF milling system under the large immersion condition (a/D=1.0) are presented in Fig. 4. As shown in Fig. 4, the SBM takes the least time to generate SLDs; however, the SLDs obtained by the 1st SDM, 2nd SDM, PCSDM, and USBM are closer to the reference than those obtained by the SBM. In Fig. 4, when n=20 and 30, the increment of computational time between the SBM and USBM is one second and two seconds, respectively, which indicates that the efficiency of the SBM and USBM is almost the same. The comparison shows that the USBM can predict the milling stability for the two-DOF milling system under the large immersion condition accurately and efficiently.
The SLDs obtained by the 1st SDM, SBM, 2nd SDM, PCSDM, and USBM with n=20 and 30 for the 2-DOF milling system under low immersion condition (a/D=0.05) are presented in Fig. 5. As shown in Fig. 5, the prediction accuracy of the PCSDM and the proposed USBM is higher than that of the 1st SDM, 2nd SDM, and SBM. Besides, the efficiency of the proposed USBM is higher than the other methods.
In general, the proposed USBM has good accuracy and efficiency for both one-DOF and two-DOF milling systems and large and low radical immersion conditions.

Conclusions
This study proposed an updated Simpson-based method for stability prediction in milling. In this study, the two-step Simpson method is employed to solve the DDE in the framework of SDM. The main conclusions sum up as follows: (1) The proposed USBM converges faster than the benchmark methods. When the discrete number n is identical, the local discretization error obtained by the USBM is less than that obtained by the benchmark methods.
(2) For the one-DOF milling system under the large immersion condition (a/D=1.0), the SLDs obtained by the USBM and PCSDM are closer to the reference than those obtained by the 1st SDM, SBM, and 2nd SDM. Meanwhile, the SBM and USBM consume less time than the 1st SDM, 2nd SDM, and PCSDM to generate SLDs. Generally, the proposed USBM has advantages in both accuracy and efficiency.
(3) Under the low immersion condition (a/D=0.05), the proposed USBM is superior in prediction accuracy and computation efficiency for the one-DOF milling system.
(4) As for the two-DOF milling system, the prediction accuracy of the PCSDM and the proposed USBM is higher than that of the 1st SDM, 2nd SDM, and SBM. Besides, the efficiency of the proposed USBM is higher than the benchmark methods. Competing interests The authors declare that there is no conflict of interests regarding the publication of this article. Availability of data and material All data generated or analyzed during this study are included in this published article. Code availability Not applicable. Authors' contributions Zhenghu Yan and Changfu Zhang and contributed the central idea, data analysis, and manuscript writing. Jianli Jia, Baoji Ma, and Xinguang Jiang contributed to carrying out additional analyses. Dong Wang, Wei Wang, and Chenxi Yang finalizing this paper. All authors read and approved the final manuscript. Ethical Approval Not applicable. Consent to Participate Not applicable. Consent to Publish Not applicable.