Multi-modal method for chatter stability prediction and control in milling of thin-walled workpiece

Journal Pre-proof

Multi-modal method for chatter stability prediction and control in milling of thin-walled workpiece Yichao Dun , Lida Zhu , Shuhao Wang PII: DOI: Reference:

S0307-904X(19)30742-5 https://doi.org/10.1016/j.apm.2019.12.003 APM 13189

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

12 December 2018 28 November 2019 3 December 2019

Please cite this article as: Yichao Dun , Lida Zhu , Shuhao Wang , Multi-modal method for chatter stability prediction and control in milling of thin-walled workpiece, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.12.003

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights 

A numerical difference method for predicting stability in milling.



Time- and position- dependent change of workpiece modal parameters, and its effect on stability prediction.



A multi-modal scheme of numerical methods.



Comparison between analytical methods and numerical methods.



Time complexity of analytical methods and numerical methods.

1

Multi-modal method for chatter stability prediction and control in milling of thin-walled workpiece Yichao Dun1, Lida Zhu1*, ShuhaoWang School of Mechanical Engineering and Automation, Northeastern University, 110819, China 1Co-first

author: these authors contributed equally to this work.

* Corresponding author: [email protected] (L.D. Zhu) Abstract Chatter stability in milling can be predicted by analytical methods or numerical methods. The system should be considered as multi-modal in milling of thin-walled workpiece. This paper proposes a numerical difference method based on Adams-Bashforth scheme. Moreover, multi-modal scheme of numerical methods is proposed. Analytical methods and numerical methods are verified by performing a series of milling trials. The experimental results are consistent with the predicted critical stability boundaries. Moreover, a new method for analyzing the computational time of analytical methods and numerical methods, which is based on the time complexity modeling, is presented. Computational time can be expressed as exact mathematical expression. By using the expression, the rate of increase of computational time can be derived. Keywords: milling; thin-walled workpiece; stability analysis in frequency domain; stability analysis in time domain; computational time 1. Introduction End milling is used in various applications, such as aerospace, automobile and others[1]. Cutting forces occurring in machining processes are the main causes that may lead to regenerative chatter[2]. Chatter is a form of self-excited vibration, which can lead to poor surface finish, short cutter lifetime, even damage the spindle[3]. Mathematical modeling of chatter analysis in metal cutting, such as turning or milling, can be expressed as a delay differential equation[4]. Continuous machining, such as turning, may cause the Hopf bifurcation[5]. Interrupted machining, such as milling, may cause Hopf bifurcation, period doubling bifurcation or flip bifurcation[6]. Prediction of chatter stability in milling is usually shown in a stability lobes diagram. In a stability lobes diagram, the abscissa represents spindle speed and the ordinate represents critical stability cutting depth[6]. The critical stability cutting depth can be solved analytically. Altintas and Budak[7] proposed a zero order approximation method (ZOA), which uses the frequency domain method to analytically solve the critical stability cutting depth and used the arithmetic mean to approximate directional coefficients. Then, Merdol and Altintas[8] proposed a multi-frequency solution method, which uses the Fourier series to approximate directional coefficients. Iglesias et al.[9] proposed a new analytical formulae related to the parameter domains of both Hopf and period doubling type stability boundaries emerging in the regenerative mechanical model of time periodical milling processes. Zhu et al. [10] considered process damping in the analytical method. Graham et al.[11] propose a robust chatter stability model based on the ZOA. The critical stability cutting depth can also be solved numerically. Insperger and Stepan[12] proposed a semi-discretization method, which numerically solves the critical stability cutting depth by approximating the direction coefficients using the arithmetic mean and using linear interpolation to approximate the time-delay term. To simplify the algorithm, the time-delay term can also be approximated by the arithmetic mean[13]. Ding et al.[14] proposed a full-discretization method, which numerically solves the critical stability cutting depth by using linear interpolation to approximate both direction coefficients and time-delay term. Ding et al.[15] solved the critical stability cutting depth by using numerical integration method. Compean et al.[16] analyzed the chatter stability of a variable pitch, variable helix and variable rake angle cutter by using Enhanced Multistage Homotopy Perturbation Method. Urbikain et al.[17] predicted critical stability boundaries by using Chebyshev collocation method. Then, Urbikain et al.[18] analyzed the multimode chatter stability using single frequency model and Chebyshev polynomials. Yue et al.[19] summarized the current state of the art in research regarding the problems of how to

2

arrive at stable chatter prediction, chatter identification, and chatter control/suppression, with a focus on milling processes. Bayly et al.[20] discussed the stability predictions at small radial immersions. Insperger[21] compared the method in reference [12] with the method in reference [14]. Ahmadi and Ismail[22] integrated the equivalent viscous model of process damping into multi-frequency solution method and semi-discretization method. The dynamic characteristics of thin-walled workpieces should be considered when milling thin-walled workpieces. Denkena and Schmidt[23] presented an enhanced simulation model to predict form deviations in end milling processes of thin-walled structures. Biermann et al.[24] computed regenerative workpiece vibrations during the five-axis milling process by using FE model. Bravo et al.[25] took the flexibility of the workpiece into account when analyzing the stability in milling. Budak et al.[26] predicted the dynamic characteristics of the workpiece by using the FE model of the workpiece. Thevenot et al.[27] generated stability lobes by considering the dynamical behavior variation of workpiece. Zhu et al.[28][29] predicted cutting force acting on thin-walled workpiece. Arnaud et al.[30] analyzed the variation of thin-walled workpiece dynamic characteristics due to metal removal and the cutting edge position by using finite element method. Iglesias et al.[31] propose a new method for frequency response function estimation using the real milling force as the input excitation, which allows obtaining the FRF under real cutting conditions. Olvera et al.[32] studied how the mass and the lightweight connecting cable of an ultra-miniature accelerometer sensor influence the frequency response function of thin-walled aluminum workpieces. Zulaika et al. [33] provided an integrated approach for designing large milling machines，which is based on modeling the interactions between process and machine by means of a stability model of the milling process in modal coordinates. Munoa et al.[34]-[36] actively suppressed chatter in milling. Kersting and Biermann[37] proved that the modal parameters of the thin-walled workpiece will change as the tool-workpiece contact position changes. Thevenot et al. [38] proved that the modal parameters of thin-walled workpieces change with material removal. This paper presents a comparison between analytical methods and numerical methods on predicting stability in milling of thin-walled workpiece considering position- and time-dependent feature of modal parameters of thin-walled workpiece. This paper proposes a numerical difference method based on Adams-Bashforth scheme and multi-modal scheme of numerical methods. And the critical stability boundaries predicted by analytical methods and numerical methods are verified by milling experiments. Furthermore, time complexity modeling of analytical methods and numerical methods are presented and then verified by measuring elapsed time on computers. Parallel algorithms of analytical methods and numerical methods which can reach linear speedup are also presented in this paper. The robustness of the proposed method was also discussed. Section 2 introduces the mathematical models about chatter stability in milling of thin-walled workpiece and figures out time complexity of different methods. Section 3 measures modal parameters of the system then predicted stability border of the system by analytical methods and numerical methods then analyzes difference between them. Section 4 verifies the accuracy of the simulations through the milling experiment. 2. Model of chatter stability in milling of thin-walled workpiece 2.1 Dynamics modeling in milling The thin-walled workpiece end milling process is shown in Fig. 1. The cutter is simplified as a spring damper structure in the X direction and the Y direction[7]. The thin-walled workpiece is simplified as a spring damper structure in the Y direction[27].

3

Fig. 1. Model of thin-walled workpiece milling operation In Fig. 1, the thickness of the chip at the No. j tool tip can be expressed as Eq. (1). 2𝜋 2𝜋 𝐶𝑕𝑖𝑝,𝑗 (𝑡) = sin ( 𝑡* [𝑥(𝑡) − 𝑥(𝑡 − 𝑇)] + cos ( 𝑡* [𝑦(𝑡) − 𝑦(𝑡 − 𝑇)] 𝑁𝑇 𝑁𝑇

(1)

Then, the tangential cutting force can be expressed as Eq. (2). 𝐹𝑡,𝑗 = 𝐾𝑡 𝑎𝑝 𝐶𝑕𝑖𝑝,𝑗 (𝑡)

(2)

The radial cutting force can be expressed as Eq. (3). 𝐹𝑟,𝑗 = 𝐾𝑟 𝐾𝑡 𝑎𝑝 𝐶𝑕𝑖𝑝,𝑗 (𝑡)

(3)

Then, project the tangential cutting force and the radial cutting force into the X direction and the Y direction. 2𝜋 2𝜋 𝐹𝑥,𝑗 = −𝐹𝑡,𝑗 cos ( 𝑡* − 𝐹𝑟,𝑗 sin ( 𝑡* 𝑁𝑇 𝑁𝑇 𝐹𝑦,𝑗

(4)

2𝜋 2𝜋 = 𝐹𝑡,𝑗 sin ( 𝑡* − 𝐹𝑟,𝑗 cos ( 𝑡* 𝑁𝑇 𝑁𝑇

Then, the resultant force on all tool tip can be expressed as Eq. (5). 𝑁

2𝜋 2𝜋 𝐹𝑥 = ∑ [𝑔 ( 𝑡 + 𝑗 * 𝐹𝑥,𝑗 ] 𝑁𝑇 𝑁 𝑗=1 𝑁

2𝜋 2𝜋 𝐹𝑦 = ∑ [𝑔 ( 𝑡 + 𝑗 * 𝐹𝑦,𝑗 ] 𝑁𝑇 𝑁

(5)

𝑗=1

1 , 𝜑𝑠𝑡 < 𝜑 𝑚𝑜𝑑 2𝜋 < 𝜑𝑒𝑥 𝑔(φ) = { 0 𝑒𝑙𝑠𝑒 The dynamics model in milling can be expressed as a delay differential equation in Eq. (6). (6)

𝐌 ∙ 𝐬̈ (𝑡) + 𝐂 ∙ 𝐬̇ (𝑡) + 𝐊 ∙ 𝐬(𝑡) = 𝐅(𝐬(𝑡), 𝐬(𝑡 − 𝑇), 𝑡) Where: 𝑚𝑥 𝐌=( 0

0 𝑐𝑥 𝑚𝑦 *, 𝐂 = ( 0

𝑎𝑥𝑥 (𝑡) 𝐔(𝑡) = ( 𝑎𝑦𝑥 (𝑡)

0 *, 𝑘𝑦

𝑥(𝑡) 𝐬(𝑡) = ( *, 𝑦(𝑡)

𝑥(𝑡) − 𝑥(𝑡 − 𝑇) 1 𝐅(𝐬(𝑡), 𝐬(𝑡 − 𝑇), 𝑡) = 𝑎𝑝 𝐾𝑡 𝐔(𝑡) ( *, 2 𝑦(𝑡) − 𝑦(𝑡 − 𝑇)

𝑎𝑥𝑦 (𝑡) 2𝜋 2𝜋 2𝜋 2𝜋 ), 𝑎𝑥𝑥 (𝑡) = −2 cos (𝑁𝑇 𝑡) sin (𝑁𝑇 𝑡) − 2 sin (𝑁𝑇 𝑡) sin (𝑁𝑇 𝑡) 𝐾𝑟 , 𝑎𝑦𝑦 (𝑡)

𝑎𝑥𝑦 (𝑡) = −2 cos ( 𝑎𝑦𝑦 (𝑡) = 2 sin (

0 𝑘𝑥 𝑐𝑦 *, 𝐊 = ( 0

2𝜋

𝑁𝑇 2𝜋

𝑁𝑇

𝑡) cos (

𝑡) cos (

2𝜋

𝑁𝑇 2𝜋

𝑁𝑇

𝑡) − 2 sin (

𝑡) − 2 cos (

2𝜋

𝑁𝑇 2𝜋

𝑁𝑇

𝑡) cos (

𝑡) cos (

2𝜋

𝑁𝑇 2𝜋

𝑁𝑇

𝑡) 𝐾𝑟 , 𝑎𝑦𝑥 (𝑡) = 2 sin (

2𝜋 𝑁𝑇

𝑡) sin (

(7) 2𝜋 𝑁𝑇

𝑡) − 2 cos (

2𝜋 𝑁𝑇

𝑡) sin (

2𝜋 𝑁𝑇

𝑡) 𝐾𝑟 ,

𝑡) 𝐾𝑟

Multi-frequency solution method [13] transforms Eq. (6) into the frequency domain and then expands the two sides of Eq. (6) into Fourier series. Then, Eq. (6) is transformed into an equation whose form is similar to the definition of eigenvalue. Therefore, the critical stability cutting depth can be solved by solving eigenvalues of a matrix. Finally, because cutting depth

4

physically means a length, the solution with zero imaginary part and positive real part is the valid critical stability cutting depth. Bravo et al.[25] proves that the system’s transfer function in Fig. 1 is equal to cutter’s transfer function plus workpiece’s transfer function. Eq. (6) is a second-order ordinary differential equation. By left-multiplying Eq. (6) by 𝐌 −1 then substituting 𝑢(𝑡) = 𝑥̇ (𝑡) and 𝑣(𝑡) = 𝑦̇ (𝑡) into Eq. (6), Eq. (6) can be transformed into a first order ODE as follows: (8)

𝐪̇ (𝑡) = 𝐑𝐪(𝑡) + 𝐀(𝑡)𝐪(𝑡) + 𝐁(𝑡)𝐪(𝑡 − 𝑇) Where: 𝑢̇ (𝑡) 𝑣̇ (𝑡) 𝐪̇ (𝑡) = ( ), 𝑥̇ (𝑡) 𝑦̇ (𝑡) 0 0 𝐁(𝑡) =

−2𝜁𝑥 𝜔𝑛𝑥 0 𝐑=( 1 0 1

− 𝑎𝑝 𝐾𝑡 × 2

𝑘𝑥

1

2 𝜔𝑛𝑦

2

𝑘𝑦

0 0 − 𝑎𝑝 𝐾𝑡 × 0 0 (0 0

2 𝜔𝑛𝑥

0 −2𝜁𝑦 𝜔𝑛𝑦 0 1

× 𝑎𝑥𝑥 (𝑡)

2 −𝜔𝑛𝑥 0 0 0

1

− 𝑎𝑝 𝐾𝑡 ×

2 𝜔𝑛𝑥

2

𝑘𝑥

1

2 𝜔𝑛𝑦

2

𝑘𝑦

× 𝑎𝑦𝑥 (𝑡) − 𝑎𝑝 𝐾𝑡 ×

0 0

0 2 −𝜔𝑛𝑦 ), 𝐀(𝑡) = 0 0

0 0

0

0

0

0

0 (0

0 0

1 2 1 2

𝑎𝑝 𝐾𝑡 × 𝑎𝑝 𝐾𝑡 ×

2 𝜔𝑛𝑥

𝑘𝑥 2 𝜔𝑛𝑦

𝑘𝑦

× 𝑎𝑥𝑥 (𝑡) × 𝑎𝑦𝑥 (𝑡)

1 2 1 2

𝑎𝑝 𝐾𝑡 × 𝑎𝑝 𝐾𝑡 ×

0 0

2 𝜔𝑛𝑥

𝑘𝑥 2 𝜔𝑛𝑦

𝑘𝑦

0 0

× 𝑎𝑥𝑦 (𝑡) × 𝑎𝑦𝑦 (𝑡) , )

(9)

× 𝑎𝑥𝑦 (𝑡) × 𝑎𝑦𝑦 (𝑡) )

Fig. 2. Intermittent cutting process during one tooth passing period The numerical methods are based on the Floquet Theory. Floquet Theory is a stability criterion. A system is stable only if the modulus of all of eigenvalues of Floquet Matrix are less than one. Therefore, the key content of the numerical methods is to construct Floquet Matrix. Semi-discretization method [12] discretizes one tooth passing period into several equal intervals, as shown in Fig. 2 (b). Then, during each interval, by approximating 𝐀(𝑡) and 𝐁(𝑡) using the means and approximating 𝐪(𝑡 − 𝑇) using linear interpolation, the semi-discretization method transforms Eq. (8) into a series of first order non-homogeneous linear differential equations. Finally, the Floquet Matrix can be obtained by multiplying the coefficient matrices of the differential equations in all intervals. Full-discretization method [14] is similar to semi-discretization method. Full-discretization method approximates 𝐀(𝑡), 𝐁(𝑡) and 𝐪(𝑡 − 𝑇) using linear interpolation. The numerical integration method [15] discretizes the cutting part within one tooth passing period into several equal intervals, i.e. the part A in Fig. 2 (c). Then, numerical integration is performed on Eq. (8) during each interval. Then, the Floquet Matrix consists of coefficients of the integration schemes during all intervals. 2.2 A higher-order piecewise numerical difference method based on Adams-Bashforth scheme Based on the discretization mode of Numerical Integration Method, this paper uses a third-order numerical difference method to construct the Floquet Matrix. Milling is a typical intermittent cutting process. This means that one complete tooth passing period is divided into two parts. In Part A, the tooth is not cutting the workpiece. In part B, the tooth is cutting the workpiece. As shown in Fig. 2 (a), the initial condition is 𝐪(𝑡𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ) = 𝐪𝑖𝑛𝑖𝑡𝑖𝑎𝑙 . First, in part A (from 𝑡𝑖𝑛𝑖𝑡𝑖𝑎𝑙 to 𝑡0 ), because 𝑎𝑥𝑥 (𝑡) = 𝑎𝑥𝑦 (𝑡) = 𝑎𝑦𝑥 (𝑡) = 𝑎𝑦𝑦 (𝑡) = 0, Eq. (8) becomes Eq. (10) when 𝑡 = 𝑡0 [15].

𝐪(𝑡0 ) = 𝑒 𝐑(𝑡0−𝑡𝑖𝑛𝑖𝑡𝑖𝑎𝑙 ) 𝐪(𝑡𝑖𝑛𝑖𝑡𝑖𝑎𝑙 )

(10)

5

Then, part B (from 𝑡0 to 𝑡𝑚 ) is discretized into m equal intervals, as shown in Fig. 2 (b). The span of each interval is 𝜏2 . In Part B, by using the third-order Adams-Bashforth scheme, 𝐪(𝑡𝑘 ) can be expressed as Eq. (11). 23𝜏2 16𝜏2 5𝜏2 𝐪̇ (𝑡𝑘−1 ) − 𝐪̇ (𝑡𝑘−2 ) + 𝐪̇ (𝑡𝑘−3 ) 12 12 12

𝐪(𝑡𝑘 ) = 𝐪(𝑡𝑘−1 ) +

(11)

Then, by substituting Eq. (8) into Eq. (11), Eq. (12) can be obtained. 23𝜏2 16𝜏2 (𝐑𝐪(𝑡𝑘−1 ) + 𝐀𝑘−1 𝐪(𝑡𝑘−1 ) + 𝐁k−1 𝐪(𝑡𝑘−1 − 𝑇)) − (𝐑𝐪(𝑡𝑘−2 ) + 𝐀𝑘−2 𝐪(𝑡𝑘−2 ) + 𝐁k−2 𝐪(𝑡𝑘−2 − 𝑇)) 12 12

𝐪(𝑡𝑘 ) = 𝐪(𝑡𝑘−1 ) +

(12)

5𝜏2 + (𝐑𝐪(𝑡𝑘−3 ) + 𝐀𝑘−3 𝐪(𝑡𝑘−3 ) + 𝐁k−3 𝐪(𝑡𝑘−3 − 𝑇)) 12 Eq. (12) can be expressed equivalently as Eq. (13). 𝚪𝑑𝑖𝑓3𝑎,𝑘 𝐪(𝑡𝑘 ) + 𝚪𝑑𝑖𝑓3𝑏,𝑘−1 𝐪(𝑡𝑘−1 ) + 𝚪𝑑𝑖𝑓3𝑐,𝑘−2 𝐪(𝑡𝑘−2 ) + 𝚪𝑑𝑖𝑓3𝑑,𝑘−3 𝐪(𝑡𝑘−3 )

(13)

= 𝚪𝑑𝑖𝑓3𝑒,𝑘−1 𝐪(𝑡𝑘−1 − 𝑇) + 𝚪𝑑𝑖𝑓3𝑓,𝑘−2 𝐪(𝑡𝑘−2 − 𝑇) + 𝚪𝑑𝑖𝑓3𝑔,𝑘−3 𝐪(𝑡𝑘−3 − 𝑇) Where: 𝚪𝑑𝑖𝑓3𝑎,𝑘 = 𝑰, 𝚪𝑑𝑖𝑓3𝑒,𝑘 =

𝚪𝑑𝑖𝑓3𝑏,𝑘 = −𝑰 −

23𝜏2 12

𝐁k ,

23𝜏2 12

𝚪𝑑𝑖𝑓3𝑓,𝑘 = −

𝐑− 16𝜏2

23𝜏2 12

𝐀𝑘 ,

𝐁k ,

12

𝚪𝑑𝑖𝑓3𝑐,𝑘 =

𝚪𝑑𝑖𝑓3𝑔,𝑘 =

5𝜏2 12

16𝜏2 12

𝐑+

16𝜏2 12

𝐀𝑘 ,

𝚪𝑑𝑖𝑓3𝑑,𝑘 = −

5𝜏2 12

𝐑−

5𝜏2 12

𝐀𝑘 ,

(14)

𝐁k

Similarly, by using the second-order Adams-Bashforth scheme, 𝐪(𝑡2 ) can be expressed as Eq. (15). 𝐪(𝑡2 ) = 𝐪(𝑡1 ) +

3𝜏2 𝜏2 𝐪̇ (𝑡1 ) − 𝐪̇ (𝑡0 ) 2 2

(15)

By substituting Eq. (8) into Eq. (15), Eq. (16) can be obtained. 𝚪𝑑𝑖𝑓2𝑎,2 𝐪(𝑡2 ) + 𝚪𝑑𝑖𝑓2𝑏,1 𝐪(𝑡1 ) + 𝚪𝑑𝑖𝑓2𝑐,0 𝐪(𝑡0 ) = 𝚪𝑑𝑖𝑓2𝑑,1 𝐪(𝑡1 − 𝑇) + 𝚪𝑑𝑖𝑓2𝑒,0 𝐪(𝑡0 − 𝑇)

(16)

Where: 𝚪𝑑𝑖𝑓2𝑎,𝑘 = 𝑰,

𝚪𝑑𝑖𝑓2𝑏,𝑘 = −𝑰 −

3𝜏2 2

3𝜏2

𝐑−

2

𝐀𝑘 ,

𝚪𝑑𝑖𝑓2𝑐,𝑘 =

𝜏2 2

𝐑+

𝜏2 2

𝐀𝑘 ,

𝚪𝑑𝑖𝑓2𝑑,𝑘 =

3𝜏2 2

𝐁𝑘 ,

𝚪𝑑𝑖𝑓2𝑒,𝑘 = −

𝜏2 2

𝐁𝑘

(17)

Once again, by using the first-order Adams-Bashforth scheme, 𝐪(𝑡1 ) can be expressed as Eq. (18). (18)

𝐪(𝑡1 ) = 𝐪(𝑡0 ) + 𝜏2 𝐪̇ (𝑡0 ) By substituting Eq. (8) into Eq. (18), Eq. (19) can be obtained. 𝚪𝑑𝑖𝑓1𝑎,1 𝐪(𝑡1 ) + 𝚪𝑑𝑖𝑓1𝑏,0 𝐪(𝑡0 ) = 𝚪𝑑𝑖𝑓1𝑐,0 𝐪(𝑡0 − 𝑇)

(19)

Where: 𝚪𝑑𝑖𝑓1𝑎,𝑘 = 𝑰,

𝚪𝑑𝑖𝑓1𝑏,𝑘 = −𝑰 − 𝜏2 𝐑 − 𝜏2 𝐀𝑘 ,

(20)

𝚪𝑑𝑖𝑓1𝑐,𝑘 = 𝜏2 𝐁𝑘

Then, according to Eq. (10), Eq. (13), Eq. (16) and Eq. (19), the iteration scheme from (𝐪(𝑡𝑚 ), 𝐪(𝑡𝑚−1 ), ⋯ , 𝐪(𝑡0

))𝑇

to

(𝐪(𝑡𝑚 − 𝑇), 𝐪(𝑡𝑚−1 − 𝑇), ⋯ , 𝐪(𝑡0 − 𝑇)) can be expressed as follows: 𝑇

𝐪(𝑡𝑚 ) 𝐪(𝑡𝑚 − 𝑇) ⋮ ⋮ 𝚲𝑑𝑖𝑓3𝑎 ( ) = 𝚲𝑑𝑖𝑓3𝑏 ( ) 𝐪(𝑡1 ) 𝐪(𝑡1 − 𝑇) 𝐪(𝑡0 ) 𝐪(𝑡0 − 𝑇)

(21)

Where:

𝚲𝑑𝑖𝑓3𝑎

𝚪𝑑𝑖𝑓3𝑎,𝑚 𝟎 ⋮ 𝟎 = 𝟎 𝟎 ( 𝟎

𝚲𝑑𝑖𝑓3𝑏 =

𝚪𝑑𝑖𝑓3𝑏,𝑚−1 𝚪𝑑𝑖𝑓3𝑎,𝑚−1 ⋱ ⋯ ⋯ ⋯ ⋯

𝟎 𝟎 ⋮ 𝟎 𝟎 𝟎 (𝑒𝐑(𝑡0−𝑡𝑖𝑛𝑖𝑡𝑖𝑎𝑙)

𝚪𝑑𝑖𝑓3𝑐,𝑚−2 𝚪𝑑𝑖𝑓3𝑏,𝑚−2 ⋱ 𝟎 𝟎 𝟎 𝟎

𝚪𝑑𝑖𝑓3𝑑,𝑚−3 𝚪𝑑𝑖𝑓3𝑐,𝑚−3 ⋱ 𝚪𝑑𝑖𝑓3𝑎,3 𝟎 𝟎 𝟎

𝟎 𝚪𝑑𝑖𝑓3𝑑,𝑚−4 ⋱ 𝚪𝑑𝑖𝑓3𝑏,2 𝚪𝑑𝑖𝑓2𝑎,2 𝟎 𝟎

⋯ ⋱ ⋱

𝚪𝑑𝑖𝑓3𝑐,1 𝚪𝑑𝑖𝑓2𝑏,1 𝚪𝑑𝑖𝑓1𝑎,1 𝟎

𝟎

⋯ ⋱ ⋱

𝟎 ⋮ 𝟎

𝚪𝑑𝑖𝑓3𝑑,0 𝚪𝑑𝑖𝑓2𝑐,0 𝚪𝑑𝑖𝑓1𝑏,0 𝑰 )

𝚪𝑑𝑖𝑓3𝑒,𝑚−1 𝟎 ⋱ ⋱ ⋱ 𝟎

𝚪𝑑𝑖𝑓3𝑓,𝑚−2 𝚪𝑑𝑖𝑓3𝑒,𝑚−2 ⋱ ⋱ ⋱ ⋱

𝚪𝑑𝑖𝑓3𝑔,𝑚−3 𝚪𝑑𝑖𝑓3𝑓,𝑚−3 ⋱ 𝟎 𝟎 𝟎

𝚪𝑑𝑖𝑓3𝑔,𝑚−4 ⋱ 𝚪𝑑𝑖𝑓3𝑒,2 𝟎 𝟎

𝚪𝑑𝑖𝑓3𝑓,1 𝚪𝑑𝑖𝑓2𝑑,1 𝟎

𝚪𝑑𝑖𝑓3𝑔,0 𝚪𝑑𝑖𝑓2𝑒,0 𝚪𝑑𝑖𝑓1𝑐,0

𝟎

⋯

𝟎

𝟎

𝟎

𝟎

(22)

𝟎 ⋮ 𝟎

)

Then, the Floquet Matrix in third-order numerical difference method can be expressed as follows:

6

𝐏𝑑𝑖𝑓3 = 𝚲𝑑𝑖𝑓3𝑎 −1 𝚲𝑑𝑖𝑓3𝑏

(23)

Similarly, the Floquet Matrix can also be derived from the forth-order, the fifth-order and the sixth-order Adams-Bashforth scheme, simply by updating Eq. (11), (15) and (18). As shown in Fig. 3, the X axis represents 𝑛𝑠 , then 𝑇 = 𝑛𝑠 ⁄(60𝑁). The Y axis represents 𝑎𝑝 , and the Z axis represents the maximum module of 𝐏𝑑𝑖𝑓3 's eigenvalues. Then, the contour with a height of 1.0 is the critical stability boundary.

Fig. 3. Distribution of the maximum module of Floquet Matrix's eigenvalues 2.3 Direct numerical difference scheme and direct numerical integration scheme piecewise numerical difference method and piecewise numerical integration method divides one complete tooth passing period into a cutting part and a non-cutting part. In fact, the non-cutting part is a special form of the cutting part, and the cutting part is the general form of the non-cutting part. This means that in the non-cutting part, the numerical difference scheme and the numerical integration scheme are still valid. Therefore, if one tooth passing period T is discretized into m equal intervals and the span of each interval is τ, then piecewise numerical difference method and piecewise numerical integration method becomes direct numerical difference scheme and direct numerical integration scheme. In direct numerical difference scheme, the 𝚲𝑑𝑖𝑓3𝑏 in Eq. (22) becomes Eq. (24). Other equations don’t change but the meaning of τ in other equations becomes 𝑇/𝑚.

𝚲𝑑𝑖𝑓3𝑏

𝟎 𝚪𝑑𝑖𝑓3𝑒,𝑚−1 𝟎 𝟎 ⋮ ⋱ ⋱ = 𝟎 𝟎 ⋱ 𝟎 𝟎 (𝐈 𝟎

𝚪𝑑𝑖𝑓3𝑓,𝑚−2 𝚪𝑑𝑖𝑓3𝑒,𝑚−2 ⋱ ⋱ ⋱ ⋱ ⋯

𝚪𝑑𝑖𝑓3𝑔,𝑚−3 𝚪𝑑𝑖𝑓3𝑓,𝑚−3 ⋱ 𝟎 𝟎 𝟎 𝟎

𝟎 𝚪𝑑𝑖𝑓3𝑔,𝑚−4 ⋱ 𝚪𝑑𝑖𝑓3𝑒,2 𝟎 𝟎 𝟎

⋯ ⋱ ⋱

𝚪𝑑𝑖𝑓3𝑓,1 𝚪𝑑𝑖𝑓2𝑑,1 𝟎 𝟎

𝟎 ⋮ 𝟎

𝚪𝑑𝑖𝑓3𝑔,0 𝚪𝑑𝑖𝑓2𝑒,0 𝚪𝑑𝑖𝑓1𝑐,0 𝟎 )

(24)

Similar to Eq. (24), direct numerical integration scheme can also be derived by modifying the initial conditions. 2.4 Multiple modal time domain method This paper makes a simplifying assumption for multi-modal problems. The ith-order natural frequency is assumed to be independent of the jth-order natural frequency. i.e. the first-order natural frequency is assumed to be independent of the second-order natural frequency. Similarly, stiffness is also assumed to be independent. And damping ratio is also assumed to be independent. Then, by combining the modal parameters in the X direction and the modal parameters in the Y direction, Eq. (7) can be extended along main diagonal. i.e. If the modal masses in the X direction are 𝑚𝑥1 and 𝑚𝑥2 , and the modal masses in the Y direction are 𝑚𝑦1 and 𝑚𝑦2 , and the damping in the X direction is 𝑐𝑥1 and 𝑐𝑥2 , and the damping in the Y direction is 𝑐𝑦1 and 𝑐𝑦2 , and the stiffness in the X direction is 𝑘𝑥1 and 𝑘𝑥2 , and the stiffness in the Y direction is 𝑘𝑦1 and 𝑘𝑦2 , then Eq. (7) can be extended to Eq. (25).

7

𝑚𝑥1 0 0 0 𝐌= 0 0 0 ( 0 𝑘𝑥1 0 0 0 𝐊= 0 0 0 ( 0

0 𝑚𝑦1 0 0 0 0 0 0 0 𝑘𝑦1 0 0 0 0 0 0

0 0 𝑚𝑥1 0 0 0 0 0 0 0 𝑘𝑥1 0 0 0 0 0

0 0 0 𝑚𝑦2 0 0 0 0 0 0 0 𝑘𝑦2 0 0 0 0

0 0 0 0 𝑚𝑥2 0 0 0

0 0 0 0 𝑘𝑥2 0 0 0

0 0 0 0 0 𝑚𝑦1 0 0

0 0 0 0 0 𝑘𝑦1 0 0

0 0 0 0 0 0 𝑘𝑥2 0

0 0 0 0 0 0 𝑚𝑥2 0

0 𝑐𝑥1 0 0 0 0 0 0 , 𝐂= 0 0 0 0 0 0 𝑚𝑦2 ) (0

0 0 0 0 , 0 0 0 𝑘𝑦2 )

0 𝑐𝑦1 0 0 0 0 0 0

0 0 𝑐𝑥1 0 0 0 0 0

0 0 0 𝑐𝑦2 0 0 0 0

0 0 0 0 𝑐𝑥2 0 0 0

0 0 0 0 0 𝑐𝑦1 0 0

0 0 0 0 0 0 𝑐𝑥2 0

0 0 0 0 , 0 0 0 𝑐𝑦2 )

𝑥11 (𝑡) 𝑦11 (𝑡) 𝑥12 (𝑡) 𝑦12 (𝑡) 𝐬(𝑡) = , 𝑥21 (𝑡) 𝑦21 (𝑡) 𝑥22 (𝑡) (𝑦22 (𝑡))

(25)

𝑎𝑥𝑥 (𝑡) 𝑎𝑥𝑦 (𝑡) 0 0 0 0 0 0 𝑥11 (𝑡) − 𝑥11 (𝑡 − 𝑇) 𝑎𝑦𝑥 (𝑡) 𝑎𝑦𝑦 (𝑡) 0 0 0 0 0 0 𝑦11 (𝑡) − 𝑦11 (𝑡 − 𝑇) 0 0 𝑎𝑥𝑥 (𝑡) 𝑎𝑥𝑦 (𝑡) 0 0 0 0 𝑥12 (𝑡) − 𝑥12 (𝑡 − 𝑇) 0 0 𝑎𝑦𝑥 (𝑡) 𝑎𝑦𝑦 (𝑡) 0 0 0 0 𝑦12 (𝑡) − 𝑦12 (𝑡 − 𝑇) 1 𝐅(𝐬(𝑡), 𝐬(𝑡 − 𝑇), 𝑡) = 𝑎𝑝 𝐾𝑡 2 0 0 0 0 𝑎𝑥𝑥 (𝑡) 𝑎𝑥𝑦 (𝑡) 0 0 𝑥21 (𝑡) − 𝑥21 (𝑡 − 𝑇) 𝑦21 (𝑡) − 𝑦21 (𝑡 − 𝑇) 0 0 0 0 𝑎𝑦𝑥 (𝑡) 𝑎𝑦𝑦 (𝑡) 0 0 𝑥22 (𝑡) − 𝑥22 (𝑡 − 𝑇) 0 0 0 0 0 0 𝑎𝑥𝑥 (𝑡) 𝑎𝑥𝑦 (𝑡) 0 0 0 0 0 𝑎𝑦𝑥 (𝑡) 𝑎𝑦𝑦 (𝑡)) (𝑦22 (𝑡) − 𝑦22 (𝑡 − 𝑇)) ( 0 Then, the numerical methods can be applied to Eq. (6), which have been discussed in 2.2, 2.3 and Ref. [12], [14], [15]. 2.5 Time complexity analysis of these methods Fig. 4 summarizes the implementation of the numerical methods. In numerical methods, there are five factors that affect the time complexity, which are listed in Tab. 1. Start

P4: The 2nd loop.

ns = ns,min

ns <= n s,max ?

No

No

Yes ns = ns + Δns ap = ap,min

P3: The 1st loop.

Draw the stability lobes diagram by using all saved eigenvalues.

ap <= a p,max?

P2: The most time consuming steps.

Yes ap = ap + Δap

Construct Λa and Λb. (In difference method in this paper; Similar matrix in Integration method In Eq. (17), (21), (29) in Ref. [15].) or Construct D0, D1, ... Dm-1. (In semi-discretization method in Eq. (52) in Ref. [12]; In Full-discretization method in Eq. (25) in Ref. [14].)

P1.1: Construct matrix

P1.2: Calculate Floquet matrix.

P1.3: Calculate eigenvalues.

Calculate the Floquet matrix.

Calculate the eigenvalues of the Floquet matrix. Then, find and save the eigenvalue with the largest module.

Fig. 4. Flow chart of numerical methods

8

End

Tab. 1. Factors that affect the time complexity of numerical methods 𝑓𝑛𝑎 :

Number of iterations of the spindle speed. For example, in Fig. 4, 𝑓𝑛𝑎 = (𝑛𝑠,𝑚𝑎𝑥 − 𝑛𝑠,𝑚𝑖𝑛 )/Δ𝑛𝑠 + 1.

𝑓𝑛𝑏 :

Number of iterations of axial depth of cut. For example, in Fig. 4, 𝑓𝑛𝑏 = (𝑎𝑝,𝑚𝑎𝑥 − 𝑎𝑝,𝑚𝑖𝑛 )/Δ𝑎𝑝 + 1.

𝑓𝑛𝑐 :

In semi-discretization method, full-discretization method, direct numerical difference scheme and direct numerical integration scheme, they discretize one tooth passing period T into 𝑓𝑛𝑐 equal intervals. In piecewise numerical integration method and piecewise numerical difference method, they discretize the Part B in Fig. 2 (b) into 𝑓𝑛𝑐 equal intervals.

𝑓𝑛𝑑 :

Number of modal parameters in the X direction.

𝑓𝑛𝑒 :

Number of modal parameters in the Y direction. In Fig. 4, in P1.1, in numerical difference methods, matrix addition is the most time consuming step with a time complexity

of O(𝑓𝑛𝑑 2 × 𝑓𝑛𝑒 2 × 𝑓𝑛𝑐 ); in numerical integration methods, matrix multiplication is the most time consuming step with a time complexity of O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 ); in semi-discretization method and full-discretization method, if LU decomposition is used, both matrix inversion and matrix multiplication are the most time consuming step with a time complexity of O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 ). In P1.2, in numerical difference methods and numerical integration methods, if LU decomposition is used, both the time complexity of matrix inversion and matrix multiplication are O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 3 ) ; in semi-discretization method and full-discretization method, the time complexity of continuous matrix multiplication is O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 4 ). In P1.3, in all numerical methods, if generalized Schur decomposition is used, the time complexity of solving matrix eigenvalues is O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 3 ). Therefore, in P2, numerical difference methods have a time complexity of O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 3 ); numerical integration methods have a time complexity of O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 3 ); semi-discretization method has a time complexity of O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 4 ); full-discretization method has a time complexity of O(𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 4 ). Then, in P3, all numerical difference methods have a time complexity of O(𝑓𝑛𝑏 ). Then, in P4, all numerical difference methods have a time complexity of O(𝑓𝑛𝑎 ). Finally, by considering P1.1, P1.2, P1.3, P2, P3 and P4, the time complexity of each numerical method is listed in Tab. 2. Tab. 2. Time complexity of each numerical method Method

Time complexity

Semi-discretization method

O(𝑓𝑛𝑎 × 𝑓𝑛𝑏 × 𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 4 )

Full-discretization method

O(𝑓𝑛𝑎 × 𝑓𝑛𝑏 × 𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 4 )

Piecewise numerical integration method

O(𝑓𝑛𝑎 × 𝑓𝑛𝑏 × 𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 3 )

Direct numerical integration scheme

O(𝑓𝑛𝑎 × 𝑓𝑛𝑏 × 𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 3 )

Piecewise numerical difference method

O(𝑓𝑛𝑎 × 𝑓𝑛𝑏 × 𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 3 )

Direct numerical difference scheme

O(𝑓𝑛𝑎 × 𝑓𝑛𝑏 × 𝑓𝑛𝑑 3 × 𝑓𝑛𝑒 3 × 𝑓𝑛𝑐 3 )

Fig. 5 summarizes the implementation of the frequency domain analytical method. In analytical method, there are three factors that affect the time complexity, which are listed in Tab. 3.

9

Start

P4: The 2nd loop.

ns = ns,min

ns <= n s,max ?

No

No

Yes ns = ns + Δns ωc = ωc,min

P3: The 1st loop.

Draw the stability lobes diagram by using all saved eigenvalues.

ωc <= ωc,max?

P2: The most time consuming steps.

Yes

End

ωc = ωc + Δωc

P1.1: Construct matrix

Construct the matrix G(ωc,ωT) in Eq. (31) in Ref. [13].

Calculate all eigenvalues of the matrix G(ωc,ωT) in Eq. (31) in Ref. [13].

P1.2: Calculate eigenvalues.

P1.3: Filter eigenvalues.

Filter each eigenvalue one by one. First, calculate ap by using each eigenvalue. Then, select and save the eigenvalue that enables the real part of ap to be positive and the imaginary part of ap to be closest to zero.

Fig. 5. Flow chart of analytical method Tab. 3. Factors that affect the time complexity of analytical method 𝑓𝑎𝑎 :

Number of iterations of the spindle speed. For example, in Fig. 5, 𝑓𝑎𝑎 = (𝑛𝑠,𝑚𝑎𝑥 − 𝑛𝑠,𝑚𝑖𝑛 )/Δ𝑛𝑠 + 1.

𝑓𝑎𝑏 :

Number of iterations of chatter frequency. For example, in Fig. 5, 𝑓𝑎𝑏 = (𝜔𝑐,𝑚𝑎𝑥 − 𝜔𝑐,𝑚𝑖𝑛 )/Δ𝜔𝑐 + 1.

𝑓𝑎𝑐 :

The number of rows of the matrix 𝐆(𝜔𝑐 , 𝜔 𝑇 ) in Eq. (31) in Ref. [13], which is determined by the number of terms of Fourier expansion. In Fig. 5, in P1.1, multi-frequency solution method constructs a 𝑓𝑎𝑐 -order matrix Ω. The time complexity of this step is

O(𝑓𝑎𝑐 2 ). In P1.2, if generalized Schur decomposition is used, the time complexity of solving matrix eigenvalues is O(𝑓𝑎𝑐 3 ). In P1.3, multi-frequency solution method needs to traverse the 𝑓𝑎𝑐 eigenvalues. Therefore, the time complexity of this step is O(𝑓𝑎𝑐 ). Therefore, in P2, multi-frequency solution method has a time complexity of O(𝑓𝑎𝑐 3 ). Then, in P3, multi-frequency solution method has a time complexity of O(𝑓𝑎𝑏 ). Then, in P4, multi-frequency solution method has a time complexity of O(𝑓𝑎𝑎 ). Finally, by considering P1.1, P1.2, P1.3, P2, P3 and P4, the time complexity of multi-frequency solution method is O(𝑓𝑎𝑎 × 𝑓𝑎𝑏 × 𝑓𝑎𝑐 3 ). Then, the time complexity models are verified by measuring the elapsed time on several CPUs. The measured elapsed time is shown in Fig. 6. As can be seen from Fig. 6, the measured elapsed time grow linearly. Then, the measured elapsed time is linearly fitted and the results of the fitting are shown in Tab. 4. In Tab. 4, all R-squares are very close to 1, therefore, the time complexity models are consistent with the measured elapsed time.

10

(a) Measured elapsed time of semi-discretization method

(b) Measured elapsed time of full-discretization method

(c) Measured elapsed time of piecewise numerical difference

(d) Measured elapsed time of piecewise numerical

method

integration method

(e) Measured elapsed time of direct numerical difference

(f) Measured elapsed time of direct numerical integration

scheme.

scheme

(g) Measured elapsed time of multi-frequency solution method Fig. 6. Measured elapsed time on several CPUs Tab. 4. Linear fitting to the measured elapsed time Method

The slope of the fitted line

R-square

CPU

Semi-discretization method

7.38e-09

0.99

Intel Pentium [email protected]

Full-discretization method

1.92e-09

0.99

Intel Core i7 [email protected]

Piecewise numerical integration method

2.62e-08

0.98

Intel Core i7 [email protected]

Direct numerical integration scheme

2.15e-08

0.97

Intel Core i7 [email protected]

Piecewise numerical difference method

1.73e-08

0.99

Intel Core i7 [email protected]

Direct numerical difference scheme

1.94e-08

0.99

Intel Core i7 [email protected]

Multi-frequency solution method

5.74e-07

0.94

Intel Core i5 [email protected]

2.6 Parallel algorithms of these methods Parallel algorithms can be applied to analytical methods and numerical methods. A parallel algorithm which can reach linear speedup will be presented in this section. In Fig. 4 and Fig. 5, P2 run repeatedly. One iteration consists of P1.1, P1.2, and P1.3. However, the input parameters of one iteration do not depend on the results of any other iterations. This means that several iterations can run independently at the same time. Therefore, several iterations can be assigned to several threads,

11

respectively. The runtime of one iteration is almost equal to the runtime of the other iterations. Then several threads run in parallel, causing the total runtime to reduce linearly. Finally, these threads should be synchronized in order to obtain the whole stability border. Tab. 5. Measured elapsed time of multi-frequency solution method and semi-discretization method The number of

Measured elapsed time of multi-frequency solution

Measured elapsed time of semi-discretization

threads

method (seconds)

method (seconds)

1

460.35

271.35

2

218.50

131.97

3

161.80

100.90

4

134.65

89.79

Then, this paper takes multi-frequency solution method and semi-discretization method as examples to compare the effects of parallel algorithms on the runtime of the methods. The parallel algorithms are verified by measuring the elapsed time on the Intel Core i7 [email protected] CPU. The measured elapsed time is listed in Tab. 5. As can be seen from Tab. 5, elapsed time of multi-frequency solution method decreases linearly with the increase of the number of threads. Elapsed time of semi-discretization method doesn’t decrease linearly but not very bad. 3. Predicting stability lobes 3.1 Machining condition Real machining process can be predicted by simulations. Stability lobes will be generated in simulations. Relations between the stable cutting depth and rotational speed of spindle can be displayed in a stability lobes diagram. The simulations are based on a CNC machining center whose detailed information is listed in Tab. 6. The workpiece is made of Ti-6Al-4V (TC4). The milling-tool with flat end is made of tungsten steel. There are 4 flutes on the cutter. Mode of milling is down milling. Diameter of cutter is 10mm. Cutting depth in radial direction is 0.5mm. 𝐾𝑡 is the cutting force coefficient in tangential direction. 𝐾𝑟 is the ratio of radial cutting force coefficient to tangential one. Both 𝐾𝑡 and 𝐾𝑟 are obtained by the milling force experiment. The average cutting force during one tooth period in full-immersion milling can be expressed as follow [39]: 𝑁 × 𝑎𝑝 𝑁 × 𝑎𝑝 × 𝐾𝑟𝑐 × 𝑓 − × 𝐾𝑟𝑒 4 𝜋 𝑁 × 𝑎𝑝 𝑁 × 𝑎𝑝 𝐹̅𝑦 = + × 𝐾𝑡𝑐 × 𝑓 + × 𝐾𝑡𝑒 4 𝜋 𝑁 × 𝑎𝑝 𝑁 × 𝑎𝑝 ̅ {𝐹𝑧 = + 𝜋 × 𝐾𝑎𝑐 × 𝑓 − 2 × 𝐾𝑎𝑒 𝐹̅𝑥 = −

(26)

The feed rate f can be controlled in milling experiment. And then the instantaneous cutting force can be measured by a dynamometer. Then average cutting force Fx , Fy and Fz can be figured out. Corresponding to several feed rate f1 , f2 , f3 ,⋯, fn , average cutting force will be Fx1 , Fx2 , Fx3 , ⋯, Fxn . Then, by using linear regression, K rc and K re can be figured out. K tc, K te , K ac and K ae can be figured by the same method. And then K t and K r can be expressed as follow [39]: {

𝐾𝑡 = 𝐾𝑡𝑐 𝐾𝑟 = 𝐾𝑟𝑐 ⁄𝐾𝑡𝑐

(27)

In this paper, 𝐾𝑡 is 1773.10N/𝑚𝑚2. 𝐾𝑟 is 0.31. 3.2 Modal test Hammer modal test is conducted to obtain modal parameters of the cutter and the workpiece. This test system is composed of a CNC machining center, an acceleration transducer, a data acquisition card, a PC, a hammer, a milling-tool and a thin-walled workpiece. Details of these facilities are shown in Tab. 6 and Fig. 7.

12

(a)

(b)

CNC

Workpiece Spindle

Cutter Clamper

Dynamometer

(c)

(d) Acceleration

Impact

transducer

hammer

Data acquisition card Computer

(a) CNC machining center TH5650 (b) Measuring cutting force (c) Data sampling system (d) Modal test Fig. 7. Facilities to measure modal parameters and cutting force Tab. 6. Detail of facilities Facility

Detail

CNC machining center

TH5650 vertical boring milling machining center

acceleration transducer

Kistler 8778A500

PC

OS: Microsoft Windows XP software: CutPRO

hammer

Kistler 9722A500

Fig. 8. FRF of cutter in X-direction

13

Fig. 9. FRF of cutter in Y-direction The raw FRFs which are measured by modal test are presented in Fig. 8 and Fig. 9. Then, modal parameters can be identified by follow steps: a)

Natural frequency locates on the intersection point between real part of FRF curve and abscissa axis. For example, in Fig. 8, the first natural frequency ωn1 is 1070.05Hz.

b)

Two frequencies near to natural frequency should be identified. The first one locates on the wave crest nearest to natural frequency. For example, it means the ωa1 in Fig. 8. Another one locates on the wave trough nearest to natural frequency. For example, it means the ωb1 in Fig. 8.

c)

Then damping ratio ζ1 equals to (ωb1 − ωa1 )⁄(2 × ωn1 ) = (1122 − 1040)⁄(2 × 1070.05) = 0.038.

d)

Amplitude of imaginary part of FRF curve should be identified, for example, it means the Imp1 in Fig. 8. Then Stiffness equals to −1⁄(2 × ζ1 × Imp1 ) = −1⁄(2 × 0.038 × (−7.59 × 10−7 )) = 1.72 × 107 . Thus, all modal parameters of the cutter can be shown in Tab. 7. Then, FRFs which are generated by modal parameters in

Tab. 7 are presented in Fig. 8 and Fig. 9. In Fig. 8 and Fig. 9, Experiment-FRFs are almost the same to Fitting-FRFs. Tab. 7. Modal parameters of the milling-tool X-direction

Y-direction

1-order natural frequency (Hz)

1070.05

1032.95

1-order stiffness (N/m)

17199700.00

67991000.00

1-order damping ratio

3.83e-2

2.32e-2

2-order natural frequency (Hz)

1721.21

1845.99

2-order stiffness (N/m)

76377700.00

103122000.00

2-order damping ratio

1.95e-2

1.35e-2

Since modal parameters of thin-walled workpiece will change dramatically during the cutting process, modal parameters of the workpiece are obtained at several different moments during the whole machining process. In stage 0, the acceleration transducer is pasted on point A1 in Fig. 10 in order to approach the real excitation point of milling operation. And hammer will also knock at point A1. Then, in stage 1, the acceleration transducer is pasted on point A2 in Fig. 10. Hammer knocks at point A2. Then, in stage2, acceleration transducer is pasted on point A3 in Fig. 10. Hammer knocks at point A3. The chemical composition of material Ti-6Al-4V (ISO5832-3) of workpiece is listed in Tab. 9. Modal parameters of the workpiece are shown in Tab. 8. And shapes of the workpiece at those moments are shown in Fig. 10. Tab. 8. Modal parameters of the workpiece stage 0

stage 1

stage 2

1-order natural frequency (Hz)

343.54

357.22

408.73

1-order stiffness (N/m)

2106100.00

683111.00

3874170.00

1-order damping ratio

4.80e-2

4.90e-2

6.12e-2

2-order natural frequency (Hz)

N/A

N/A

504.00

2-order stiffness (N/m)

N/A

N/A

10244800.00

2-order damping ratio

N/A

N/A

2.88e-2

3-order natural frequency (Hz)

N/A

N/A

2292.14

3-order stiffness (N/m)

N/A

N/A

46962000.00

3-order damping ratio

N/A

N/A

1.46e-2

4-order natural frequency (Hz)

N/A

N/A

4623.64

4-order stiffness (N/m)

N/A

N/A

445815000.00

4-order damping ratio

N/A

N/A

1.73e-3

14

Fig. 10. Shapes of the workpiece 3.3 Stability lobes in milling process Milling modeling is a self-excited vibration model[3]. Therefore, in frequency domain, a frequency response function (FRF) represents the ratio of the vibration displacement to the cutting force. The system’s transfer function in Fig. 1 is equal to cutter’s transfer function plus workpiece’s transfer function[25]. Therefore, by observing Fig. 11 and Fig. 12, it can be roughly analyzed whether the dynamic characteristics of the milling system is mainly affected by the cutter or by the workpiece. It can be seen from Fig. 11 that when the frequency is low, the dynamic characteristics of the system are mainly affected by the workpiece; when the frequency is high, the dynamic characteristics of the system are mainly affected by the cutter. Then, as the material is removed, the dynamic characteristics of the system change. In Fig. 12, when the frequency is low, the dynamic characteristics of the system are equally affected by the cutter and the workpiece; when the frequency is high, the dynamic characteristics of the system are mainly affected by the workpiece.

(a) Amplitude of FRFs

(b) Phase of FRFs

Fig. 11. Cutter’s FRF, workpiece’s FRF and system’s FRF in Y direction in stage 0

(a) Amplitude of FRFs

(b) Phase of FRFs

Fig. 12. Cutter’s FRF, workpiece’s FRF and system’s FRF in Y direction in stage 2 The critical stability boundaries are predicted by using ZOA method, multi-frequency solution method, semi-discretization method, full-discretization method, direct numerical difference scheme, piecewise numerical difference method, direct numerical integration scheme, and piecewise numerical integration method, respectively. To describe the predicted boundaries, the following functions are defined to represent the predicted critical stability boundaries.

15

𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 𝑎𝑝,𝑚𝑎𝑥 {𝑎𝑝,𝑚𝑎𝑥

= 𝑕1 (𝑛𝑠 ) = 𝑕2 (𝑛𝑠 ) = 𝑕3 (𝑛𝑠 ) = 𝑕4 (𝑛𝑠 ) = 𝑕5 (𝑛𝑠 ) = 𝑕6 (𝑛𝑠 ) = 𝑕7 (𝑛𝑠 ) = 𝑕8 (𝑛𝑠 ) = 𝑕9 (𝑛𝑠 ) = 𝑕10 (𝑛𝑠 ) = 𝑕11 (𝑛𝑠 ) = 𝑕12 (𝑛𝑠 ) = 𝑕13 (𝑛𝑠 ) = 𝑕14 (𝑛𝑠 )

, Means the critical stability cutting depth predicted by using multi − frequency solution method. , Means the critical stability cutting depth predicted by using ZOA. , Means the critical stability cutting depth predicted by using semi − discretization method. , Means the critical stability cutting depth predicted by using full − discretization method. , Means the critical stability cutting depth predicted by using direct numerical difference scheme 3rd order. , Means the critical stability cutting depth predicted by using direct numerical difference scheme 4th order. , Means the critical stability cutting depth predicted by using direct numerical difference scheme 5th order. , Means the critical stability cutting depth predicted by using direct numerical difference scheme 6th order. , Means the critical stability cutting depth predicted by using piecewise numerical difference method 3rd order. , Means the critical stability cutting depth predicted by using piecewise numerical difference method 4th order. , Means the critical stability cutting depth predicted by using piecewise numerical difference method 5th order. , Means the critical stability cutting depth predicted by using piecewise numerical difference method 6th order. , Means the critical stability cutting depth predicted by using direct numerical integration scheme. , Means the critical stability cutting depth predicted by using piecewise numerical integration method.

(28)

(b) Absolute differences among numerical solutions and analytical solutions

(c) Relative differences among numerical solutions and analytical solutions

(a) Stability lobes

(d) Absolute differences among numerical solutions

(e) Relative differences among numerical solutions

Fig. 13. Predicting critical stability boundaries in stage 0 using numerical methods and analytical methods

The critical stability boundaries during stage 0 are shown in Fig. 13 (a). The area above the boundary is unstable, while the area below the boundary is stable. Some general characteristics can be obtained from Fig. 13 (a). Relationship between rotational speed of spindle and border of the stable cutting depth is not monotonous. When rotational speed of spindle is growing continuously, the critical stable cutting depth increases and decreases alternately. This trend is like a wave. When rotational speed of spindle is not very high, overall, the critical stable cutting depth tends to increase with the growth of rotational speed of spindle. With the growth of rotational speed spindle, lobes tend to become sparser. In other words, Lobes will be very dense when rotational speed of spindle is very low. The absolute differences among the critical stability boundaries predicted by the numerical methods, which means 𝑕3 (𝑛𝑠 ) minus 𝑕14 (𝑛𝑠 ), or 𝑕4 (𝑛𝑠 ) minus 𝑕14 (𝑛𝑠 ), etc., are shown in Fig. 13 (d). The relative differences among the critical stability boundaries predicted by the numerical methods, which means [𝑕3 (𝑛𝑠 ) − 𝑕14 (𝑛𝑠 )]⁄𝑕14 (𝑛𝑠 ), or [𝑕4 (𝑛𝑠 ) − 𝑕14 (𝑛𝑠 )]⁄𝑕14 (𝑛𝑠 ), etc., are shown in Fig. 13 (e). As can be seen from Fig. 13 (d) and (e), the critical stability boundaries predicted by the numerical methods are very consistent. In Fig. 13 (d), the absolute differences are substantially less than 0.06 mm. And in Fig. 13 (e), the relative

16

differences are substantially less than 3%. The absolute differences between the critical stability boundary predicted by the numerical method and the stability boundaries predicted by the analytical methods, which means 𝑕1 (𝑛𝑠 ) minus 𝑕14 (𝑛𝑠 ), or 𝑕2 (𝑛𝑠 ) minus 𝑕14 (𝑛𝑠 ), are shown in Fig. 13 (b). The relative differences between the critical stability boundary predicted by the numerical method and the stability

boundaries predicted by the analytical methods, which means [𝑕1 (𝑛𝑠 ) − 𝑕14 (𝑛𝑠 )]⁄𝑕14 (𝑛𝑠 ), or [𝑕2 (𝑛𝑠 ) − 𝑕14 (𝑛𝑠 )]⁄𝑕14 (𝑛𝑠 ), are shown in Fig. 13 (c). As can be seen from Fig. 13 (b) and (c), when the spindle speed is smaller than 10000 rpm, the critical stability boundary predicted by the numerical method is also consistent with the critical stability boundaries predicted by the analytical methods. In Fig. 13 (b) and (c), when the spindle speed is greater than 10000 rpm, the difference between the critical stability boundary predicted by ZOA and the critical stability boundary predicted by piecewise numerical integration method begins to increase, but this difference is caused by the low approximation precision. Multi-frequency solution method uses a higher order Fourier series expansion than ZOA. Therefore, as can be seen from Fig. 13 (b) and (c), when the approximation precision is improved, the numerical methods and the analytical methods are still consistent.

(b) Absolute differences among numerical solutions and analytical solutions

(c) Relative differences among numerical solutions and analytical solutions

(a) Stability lobes

(d) Absolute differences among numerical solutions

(e) Relative differences among numerical solutions

Fig. 14. Predicting critical stability boundaries in stage 2 using numerical methods and analytical methods

The critical stability boundaries during stage 2 are shown in Fig. 14 (a). The absolute differences and relative differences among the critical stability boundaries predicted by the numerical methods are shown in Fig. 14 (d) and (e), respectively. The critical stability boundaries predicted by the numerical methods are also very consistent in stage 2. In Fig. 14 (b) and (c), the differences between the critical stability boundaries predicted by the analytical methods and the critical stability boundaries predicted by the numerical methods tend to increase. However, in Fig. 14 (a), the critical stability boundaries predicted by analytical methods is very similar to the contour of the critical stability boundary predicted by numerical methods. A Stability Lobes Diagram physically means the maximum value that 𝑎𝑝 can reach at a given 𝑛𝑠 , if no chatter occurs. This means that for a milling stability prediction algorithm, 𝑛𝑠 is the input to the algorithm and 𝑎𝑝,𝑚𝑎𝑥 is the output of the algorithm. Therefore, the sensitivity of 𝑎𝑝,𝑚𝑎𝑥 to 𝑛𝑠 can be used to measure the robustness of an algorithm. A milling stability prediction algorithm is robust if its sensitivity is always finite. The sensitivity of a milling stability prediction algorithm can be defined as 𝑑𝑎𝑝,𝑚𝑎𝑥 /𝑑𝑛𝑠 [40]. Sensitivity of piecewise numerical difference method in stage 0 is shown in Fig. 15 (a). The sensitivity is

17

approximated by using central difference scheme. Error of central difference scheme is shown in Fig. 15 (b). As can be seen from Fig. 15 (a), the sensitivity of piecewise numerical difference method is finite and its range is (−0.1, 0.05). Therefore, the impact of small changes in 𝑛𝑠 on 𝑎𝑝,𝑚𝑎𝑥 is finite. Therefore, the piecewise numerical difference method is robust.

(a) Sensitivity of piecewise numerical difference

(b) Error of central difference

method in stage 0 Fig. 15. Sensitivity of piecewise numerical difference method In the milling stability prediction algorithm proposed in this paper, the Adams-Bashforth scheme is used to approximate the solution of the differential equation. For example, if the 3rd-order Adams-Bashforth scheme is applied to the algorithm, the local truncation error is 𝑂[60(𝜑𝑒𝑥 − 𝜑𝑠𝑡 )⁄(2𝜋𝑛𝑠 𝑓𝑛𝑐 )]4 . This means that as 𝑛𝑠 increases or 𝑓𝑛𝑐 increases, the local truncation error will decrease. Then, the reduction of the local truncation error will reduce the error of the eigenvalues of 𝐏𝑑𝑖𝑓3 in Eq. (23), thereby reducing the error of the critical stability boundary prediction. In Fig. 16 (a), as 𝑛𝑠 decreases or 𝑓𝑛𝑐 decreases, the error of the critical stability depth of cut increases. In an extreme case, when 𝑛𝑠 is less than 2750rpm and 𝑓𝑛𝑐 is less than 25, the algorithm has failed. If 𝑓𝑛𝑐 increases to 200% (i.e. increases from 25 to 50), then critically expired 𝑛𝑠 is reduced to 50% (i.e. from 2750 to 1375), but 𝑛𝑠 ∙ 𝑓𝑛𝑐 is still 68750. Since the Fourier series is used in multi-frequency solution method, there is also a truncation error when the multi-frequency solution method is implemented using a computer. As can be seen from Fig. 16 (b), the error can be reduced by reducing Δ𝜔 in Fig. 5 or increasing 𝑓𝑎𝑐 .

(a) Lobes predicted by the 3rd-order

(b) Lobes predicted by multi-frequency solution method in stage 0

Adams-Bashforth scheme piecewise numerical difference method Fig. 16. The effect of the truncation error on the predicted lobes.

18

(a) 3D view of the changing trend

Fig. 17. The critical stability boundary is changing with the material removal The critical stability boundary varies with the material removal. The change in the critical stability boundary from stage 0 to stage 1 is shown in Fig. 17. As can be seen in Fig. 17 (a), with the material removal, the critical stability boundary tends to decrease first and then tends to rise. This trend can be seen more clearly in Fig. 17 (b) and Fig. 17 (c). The height difference between stage 1 and stage 0 is almost always negative. And the height difference between stage 2 and stage 1 is almost always positive. The removal of the material can reduce the thickness of the thin-walled workpiece, which will cause a decrease in the stiffness of the workpiece and then cause a drop in the critical stability boundary. On the other hand, with the material removal, the tool-tip is closer to the bottom of the workpiece. The stiffness of the bottom of the workpiece is greater, which causes an increase in the critical stability boundary. Under the combined influence of these two factors, the trend in Fig. 17 (a) appears. 4. Verification 4.1 Experiment setup An experiment is conducted to verify the lobs in Fig. 13 (a). The workpiece in verification experiment is the same to the workpiece in simulation. The chemical composition of material Ti-6Al-4V (ISO5832-3) of workpiece is listed in Tab. 9. Shape of the workpiece is presented at (a) in Fig. 10. And the way to clamp workpiece is shown in (b) in Fig. 7. Six down milling test is performed on the CNC machining center which has been used in simulation. The machining center is shown in (a) in Fig. 7. The cutter with flat end is made of tungsten steel. There are 4 flutes on the cutter. Diameter of cutter is 10mm. The cutter is shown in (b) in Fig. 7. Signal of Cutting force in Y-direction is measured to infer if chatter occurs. Signal of force is measured by Kistler dynamometer 9257B with the sampling rate 3000Hz. The dynamometer is shown in (b) in Fig. 7. The force signal will be transmitted to the computer. By analyzing the data with the computer, cutting forces in X-direction, Y-direction and Z-direction at any moment can be obtained and saved as a file. Only cutting forces in Y-direction are useful. Tab. 9. Chemical composition of Ti-6Al-4V (ISO5832-3) Element

Ti

Al

V

Fe

C

N

H

O

Others

Mass (%)

Base

5.50-6.75

3.50-4.50

< 0.30

< 0.08

< 0.05

< 0.015

< 0.20

< 0.40

Six cases being tested is presented in Fig. 18 (d). In case A, rotational speed of spindle is 2600 rpm. Cutting depth in direction of spindle is 1mm. Cutting depth in radial direction is 0.5mm. The cutter feed at 0.005 mm/tooth. In case B, rotational speed of spindle is 2600 rpm. Cutting depth in direction of spindle is 2mm. Cutting depth in radial direction is 0.5mm. The cutter feed at 0.005 mm/tooth. In case C, rotational speed of spindle is 2600 rpm. Cutting depth in direction of spindle is 3.2mm. Cutting depth in radial direction is 0.5mm. The cutter feed at 0.005 mm/tooth. In case D, rotational speed of spindle is 2600 rpm. Cutting depth in direction of spindle is 4.5mm. Cutting depth in radial direction is 0.5mm. The cutter feed at 0.005 mm/tooth. In case E, rotational speed of spindle is 3100 rpm. Cutting depth in direction of spindle is 0.5mm. Cutting depth in radial direction is 0.5mm. The cutter feed at 0.005 mm/tooth. In case F, rotational speed of spindle is 3100 rpm. Cutting depth in direction of spindle is 2mm. Cutting depth in radial direction is 0.5mm. The cutter feed at 0.005 mm/tooth.

19

4.2 Experiment result analysis Fast Fourier Transform is applied to the experimentally measured milling forces so that all the frequencies can be separated. FFT in these cases are presented in Fig. 18 (a), Fig. 18 (b), Fig. 18 (c), Fig. 18 (g), Fig. 18 (h) and Fig. 18 (i). In case A, cutting depth is below the critical stability boundaries predicted by analytical methods and numerical methods. So, chatter frequency should not be found. In Fig. 18 (a), the main frequency of cutting force is 433Hz. It is 10 times of rotational frequency of spindle. No chatter frequency is found. The stability in case A is predicted successfully by analytical methods and numerical methods. In case B, cutting depth is below the critical stability boundaries predicted by analytical methods and numerical methods. So, chatter frequency should not be found. In Fig. 18 (b), the master frequency of cutting force is 347Hz. It is 2 times of tooth passing frequency. No chatter frequency is found. The stability in case B is predicted successfully by analytical methods and numerical methods.

(a) FFT in case A

(b) FFT in case B

(c) FFT in case C stable cutting trace

(e) Surface in case C chatter wave

(d) Cases to be verified

(g) FFT in case F

(f) Surface in case D

(h) FFT in case E

(i) FFT in case D

Fig. 18. Analysis of experimental results In Fig. 18 (c), two master frequencies of cutting force are 347Hz and 433Hz. The one is 2 times of tooth passing frequency. And another one is 10 times of rotational frequency of spindle. No chatter frequency is found. This result is closer to the critical stability boundary predicted by numerical methods. There is an error in the prediction by ZOA, but when the multi-frequency

20

solution method with high order approximation ability is used, the error disappears. This still means that both numerical methods and analytical methods are valid. Surface of the workpiece after the milling operation is shown in Fig. 18 (e). The surface is smooth. In case D, cutting depth is above the critical stability boundaries predicted by analytical methods and numerical methods. So, a chatter frequency should be found. In Fig. 18 (i), a frequency at 429Hz is found. It is not several times of rotational frequency of spindle, tooth passing frequency or natural frequencies of system. Thus, it is the chatter frequency. The stability in case D is predicted successfully by analytical methods and numerical methods. Surface of the workpiece after the milling operation is shown in Fig. 18 (f). Texture due to chatter can be observed. In case E, cutting depth is below the critical stability boundaries predicted by analytical methods and numerical methods. So, chatter frequency should not be found. In Fig. 18 (h), the master frequency of cutting force is 414Hz. It is 2 times of tooth passing frequency. No chatter frequency is found. The stability in case E is predicted successfully by analytical methods and numerical methods. In case F, cutting depth is above the critical stability boundaries predicted by analytical methods and numerical methods. So, a chatter frequency should be found. In Fig. 18 (g), a frequency at 410Hz is found. It is not several times of rotational frequency of spindle, tooth passing frequency or natural frequencies of system. Thus, it is the chatter frequency. The stability in case F is predicted successfully by analytical methods and numerical methods. Through the above analysis, the experimentally measured critical stability boundary is shown in Fig. 18 (d). It can be seen that the experimentally measured critical stability boundary, critical stability boundaries predicted by analytical methods and critical stability boundaries predicted by numerical methods are consistent. Therefore, both numerical methods and analytical methods are valid. 5. Conclusion This paper proposes a numerical difference method based on Adams-Bashforth scheme for stability prediction in milling. Numerical simulations and experimental verifications were conducted. The results can be summarized as follows: 1.

When using a computer to implement the algorithm, both the numerical difference method and the multi-frequency solution method are affected by the truncation error. This effect will decrease as 𝑛𝑠 increases or 𝑓𝑛𝑐 increases or Δ𝜔 decreases or 𝑓𝑎𝑐 increases.

2.

When the workpiece is very thin, modal parameters of the system are determined by cutter as well as workpiece. It means that the stability of thin-walled workpiece milling process can hardly be improved by just increasing stiffness of the cutter. The critical stability boundaries will change with material removal.

3.

In the initial stage of thin-wall workpiece milling, the critical stability boundaries predicted by numerical methods are the same as the critical stability boundaries predicted by analytical methods, and they are all valid. As the material is removed, numerical methods and analytical methods begin to differ, but they are still similar in general.

4.

Reducing 𝑓𝑎𝑐 is an effective way to reduce computational time of multi-frequency solution method, since computational time increases with the cubic of 𝑓𝑎𝑐 . Reducing 𝑓𝑛𝑐 , 𝑓𝑛𝑑 and 𝑓𝑛𝑒 is an effective way to reduce computational time of semi-discretization method and full-discretization method, since computational time increases with the fourth power of 𝑓𝑛𝑐 . Reducing 𝑓𝑛𝑐 , 𝑓𝑛𝑑 and 𝑓𝑛𝑒 is an effective way to reduce computational time of numerical integration method and numerical difference method, since computational time increases with the cubic of 𝑓𝑛𝑐 , 𝑓𝑛𝑑 and 𝑓𝑛𝑒 .

5.

Computational time of parallel multi-frequency solution method decreases linearly with the increase of the number of threads. Computational time of parallel semi-discretization method doesn’t decrease linearly but not very bad. The following tips will help improve performance on time consumption:

1.

𝚲𝑑𝑖𝑓3𝑎 and 𝚲𝑑𝑖𝑓3𝑏 are typical sparse matrices. Therefore, when implementing the algorithm proposed in this paper, some algorithms optimized for sparse matrices can be applied.

2.

If 𝑛𝑠 is high, then 𝑓𝑛𝑐 can be small, which can greatly improve performance on time consumption.

Competing Interests The authors declare that there are no competing interests regarding the publication of this paper.

21

Acknowledgments This work was supported by the National Natural Science Foundation of China (51975112) and Fundamental Research Funds for Central Universities(N180305032) and Supported by Liao Ning Revitalization Talents Program (XLYC1807063). Appendix A Nomenclature list. 𝐶𝑕𝑖𝑝,𝑗 (𝑡):

The thickness of the chip on the No. j tool tip

𝑁:

Number of flutes on cutter.

(mm). 𝑇:

Tooth passing period (s).

𝑡:

Time /s

𝑥(𝑡):

Displacement in the X direction (mm).

𝑦(𝑡):

Displacement in the Y direction (mm).

𝐹𝑡,𝑗 :

The tangential cutting force on the No. j tool

𝐾𝑡 :

Cutting

tip (N). 𝑎𝑝 :

Cutting depth in direction of spindle (mm).

force

coefficient

in

tangential

direction (MPa). 𝐹𝑟,𝑗 :

The radial cutting force on the No. j tool tip (N).

𝐾𝑟 :

Ratio of radial cutting force coefficient to

𝐹𝑥,𝑗 :

tangential one. 𝐹𝑦,𝑗 :

The Y direction cutting force on the No. j tool

tip (N). 𝐹𝑥 :

tip (N). 𝐹𝑦 :

the resultant force on cutter in Y direction

The X direction cutting force on the No. j tool

the resultant force on cutter in X direction (N).

𝑔(φ):

(N).

A piecewise function that describes the cutting state.

𝜑𝑠𝑡 :

Start angle of immersion.

𝜑𝑒𝑥 :

Exit angle of immersion.

𝐌:

Modal mass matrix.

𝐬(𝑡):

Displacement matrix.

𝐂:

Modal damping matrix.

𝐊:

Modal stiffness matrix.

𝑚𝑥 :

Modal mass in X-direction.

𝑐𝑥 :

Modal damping in X-direction.

𝑚𝑦 :

Modal mass in Y-direction.

𝑐𝑦 :

Modal damping in Y-direction.

𝑘𝑥 :

Modal stiffness in X-direction.

𝐅(𝐬(𝑡), 𝐬(𝑡 −

Dynamic cutting force matrix.

𝑘𝑦 :

Modal stiffness in Y-direction.

𝐔(𝑡):

Directional coefficient matrix.

𝑎𝑥𝑥 (𝑡), 𝑎𝑥𝑦 (𝑡),

Directional coefficients.

𝑢(𝑡):

The first order differential of 𝑥(𝑡).

𝑣(𝑡):

The first order differential of 𝑦(𝑡).

𝜔𝑛𝑥 :

Natural frequency in the X direction (Hz).

𝜔𝑛𝑦 :

Natural frequency in the Y direction (Hz).

𝑘𝑥 :

Stiffness in the X direction (N/m).

𝑘𝑦 :

Stiffness in the Y direction (N/m).

𝜁𝑥 :

Damping ratio in the X direction.

𝜁𝑦 :

Damping ratio in the Y direction.

𝑡𝑖𝑛𝑖𝑡𝑖𝑎𝑙 :

The initial moment of the Tooth passing

𝑇), 𝑡):

𝑎𝑦𝑥 (𝑡), 𝑎𝑦𝑦 (𝑡):

period. 𝑡0 :

The initial moment of the discretization.

𝑚:

Number of intervals.

𝑡𝑚 :

The last moment of the discretization.

𝜏2 :

The span of an interval (s).

𝐏𝑑𝑖𝑓3 :

The Floquet Matrix.

𝑓𝑛𝑎 :

Number of iterations of the spindle speed.

𝑓𝑛𝑏 :

Number of iterations of axial depth of cut.

𝑓𝑛𝑐 :

Number of intervals.

𝑓𝑛𝑑 :

Number of modal parameters in the X

𝑓𝑛𝑒 :

Number of modal parameters in the Y

direction. Synonymous with 𝑁𝑚𝑥 .

direction. Synonymous with 𝑁𝑚𝑦 .

𝑓𝑎𝑎 :

Number of iterations of the spindle speed.

𝑓𝑎𝑏 :

Number of iterations of chatter frequency.

𝑓𝑎𝑐 :

The number of rows of the matrix 𝐆(𝜔𝑐 , 𝜔 𝑇 )

𝐹̅𝑥 :

Average cutting force in the X direction.

in Eq. (31) in Ref. [13]. 𝐹̅𝑦 :

Average cutting force in the Y direction.

𝐹̅𝑧 :

Average cutting force in the Z direction.

𝑓:

Feed rate in X-direction (mm / tooth).

𝐾𝑡𝑐 :

Tangential cutting force coefficients (MPa).

22

𝐾𝑡𝑒 :

Edge tangential cutting force coefficients

𝐾𝑟𝑐 :

Radial cutting force coefficients (MPa).

(MPa). 𝐾𝑟𝑒 :

Edge radial cutting force coefficients (MPa).

𝐾𝑎𝑐 :

Axial cutting force coefficients (MPa).

𝐾𝑎𝑒 :

Edge axial cutting force coefficients (MPa).

𝑛𝑠 :

Spindle speed (rpm).

𝜔𝑐 :

Chatter frequency.

𝑎𝑝,𝑚𝑎𝑥 :

Critical stability cutting depth in direction of spindle (mm).

References [1].

J.H. Shaik, J. Srinivas, Analytical prediction of chatter stability of end milling process using three-dimensional cutting force model, Journal of the Brazilian Society of Mechanical Sciences and Engineering. 39 (2017) 1633–1646. doi:10.1007/s40430-016-0567-x.

[2].

M. Wan, Y.C. Ma, J. Feng, W.H. Zhang, Study of static and dynamic ploughing mechanisms by establishing generalized model

with

static

milling

forces,

International

Journal

of

Mechanical

Sciences.

114

(2016)

120–131.

doi:10.1016/j.ijmecsci.2016.05.010. [3].

Y.C. Ma, M. Wan, W.H. Zhang, Time domain simulation of milling chatter stability, Materials Science Forum. 836–837 (2016) 94–98. doi:10.4028/www.scientific.net/MSF.836-837.94.

[4].

J. Munoa, X. Beudaert, Z. Dombovari, Y. Altintas, E. Budak, C. Brecher, G. Stepan, Chatter suppression techniques in metal cutting, CIRP Annals - Manufacturing Technology. 65 (2016) 785–808. doi:10.1016/j.cirp.2016.06.004.

[5].

Xie Q , Zhang Q , Han J, Hopf bifurcation for delay differential equation with application to machine tool chatter, Applied Mathematical Modelling, 36.8 (2012) 3803-3812. Doi:10.1016/j.apm.2011.11.011.

[6].

J. Gradišek, M. Kalveram, T. Insperger, K. Weinert, G. Stépán, E. Govekar, I. Grabec, On stability prediction for milling, International Journal of Machine Tools and Manufacture. 45 (2005) 769–781. doi:10.1016/j.ijmachtools.2004.11.015.

[7].

Y. Altintaş, E. Budak, Analytical Prediction of Stability Lobes in Milling, CIRP Annals - Manufacturing Technology. 44 (1995) 357–362. doi:10.1016/S0007-8506(07)62342-7.

[8].

S.D. Merdol, Y. Altintas, Multi frequency solution of chatter stability for low immersion milling, Journal of Manufacturing Science and Engineering, Transactions of the ASME. 126 (2004) 459–466. doi:10.1115/1.1765139.

[9].

A. Iglesias, J. Munoa, J. Ciurana, Z. Dombovari, G. Stepan, Analytical expressions for chatter analysis in milling operations with one dominant mode, Journal of Sound and Vibration. 375 (2016) 403–421. doi:10.1016/j.jsv.2016.04.015.

[10]. L. Zhu, B. Liu, H. Chen, Research on chatter stability in milling and parameter optimization based on process damping, JVC/Journal of Vibration and Control. 24 (2018) 2642–2655. doi:10.1177/1077546317692159. [11]. E. Graham, M. Mehrpouya, S.S. Park, Robust prediction of chatter stability in milling based on the analytical chatter stability, Journal of Manufacturing Processes. 15 (2013) 508–517. doi:10.1016/j.jmapro.2013.08.005. [12]. T. Insperger, G. Stépán, Updated semi-discretization method for periodic delay-differential equations with discrete delay, International Journal for Numerical Methods in Engineering. 61 (2004) 117–141. doi:10.1002/nme.1061. [13]. Y. Altintas, G. Stepan, D. Merdol, Z. Dombovari, Chatter stability of milling in frequency and discrete time domain, CIRP Journal of Manufacturing Science and Technology. 1 (2008) 35–44. doi:10.1016/j.cirpj.2008.06.003. [14]. Ding, Y. , Zhu, L. M. , Zhang, X. J. , & Ding, H. . A full-discretization method for prediction of milling stability. International Journal of Machine Tools & Manufacture, 50(5) (2010) 502-509. doi:10.1016/j.ijmachtools.2010.01.003. [15]. Y. Ding, L. Zhu, X. Zhang, H. Ding, Numerical integration method for prediction of milling stability, Journal of Manufacturing Science and Engineering, Transactions of the ASME. 133 (2011). doi:10.1115/1.4004136. [16]. F.I. Compeán, D. Olvera, F.J. Campa, L.N. López De Lacalle, A. Elías-Zúñiga, C.A. Rodríguez, Characterization and stability analysis of a multivariable milling tool by the enhanced multistage homotopy perturbation method, International Journal of Machine Tools and Manufacture. 57 (2012) 27–33. doi:10.1016/j.ijmachtools.2012.01.010. [17]. G. Urbikain, L.N. López De Lacalle, F.J. Campa, A. Fernández, A. Elías, Stability prediction in straight turning of a flexible workpiece by collocation method, International Journal of Machine Tools and Manufacture. 54–55 (2012) 73–81.

23

doi:10.1016/j.ijmachtools.2011.11.008. [18]. G. Urbikain, F.J. Campa, J.J. Zulaika, L.N. López De Lacalle, M.A. Alonso, V. Collado, Preventing chatter vibrations in heavy-duty turning operations in large horizontal lathes, Journal of Sound and Vibration. 340 (2015) 317–330. doi:10.1016/j.jsv.2014.12.002. [19]. C. YUE, H. GAO, X. LIU, S.Y. LIANG, L. WANG, A review of chatter vibration research in milling, Chinese Journal of Aeronautics. 32 (2019) 215–242. doi:10.1016/j.cja.2018.11.007. [20]. P. V. Bayly, B.P. Mann, T.L. Schmitz, D.A. Peters, G. Stepan, T. Insperger, Effects of radial immersion and cutting direction on chatter instability in end-milling, ASME International Mechanical Engineering Congress and Exposition, Proceedings. (2002) 351–363. doi:10.1115/IMECE2002-39116. [21]. T. Insperger, Full-discretization and semi-discretization for milling stability prediction: Some comments, International Journal of Machine Tools and Manufacture. 50 (2010) 658–662. doi:10.1016/j.ijmachtools.2010.03.010. [22]. K. Ahmadi, F. Ismail, Stability lobes in milling including process damping and utilizing Multi-Frequency and Semi-Discretization Methods, International Journal of Machine Tools and Manufacture. 54–55 (2012) 46–54. doi:10.1016/j.ijmachtools.2011.11.007. [23]. B. Denkena, C. Schmidt, Experimental investigation and simulation of machining thin-walled workpieces, Production Engineering. 1 (2007) 343–350. doi:10.1007/s11740-007-0017-9. [24]. D. Biermann, P. Kersting, T. Surmann, A general approach to simulating workpiece vibrations during five-axis milling of turbine blades, CIRP Annals - Manufacturing Technology. 59 (2010) 125–128. doi:10.1016/j.cirp.2010.03.057. [25]. U. Bravo, O. Altuzarra, L.N. López De Lacalle, J.A. Sánchez, F.J. Campa, Stability limits of milling considering the flexibility of the workpiece and the machine, International Journal of Machine Tools and Manufacture. 45 (2005) 1669–1680. doi:10.1016/j.ijmachtools.2005.03.004. [26]. E. Budak, L.T. Tunç, S. Alan, H.N. Özgüven, Prediction of workpiece dynamics and its effects on chatter stability in milling, CIRP Annals - Manufacturing Technology. 61 (2012) 339–342. doi:10.1016/j.cirp.2012.03.144. [27]. V. Thevenot, L. Arnaud, G. Dessein, G. Cazenave-Larroche, Integration of dynamic behaviour variations in the stability lobes method: 3D lobes construction and application to thin-walled structure milling, International Journal of Advanced Manufacturing Technology. 27 (2006) 638–644. doi:10.1007/s00170-004-2241-1. [28]. L. Zhu, B. Liu, X. Wang, Z. Xu, Research on Cutting Force of Turn-Milling Based on Thin-Walled Blade, Advances in Materials Science and Engineering. 2016 (2016). doi:10.1155/2016/2638261. [29]. L. Zhu, Z. Yang, Z. Li, Investigation of mechanics and machinability of titanium alloy thin-walled parts by CBN grinding head,

International

Journal

of

Advanced

Manufacturing

Technology.

100

(2019)

2537–2555.

doi:10.1007/s00170-018-2795-y. [30]. L. Arnaud, O. Gonzalo, S. Seguy, H. Jauregi, G. Peigné, Simulation of low rigidity part machining applied to thin-walled structures,

International

Journal

of

Advanced

Manufacturing

Technology.

54

(2011)

479–488.

doi:10.1007/s00170-010-2976-9. [31]. A. Iglesias, J. Munoa, C. Ramírez, J. Ciurana, Z. Dombovari, FRF Estimation through Sweep Milling Force Excitation (SMFE), Procedia CIRP. 46 (2016) 504–507. doi:10.1016/j.procir.2016.04.019. [32]. D. Olvera, A. Elías-Zúñiga, O. Martínez-Romero, L.N. López de Lacalle, H. Martínez-Alfaro, H.R. Siller, M.W. Pineda, Improved predictions of the stability lobes for milling cutting operations of thin-wall components by considering ultra-miniature accelerometer mass effects, International Journal of Advanced Manufacturing Technology. 86 (2016) 2139–2146. doi:10.1007/s00170-015-8287-4. [33]. J.J. Zulaika, F.J. Campa, L.N. Lopez De Lacalle, An integrated processmachine approach for designing productive and lightweight milling machines, International Journal of Machine Tools and Manufacture. 51 (2011) 591–604. doi:10.1016/j.ijmachtools.2011.04.003. [34]. J. Munoa, A. Iglesias, A. Olarra, Z. Dombovari, M. Zatarain, G. Stepan, Design of self-tuneable mass damper for modular fixturing systems, CIRP Annals - Manufacturing Technology. 65 (2016) 389–392. doi:10.1016/j.cirp.2016.04.112.

24

[35]. J. Munoa, X. Beudaert, K. Erkorkmaz, A. Iglesias, A. Barrios, M. Zatarain, Active suppression of structural chatter vibrations using machine drives and accelerometers, CIRP Annals - Manufacturing Technology. 64 (2015) 385–388. doi:10.1016/j.cirp.2015.04.106. [36]. J. Munoa, I. Mancisidor, N. Loix, L.G. Uriarte, R. Barcena, M. Zatarain, Chatter suppression in ram type travelling column milling machines using a biaxial inertial actuator, CIRP Annals - Manufacturing Technology. 62 (2013) 407–410. doi:10.1016/j.cirp.2013.03.143. [37]. P. Kersting, D. Biermann, Modeling workpiece dynamics using sets of decoupled oscillator models, Machining Science and Technology. 16 (2012) 564–579. doi:10.1080/10910344.2012.731948. [38]. V. Thevenot, L. Arnaud, G. Dessein, G. Cazenave-Larroche, Influence of material removal on the dynamic behavior of thin-walled

structures

in

peripheral

milling,

Machining

Science

and

Technology.

10

(2006)

275–287.

doi:10.1080/10910340600902082. [39]. Altintas, Y. Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design, second ed., Cambridge University Press, New York, 2012. [40]. Kurdi, M. H., Haftka, R. T., Schmitz, T. L., Mann, B. P., A robust semi-analytical method for calculating the response sensitivity

of

a

time

delay

system,

Journal

https://doi.org/10.1115/1.2981093.

25

of

Vibration

and

Acoustics.

130

(2008)

064504.