A novel Chebyshev-wavelet-based approach for accurate and fast prediction of milling stability

Journal Pre-proof A novel Chebyshev-wavelet-based approach for accurate and fast prediction of milling stability Chengjin Qin, Jianfeng Tao, Haotian Shi, Dengyu Xiao, Bingchu Li, Chengliang Liu PII:

S0141-6359(19)30557-4

DOI:

https://doi.org/10.1016/j.precisioneng.2019.11.016

Reference:

PRE 7067

To appear in:

Precision Engineering

Received Date: 27 July 2019 Revised Date:

16 November 2019

Accepted Date: 21 November 2019

Please cite this article as: Qin C, Tao J, Shi H, Xiao D, Li B, Liu C, A novel Chebyshev-wavelet-based approach for accurate and fast prediction of milling stability, Precision Engineering (2019), doi: https:// doi.org/10.1016/j.precisioneng.2019.11.016. 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.

A novel Chebyshev-wavelet-based approach for accurate and fast prediction of milling stability Chengjin Qin1, Jianfeng Tao1*, Haotian Shi1, Dengyu Xiao1, Bingchu Li1, 2, Chengliang Liu1 1

State Key Laboratory of Mechanical System and Vibration, School of Mechanical Engineering,

Shanghai Jiao Tong University, Shanghai 200240, China 2

School of Mechanical Engineering, University of Shanghai for Science and Technology, Shanghai

200093, China

*Corresponding author, E-mail: [email protected]

Abstract： ：Currently, semi-analytical stability analysis methods for milling processes focus on improving prediction accuracy and simultaneously reducing computing time. This paper presents a Chebyshev-wavelet-based method for improved milling stability prediction. When including regenerative effect, the milling dynamics model can be concluded as periodic delay differential equations, and is re-presented as state equation forms via matrix transformation. After divide the period of the coefficient matrix into two subintervals, the forced vibration time interval is mapped equivalently to the definition interval of the second kind Chebyshev wavelets. Thereafter, the explicit Chebyshev–Gauss–Lobatto points are utilized for discretization. To construct the Floquet transition matrix, the state term is approximated by finite series second kind Chebyshev wavelets, while its derivative is acquired with a simple and explicit operational matrix of derivative. Finally, the milling stability can be semi-analytically predicted using Floquet theory. The effectiveness and superiority of the presented approach are verified by two benchmark milling models and comparisons with the representative existing methods. The results demonstrate that the presented approach is highly accurate, fast and easy to implement. Meanwhile, it is shown that the presented approach achieves high stability prediction accuracy and efficiency for both large and low radial-immersion milling operations. Key words ： Milling stability prediction; periodic delay differential equations; second kind Chebyshev wavelets; Chebyshev-wavelet-based approach; Floquet theory

1 Introduction Due to the flexibility of the machining system, an unexpected unstable instability, known as regenerative chatter, is always unavoidable during the machining process [1]. It arises from a self-excitation mechanism, and is commonly concluded as delay differential equations (DDEs)

1

with time-periodic coefficients [2, 3]. Chatter has various negative impacts on machining processes, including poor surface quality, accelerated tool wear, harsh noise and reduced machining efficiency [4, 5]. For these reasons, researches on the prediction [6, 7], monitoring [8, 9], and suppression [10] of milling chatter have attracted great attention from academia and industry. In order to avoid the occurrence of chatter, selecting the appropriate machining parameter combinations based on the stability analysis is of vital importance. Therefore, the need for accurate and efficient milling stability analysis is ever increasing. Considering that periodic delay

differential equations have infinite dimensional state space, researchers have introduced the numerical simulation methods for milling stability prediction over the early decades. Davies et al. [11] presented a time-domain solution-based stability prediction technique for low radial immersion milling operations. To predict forces, surface finish, and chatter stability under very small radial immersion condition, Campomanes and Altintas [12] proposed an improved time domain method, in which non-linear effects that are difficult to model could be accurately accounted for. Based on a numerical discretization scheme along the variable cutting edge profile, Urbikain et al. [13] proposed a comprehensive dynamic force model for milling operations with a barrel shape, in which stability maps with much more information could be obtained. Li and Liu [14] developed a predictive time-domain model for stability analysis of end milling, which could also applied to complicated milling processes with complicated tool geometry or time-varying cutting conditions. Although the numerical simulation methods achieve good adaptability for complicated milling stability problems, they require large number of iterative operations to complete the calculation, leading to extremely long computation time. Fortunately, the analytical and semi-analytical methods make up for this deficiency. By transforming the original DDEs into the frequency-domain describe, Altintas and Budak [15] aquried the analytical expression for the critical depth of cut, known as the zeroth order approximation method. This method is extremely advantageous in terms of computational efficiency. However, since only the zero-order term of Fourier expansion is retained, it cannot be competent for small radial-immersion milling problems. In order to improve the prediction accuracy, Merdol and Altintas [16] proposed the so-called multi-frequency method. From another point of view, scholars have proposed the Floquet-theory-based semi-analytical methods for milling stability analysis. For instance, Insperger and co-workers [17, 18] proposed the robust semi-discretization method (0st SDM) for stability prediction of periodic delay-differential equations. To further obtain higher accuracy, higher semi-discretization methods were presented and discussed by Insperger et al. [19]. It was found that the first-order semi-discretization method (1st SDM) was enough when approximating the time-periodic coefficients by piecewise constant functions. The SDM has been experimentally validated and widely used to predict the stability of various milling processes, thus it is commonly utilized as the benchmark method for reference. Bayly et al. [20] introduced the so-called temporal finite element analysis method to overcome stability analysis problems for interrupted cutting. Butcher et al. [21] presented the Chebyshev

2

collocation method for milling analysis, in which the optimal stable immersion levels were also analyzed. Over the past ten years, scholars have made great efforts to improve both the calculation accuracy and speed. By separating the original coefficient matrix into a constant matrix and a periodic matrix, Ding et al. [22] presented the full-discretization method (FDM) for stability prediction of milling. Compared with SDMs, the FDM has been shown to improve computational speed without sacrificing accuracy. Later, to increase the calculation speed, Liu et al. [23] proposed an efficient FDM, which has a relatively compact algorithm structure. The complete discretization scheme (CDS) was recommended by Li and co-workers [24], in which the exponential matrix is eliminated by numerical differentiation for higher efficiency. Ding et al. [25] introduced the numerical integration method (NIM) for stability analysis of milling. Niu et al. [26] presented the generalized Runge–Kutta method (GRKM). It achieves high computational accuracy, but the relatively complicated structure leads to a decrease in calculation speed. To improve the computational accuracy and speed, Qin et al. [27, 28] proposed the so-called holistic-discretization methods (HDMs). On the basis of Runge-Kutta methods, Li et al. [29] proposed an enhanced complete discretization method (RKCDM) for the prediction of milling stability. Olvera and co-workers [30, 31] suggested the enhanced multistage homotopy perturbation method (EMHPM) for the prediction of milling stability lobes. The second-order semi-discretization method (2nd SDM) was proposed by Jiang and co-workers [32], in which the Newton interpolation polynomials were applied. On the other hand, with the help of higher interpolation methods, Quo et al. [33] and Ozoegwu et al. [34] presented higher-order FDMs. Unfortunately, these methods sacrifice computational speed while improving prediction accuracy. Recently, two updated versions of FDMs (UFDM and 3rd UFDM) were recommended by Tang et al. [35] and Yan et al. [36], in which both the state term and the delay term were approximated with higher-order interpolation. Qin et al. [37] proposed an Adams-Moulton-based method for stability prediction of milling operations. To improve calculation accuracy and reduce computational time, Qin et al. [38] further developed the Adams-Simpson-based method by using the predictor-corrector scheme. From the perspective of numerical differentiation, zhang et al. [39] proposed the Simpson-based method (SBM). Lu et al. [40] recommended the spline-based method for stability analysis of milling processes. It was shown that the spline-based method has higher accuracy and efficiency for milling with variable pitch cutters at low radial immersion cases compared with [22] and [25]. By using the finite difference method and extrapolation method, Zhang et al. [41] proposed the numerical differentiation method for stability prediction of two degrees of freedom (DOF) high speed milling. To obtain higher computation efficiency, Dong et al. [42] proposed a reconstructed version of SDM by using the Shannon standard orthogonal basis, which was further extended to milling with variable spindle speed [43]. Recently, wavelets that possess many useful properties have found their way into many different engineering and science fields [44]. As very well-localized functions, wavelets not only allows accurate representation of functions, but also establishes links with fast numerical algorithms.

3

Besides the signal analysis and image processing, there is a growing research interest in using wavelets for approximating solution of differential equations with great computational complexity [45]. Among various wavelet families, the Chebyshev wavelets have demonstrated the potential to obtain both high rapidity of convergence and easy implement-ability [46, 47]. To improve the stability prediction accuracy and reduce the computational cost, this work proposes a novel Chebyshev-wavelet-based method (CWM) for milling stability analysis. The highlight of the proposed CWM is that it approximates the state term by the second kind Chebyshev wavelets, and obtains the corresponding derivative term by a simple and explicit operational matrix of derivative. By converting and discretizing periodic DDE into a system of algebraic equations, the Floquet transition matrix can be finally constructed and utilized for stability prediction of milling. The rest of this paper is organized as follows. Section 2 gives a brief description for the second kind Chebyshev wavelets. The Chebyshev-wavelet-based method is developed step by step in Section 3. Section 4 validates computational accuracy, speed and applicability of the proposed method by benchmark milling examples and making comparisons with the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM. Finally, the conclusion is drawn in Section 5.

2 The second kind Chebyshev wavelets As we all know, wavelets are a family of functions consisting of dilation and translation of a single function ψ(t) known as the mother wavelet. The second kind Chebyshev wavelets ψnm(t)=ψ(k,n,m,t) are defined on the interval [0, 1] and have four arguments, namely

 2( k +3)/ 2 U m (2k +1 t − 2n − 1)  ψ nm (t ) =  π 0 

n n +1 ≤t < k k 2 2 otherwise

（1）

where n=0,1,…,2k–1 and m=0,1,…,M. k is a non-negative integer. Um(t) is the second kind Chebyshev polynomials of degree m, and is defined by

U m (t ) =

sin [ ( m + 1)θ ] sin(θ )

, t = cos(θ )

（2）

On the interval [–1, 1], the second kind Chebyshev polynomials are orthogonal with respect to the

weight function, and can be alternatively generated by using the following recurrence relation

U m+1 (t ) = 2tU m (t ) − U m−1 (t ), m = 1, 2,L

（3）

where U0(t)=1, and U1(t)=2t. To intuitively recognize and better understand the second kind Chebyshev wavelet, Fig. 1 illustrates the curve of its value over time with k=0, n=0, and M=5. Utilizing the orthonormality, a function g(t) defined over the interval [0,1] can be approximated by finite series second kind Chebyshev wavelets as

4

2 k −1 M

g (t )= ∑ ∑ cnmψ nm (t )

（4）

n =0 m = 0

where cnm denote the unknown coefficients to be determined.

Fig.1 The second kind Chebyshev wavelets with k=0, n=0, and M=5

3 Mathematical model and algorithm description For milling operations, mathematical model of regenerative chatter can be formulated as the following delay differential equations with periodic coefficients

&&(t ) + Cq& (t ) + Kq (t ) = − a p K c (t ) [q (t ) − q (t − T ) ] Mq

（5）

where M, C and K denote the modal matrices of the milling system. Meanwhile, q(t) represents the tool modal displacement vector, and ap denotes the depth of cut. The coefficient matrix Kc(t) is time-periodic, satisfying Kc(t)= Kc(t+T). Here, Τ equals to the tooth-passing period, i.e., T=60/(ΩN), where N is the number of tool teeth, and Ω represents the spindle speed (rpm). To construct the Floquet transition matrix, the proposed method converts the periodic DDE into a system of algebraic equations, and is implemented by the following steps. With the aid of matrix transformation, the dynamics model of milling operations is re-expressed as state equation forms. Then it divides the period of the coefficient matrix into two subintervals based on its value. Moreover, the forced vibration time period is converted to the definition interval of the Chebyshev wavelets, and the explicit Chebyshev–Gauss–Lobatto points are adopted for discretization. On this basis, it approximates the state term by the Chebyshev wavelets, and obtains the corresponding derivative term by the explicit operational matrix of derivative. Finally, the Floquet transition matrix can be obtained for prediction of the milling stability

5

By applying the matrix transformation p(t)=Mq(t)+Cq(t) and x(t)=[q(t), p(t)]T , the milling dynamics equation can be re-expressed into the following state equation forms

x& (t ) = ( Ac + A(t ) ) x(t ) − A(t )x(t − T )

（6）

with

 −M −1C / 2 Ac =  −1  CM C / 4 − K

0 0  M −1  A ( ) = ， t    −CM −1 / 2   − a p K c (t ) 0 

（7）

It is seen that the time-periodic coefficient matrix Kc(t) determines the vibration form of the system described by Eq. (6). When the matrix Kc(t) is a zero matrix, it experiences free vibration process. On the contrary, if matrix Kc(t) contains any non-zero elements, it experiences forced vibration process. Consequently, we divide the period of the coefficient matrix Τ into two sub-intervals, i.e., the free vibration time interval and the forced vibration time interval. Denote the forced vibration time period by To=T−Tr, where Tr is the free vibration time period. Introducing transformation t=Toη+Tr, η∈[0, 1], the forced vibration interval [Tr, T] can be converted to the interval [0, 1], and Eq. (6) can be equivalently re-expressed as

x& (η ) = To ( Ac + A(η ) ) x(η ) − A(η )x(η − T / To )

（8）

For discretization, the explicit Chebyshev–Gauss–Lobatto points are utilized, namely

ti = 0.5 1− cos ( iπ / md )  , i = 0,L, md

（9）

where md=2k(M+1)–1. During the free vibration time interval [0, Tr], Eq. (8) has an analytical solution. Specifically, the state term x(t0) can be analytically deduced as

x (t0 ) = e A cTr x(tmd − T / To )

（10）

At the discrete time points ti, the state equation Eq. (8) should satisfy

 ( A c + A(t0 ) ) x(t0 ) − A(t0 )x(t0 − T / To )   x& (t0 )     &  A c + A(t1 ) ) x(t1 ) − A(t1 )x(t1 − T / To )  ( x ( t )  1   =T   M  o M     &   x ( t ) A + A ( t ) x ( t ) − A ( t ) x ( t − T / T ) m d   md md md md o   c

(

（11）

)

Subsequently, it approximates the state term by the second kind Chebyshev wavelets, and obtains the corresponding derivative term by the explicit operational matrix of derivative. At discrete time points ti, the scalar function x(ti) is approximated by 2k −1 M

x (ti )= ∑ ∑ bnmψ nm (ti ) n =0 m =0

6

（12）

For simplicity, let bi=bnm, ψi(ti)=ψnm(ti), i=nM+m, then the unknown coefficients bi can be

obtained as follows

 b0  b1 b=   M   bmd

  ψ 0 (t0 ) ψ 1 (t0 ) L ψ md (t0 )     ψ ψ ψ ( t ) ( t ) L ( t ) m 0 1 1 1 1   d =    M M M M     ψ 0 (tmd ) ψ 1 (tmd ) L ψ md (tmd ) 

−1

 x(t0 )   x (t0 )       x(t1 )  =Ψ −1  x (t1 )   M   M       x(tmd )   x(tmd ) 

（13）

On the basis of the second kind Chebyshev wavelet operational matrix of derivative, one can obtain the derivative of the state function x(ti) by

( x&(t ), x&(t ),L, x&(t )) = b D( ψ(t ), ψ(t ),L, ψ(t )) =b DΨ T

0

1

md

T

0

1

T

md

（14）

Here the operational matrix of derivative D is a (md+1)×(md+1) matrix and can be deduce as

H 0 L 0    0 H L 0  D= M M O M    0 0 L H

（15）

where H is a (M+1)×(M+1) matrix and its (i, j)th element Hij is given by

2k +2 j Hij =  0

i = 2,L,(M +1); j = 1,L,(i −1) and (i + j) odd otherwise

（16）

By combining Eq. (11) and Eq. (12), and employing the Kronecker product ⊗, one can obtain a vector form of the derivative x& (ti ) by

 x& (t0 )   x(t0 )   x(t0 )   &       x(t1 )  = ΨDT Ψ −1 ⊗ I  x(t1 )  = G  x(t1 )  ) s× s  M   M   M  (       & x ( t ) x ( t ) x ( t ) m m m  d   d   d 

（17）

in which s represents the dimension of vector x(t). Combining the Eqs. (10), (11) and (17), a discrete map can be obtained by

 x(t0 − T / To )   x(t0 )      x(t1 )  x(t1 − T / To )    E =F    M  M      x(tmd )   x(tmd − T / To )  with

7

（18）

0 0  0 A(t1 ) Q = − M M  0 0 E= 0  0 R = Q+  M  0

  L 0  O M   L A(tmd )  L

0

1 Gu + Q To 0 L e AcTr   0 L 0  M O M   0 L 0 

（19）

（20）

（21）

It should be pointed out that the constant matrix Gu is obtained by modifying the first n rows of the matrix G as [In×n, 0n×n , …, 0n×n]. Finally, the Floquet transition matrix Ф with the proposed method can be constructed by

Φ = E − 1F

（22）

According to the Floquet theory, when the spectral radius κ(Ф) of the Floquet transition matrix is larger than one, the milling process is unstable.

4 Numerical results and validation To verify the effectiveness and superiority of the proposed method, comparisons among the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the proposed CWM are conducted. Herein, the benchmark milling operations in Ref. [18] are utilized for algorithm validation, which have been experimentally validated. For equal and complete comparison, all the algorithms employ the same program structure in MATLAB 2017a and run on the same desktop computer (Intel Core i3-6500 3.3 GHz, 4.0 GB of DDR3L RAM, Windows 10 OS).

4.1 Benchmark milling operations The single and two DOF milling operations in Ref. [18] have been verified by experiments and are commonly used as the benchmark examples for algorithm validation. By applying the matrix transformation described in Section 3, the single DOF milling dynamics equation in Ref. [18] can be re-expressed into the state equation forms Eq. (6), with

−ζωn  Ac =  2 2 2  mt ζ ωn − mt ωn

0 0 1/ mt    , A (t ) =  − a h (t ) 0  −ζωn   p xx 

8

（23）

Similarly, for the two DOF milling, the corresponding coefficient matrices can be obtained as −ζωn   0 Ac =  2  (ζ −1)ωn 2 mt  0 

0

1/ mt 0

−ζωn

−ζωn

0 (ζ − 1)ωn mt 2

2

0

0 0 0     0 0 1/ mt  , A(t ) =   −a p hxx (t ) −a p hxy (t ) 0    −ζωn   −a p hyx (t ) −a p hyy (t )

0 0  0 0  （24） 0 0  0 0 

In Eq. (23) and Eq. (24), mt, ζ and ωn are the modal parameters of the system. For the two-DOF milling operations, the modal parameters are assumed to be equal at the x and y directions. Herein, hxx(t), hxy(t), hyx(t) and hyy(t) are the elements of the coefficient matrix Kc(t), and can be deduced as N  h ( t ) = [ K t cos(φk (t )) + K n sin(φk (t ))]σ (φk (t )) sin(φk (t )) ∑  xx k =1  N   hxy (t ) = ∑ [ K t cos(φk (t )) + K n sin(φk (t ))]σ (φk (t )) cos(φk (t ))  k =1  N  h (t ) = [ − K sin(φ (t )) + K cos(φ (t ))]σ (φ (t )) sin(φ (t )) ∑ t k n k k k  yx k =1  N  h (t ) = [ − K sin(φ (t )) + K cos(φ (t ))]σ (φ (t )) cos(φ (t )) ∑ t k n k k k  yy k =1

（25）

where Kt and Kn are the cutting force coefficients. The angular position of the kth tooth φk(t) is calculated by φk=(2πΩ/60)t+2π (k–1)/N. The window function, σ(φk(t)) is utilized to determine whether the kth tool teeth is cutting the part, i.e., σ(φk(t))=1 when φst <φk(t)< φex. Herein, φst and φex are the start and exit immersion angles. For down-milling operations, φst=arcos(2a/D–1), and φex=π; for up-milling operations, φst=0, and φex=arcos(1–2a/D); where a/D is the radial immersion ratio. The system parameters utilized for validation are the same with those in Ref. [18]: down-milling, N=2, mt=0.03993 kg, ζ=0.011, ωn=922×2π rad/s, Kt=600 MPa and Kn=200 MPa. Fig.2 illustrates the real part and the imaginary part of frequency response function (FRF).

Fig.2 The frequency response function (FRF) of the benchmark milling operations

9

4.2 Convergence rate verification To evaluate and verify the convergence rate of the proposed method, comparisons with the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM are conducted. As an important evaluation index of algorithm accuracy, the convergence rate is illustrated by expressing the difference between the approximate and the exact value as a function of the time intervals md. Insperger et al. [19] has proved that the local discretization error of the 1st SDM is O(τ3), where τ represents discrete step. According to [36], as an enhanced version of FDMs, the 3rd UFDM achieved higher convergence rate than NIM and its approximation order is higher than O(τ3). For the GRKM, Niu et al. [26] found that its approximation order is larger than O(τ4). In addition, it has been shown in [38] that the ASM obtained higher convergence rate than AMM and its local discretization error is higher than O(τ4). Following a similar way in [19], it was found that the local discretization error of the EMPHM in [30] is also O(τ3). In the proposed method, the state term is approximated by the second kind Chebyshev wavelets, which have spectral convergence accuracy. The discrete step size is not uniform, and the time intervals md equals to number of the Chebyshev wavelets. Hence, the proposed CWM is expected to be spectral convergent, and obtain

higher convergence rate than the other methods. Fig. 3 illustrates the convergence rates of the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM with radial immersion ratio a/D=1.0. The reference spectral radius |κ0| is calculated by the CWM with md =700, and can be considered as the exact one. And the spectral radius calculated by the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM with different time intervals md is denoted as |κ|. For verification, different spindle speeds and depths of cut are used: Ω=5000 rpm, ap=1.2 mm and 2.4 mm; Ω=8000 rpm, ap=1.3 mm and 2.6 mm. In Fig. 3, the longitudinal coordinates are represented in the logarithmic form for ease of comparison. It is seen that as the time intervals md increases, the spectral radius calculated by the CWM converges to the exact spectral radius much faster than those calculated by the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM. Hence, the proposed method achieves higher convergence rate and computational accuracy than the other methods.

10

Fig.3 Comparisons of convergence rates among the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM with the radial immersion ratio a/D=1.0 Fig. 4 presents the convergence rates of the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the proposed method with the radial immersion ratio a/D=0.05. The system parameters are the same as in the previous case. The reference spectral radius |κ0| is also calculated by the CWM with md =700. Herein the parameter combinations of spindle speeds and depths of cut are selected as follows: Ω=6000 rpm, ap=1.5 mm and 3.5 mm; Ω=10000 rpm, ap=1.8 mm and 3.8 mm. As shown in Fig. 4, the convergence rates of the proposed method for the radial immersion ratio a/D=0.05 is higher than those of the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM. Consequently, the comparison results show that the proposed method obtains higher convergence rate than the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM for both large and low radial immersion ratios.

11

Fig.4 Comparisons of convergence rates among the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM with the radial immersion ratio a/D=0.05

4.3 Milling stability analysis and comparison Herein, to validate the effectiveness of the proposed method and examine the calculation speed, the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM are utilized to predict the milling stability for the benchmark milling operations in Ref. [18]. The stability lobes diagrams for the single-DOF milling model with radial immersion ratios a/D=1.0 and the computation time of these methods are presented in Fig. 5. The time intervals md of the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM are set as 24 and 32, respectively. In Fig. 5, the reference stability boundaries demoded by the red line are calculated by the CWM with md=300. The stability lobes diagrams are constructed over a 200×200 sized grid, and the domain of machining parameters are set as follows: the depths of cut ap∈ [0 mm, 5 mm], and the spindle speed Ω∈[5 krpm, 25 krpm]. It is seen that the accuracy of the stability lobes diagrams calculated by the proposed method is better than those calculated by the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM. Hence, the proposed method obtains higher computational accuracy than the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM. Furthermore, the result shows that the proposed CWM achieves higher calculation speed than the other methods. Compared with the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM, the CWM can save the computation time by almost 97%, 83.7~88.9%, 70%~78%, 65.8~72.6% and 76%~77.8%, respectively.

12

13

Fig. 5 Stability lobes diagrams for the single-DOF milling model with full radial immersion ratio a/D=1.0 and the computation time of the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM The stability lobes diagrams for the single-DOF milling model with half radial immersion ratio a/D=0.5 and the computation time of these methods are presented in Fig. 6. The time intervals md of the 1st SDM, the GRKM and the CWM are set as 12 and 16, respectively. Also, the reference stability boundaries are calculated by the CWM with md=300, and the stability lobes diagrams are constructed over a 200×200 sized grid. The domain of machining parameters are now set as follows: the depths of cut ap∈ [0 mm, 3.5 mm], and the spindle speed Ω∈[5 krpm, 25 krpm]. As shown in Fig. 6, the proposed approach obtains both higher prediction accuracy and calculation speed. The accuracy of the stability lobes diagrams with the proposed approach is better than those calculated by the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM, which demonstrates that the presented CWM achieves a convergence rate than the other methods. Simultaneously, the result show that under the same calculation parameters, the computation time of the CWM can be reduced by approximately 98% compared with the 1st SDM and by 89.3%~90.5% compared with the 3rd UFDM. In addition, the CWM can save the computation time by 77%~81% compared with the GRKM, by 72.9%~75.7% compared with the ASM and by 77% compared with the EMHPM.

14

15

Fig. 6 Stability lobes diagrams for the single-DOF milling model with full radial immersion ratio a/D=0.5 and the computation time of the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM To further verify the effectiveness of the proposed method, the stability lobes diagrams for the two-DOF milling model and the computation time of the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM are presented in Fig. 7. Now, the immersion ratio a/D is set as 1.0 and the stability lobes diagrams are constructed over a 200×100 sized grid. The reference stability boundaries are demoded by the red line and are calculated by the CWM with md=300. Besides, the number of time intervals md is set as 12 and 16, respectively. According to Fig. 7, it is seen that the accuracy of the stability lobes diagrams calculated by the proposed method is better than those calculated by the other methods. Hence, the proposed method obtains higher computational accuracy than the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM. Furthermore, the result shows that the proposed CWM achieves higher calculation speed than the other methods. Compared with the 1st SDM, the computation time of the proposed method can be saved by almost 90.6%~91.7%. Compared with the 3rd UFDM, the computational time of the CWM can be reduced by approximately 63.3%~65.5%. Besides, compared with the GRKM, ASM and the EMPHM, the proposed method can reduce the computation time by 45.6%~46.7%, 39.8%~42.0% and 47.8%~54.2%, respectively. In summary, the results that the proposed method achieves both higher computational accuracy and calculation speed for stability prediction of both single and two DOF milling operations.

16

17

Fig.7 Stability lobes diagrams for the two-DOF milling model with full radial immersion ratio a/D=0.5 and the computation time of the 1st SDM, the 3rd UFDM, the GRKM, the ASM, the EMHPM and the CWM

5 Conclusion In this work, a novel Chebyshev-wavelet-based method is developed for accurately and efficiently predicting the milling stability. To obtain the Floquet matrix for stability analysis, it discretizes the forced vibration time interval into the explicit Chebyshev–Gauss–Lobatto points, and approximates the state term of the state equation by the second kind Chebyshev wavelets. Simultaneously, a simple and explicit operational matrix of derivative is presented to obtain the derivative term of the state equation. It is a new approach of using the wavelets to semi-analytically predict the milling stability governed by DDEs. The analysis and comparison of convergence rate show that the proposed CWM is spectral convergent and its convergence rate is much higher than the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM. Moreover, when making stability prediction for the benchmark milling operations, it is seen that the proposed approach achieves higher calculation speed than the other methods. When compared with the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM, the proposed method can reduce the computation time by almost 97%~98%, 83.7~90.5%, 70%~81%, 65.8%~75.7% and 76%~77.8% for the single DOF milling operation, and by approximately 90.6%~91.7%, 63.3%~65.5%, 45.6%~46.7%, 39.8%~42.0% and 47.8%~54.2% for the two DOF milling operation. The main contribution of this study is to propose an alternative for stability prediction of milling processes with a single discrete time delay. In the future, efforts will be made to extend the proposed method to more complex milling processes and complex tools (e.g., serrated end mills, variable helix/pitch end mills).

18

Acknowledgements This work was partially supported by the National Key R&D Program of China (Grant No. 2018YFB1702503), the Science and Technology Planning Project of Guangdong Province (Grant No. 2017B090914002) and the China Postdoctoral Science Foundation (Grant No. 2019M661496）.

Conflict of interest The authors declare that they have no conflict of interest.

References 1.

Dai Y, Li H, Xing X, Hao B. Prediction of chatter stability for milling process using precise integration method. Precis Eng 2018;52:152–157

2.

Yue C, Gao H, Liu X, Liang SY, Wang L. A review of chatter vibration research in milling. Chin J Aeronaut 2019;32(2):215–242

3.

Altintas Y, Stepan G, Merdol D, Dombovari Z. Chatter stability of milling in frequency and discrete time domain. CIRP J Manuf Sci Technol 2008;1(1): 35–44

4.

Quintana G, Ciurana J Stepan. Chatter in machining processes: a review. Int J Mach Tools Manuf 2011;51(5): 363–376

5.

Tao JF, Qin CJ, Xiao DY, Shi HT, Ling X, Li BC, Liu CL. Timely chatter identification for robotic drilling using a local maximum synchrosqueezing-based method. J Intell Manuf, 2019; Doi: 10.1007/s10845-019-01509-5

6.

Dong XF, Qiu ZZ.Stability analysis in milling process based on updated numerical integration

method. Mech Syst Signal Process, 2019; Doi: 10.1016/j.ymssp.2019.106435 7.

Huang T, Zhang X, Ding H. A novel approach with smallest transition matrix for milling stability prediction. Nonlinear Dynamics 2017;90(1):95–104

8.

Tao JF, Qin CJ, Liu CL. A synchroextracting-based method for early chatter identification of robotic drilling process. Int J Adv Manuf Technol 2019;100:273–285

9.

Tao JF, Qin CJ, Xiao DY, Shi HT, Liu CL. A pre-generated matrix-based method for real-time robotic drilling chatter monitoring. Chin J Aeronaut 2019; Doi: 10.1016/j.cja.2019.09.001

10. Munoa J, Beudaert X, Dombovari Z, Altintas Y, Budak E, Brecher C, Stepan G. Chatter suppression techniques in metal cutting.CIRP Ann Manuf Techn 2016;65(2):785–808 11. Davies MA, Pratt JR, Dutterer B, Burns TJ. Stability prediction for low radial immersion

19

milling. J Manuf Sci E-T ASME 2002;124: 217–225 12. Campomanes ML, Altintas Y. An improved time domain simulation for dynamic milling at small radial immersions. J Manuf Sci E-T ASME 2003;125: 416–422 13. Urbikain G, Olvera D, López de Lacalle LN. Stability contour maps with barrel cutters considering the tool orientation. Int J Adv Manuf Technol 2017;89(9–12): 2491–2501 14. Li ZQ, Liu Q. Solution and analysis of chatter stability for end milling in the time-domain. Chin J Aeronaut 2008;21: 169–178 15. Altintas Y, Budak E. Analytical prediction of stability lobes in milling. CIRP Ann 1995;44(1): 357–362 16. Merdol SD, Altintas Y. Multi frequency solution of chatter stability for low immersion milling. J Manuf Sci Eng 2004;126(3): 459–466 17. Insperger T, Stepan G. Semi-discretization method for delayed systems. Int J Numer Methods Biomed Eng 2002;55(5): 503–518 18. Insperger T, Stepan G. Updated semi-discretization method for periodic delay-differential equations with discrete delay. Int J Numer Methods Biomed Eng 2004;61(1): 117–141 19. Insperger T, Stepan G, Turi J. On the higher-order semidiscretizations for periodic delayed systems. J Sound Vib 2008;313(1): 334–341 20. Bayly PV, Halley JE, Mann BP, Davies MA. Stability of interrupted cutting by temporal finite element analysis. Trans ASME J Manuf Sci Eng 2003;125: 220–225 21. Butcher EA, Bobrenkov OA, Bueler E, Nindujarla P. Analysis of milling stability by the Chebyshev collocation method: Algorithm and optimal stable immersion levels. J Comput Nonlin Dyn 2009;4: 031003 22. Ding Y, Zhu LM, Zhang XJ, Ding H. A full-discretization method for prediction of milling stability. Int J Mach Tools Manuf 2010;50(5): 502–509 23. Liu Y, Zhang DH, Wu BH. An efficient full-discretization method for prediction of milling stability. Int J Mach Tools Manuf 2012; 63:44–48 24. Li MZ, Zhang GJ, Huang Y. Complete discretization scheme for milling stability prediction. Nonlinear Dyn 2013;71: 187–199 25. Ding Y, Zhu LM, Zhang XJ, Ding H. Numerical integration method for prediction of milling stability. J Manuf Sci Eng 2011;133(3): 031005 26. Niu JB, Ding Y, Zhu LM, Ding H. Runge–Kutta methods for a semi-analytical prediction of milling stability. Nonlinear Dyn 2014;76(1):289–304

20

27. Qin CJ, Tao JF, Liu CL. A predictor-corrector-based holistic-discretization method for accurate and efficient milling stability analysis. Int J Adv Manuf Technol 2018;96:2043–2054 28. Qin CJ, Tao JF, Liu CL. A novel stability prediction method for milling operations using the holistic-interpolation scheme. Proc IMechE Part C: J Mechanical Engineering Science

2019;233(13):4463–4475 29. Li ZQ, Yang ZK, Peng YR, Zhu F, Ming XZ. Prediction of chatter stability for milling process using Runge-Kutta-based complete discretization method. Int J Adv Manuf Technol 2016;86(1–4):943–952 30. Olvera D, Elías-Zúñiga A, Martínez-Alfaro H, López de Lacalle LN, Rodríguez CA, Campa FJ. Determination of the stability lobes in milling operations based on homotopy and simulated annealing techniques. Mechatronics 2014;24: 177–185 31. Compeán FI, Olvera D, Campa FJ, López de Lacalle LN, Elías-Zúñiga A, Rodríguez CA. Characterization and stability analysis of a multi variable milling tool by the enhanced multistage homotopy perturbation method. Int J Mach Tools Manuf 2012;57: 27–33 32. Jiang SL, Sun YW, Yuan XL, Liu WR. A second-order semi-discretization method for the efficient and accurate stability prediction of milling process. Int J Adv Manuf Tech 2017;92(1–4):583–595 33. Quo Q, Sun YW, Jiang Y. On the accurate calculation of milling stability limits using third-order full-discretization method. Int J Mach Tools Manuf 2012;62: 61–66 34. Ozoegwu CG, Omenyi SN, Ofochebe SM. Hyper-third order full-discretization methods in milling stability prediction. Int J Mach Tools Manuf 2015;92: 1–9 35. Tang X, Peng F, Yan R, Gong Y, Li Y, Jiang L. Accurate and efficient prediction of milling stability

with

updated

fulldiscretization

method.

Int

J

Adv

Manuf

Technol

2017;88(9–12):2357–2368 36. Yan Z, Wang X, Liu Z, Wang D, Jiao L, Ji Y. Third-order updated full-discretization method for milling stability prediction. Int J Adv Manuf Technol 2017;92(5–8):2299–2309 37. Qin CJ, Tao JF, Li L, Liu CL. An Adams-Moulton-based method for stability prediction of milling processes. Int J Adv Manuf Technol 2017;89(9–12): 3049–3058 38. Qin CJ, Tao JF, Liu CL. Stability analysis for milling operations using an Adams-Simpson-based method. Int J Adv Manuf Technol 2017;92(1–4): 969–979 39. Zhang Z, Li HG, Meng G, Liu C. A novel approach for the prediction of the milling stability based on the Simpson method. Int J Mach Tools Manuf 2015;99: 43–47

21

40. Lu YA, Ding Y, Peng ZK, Chen ZZC, Zhu LM. A spline-based method for stability analysis of milling processes. Int J Adv Manuf Technol 2016;89:2571–2586 41. Zhang XJ, Xiong CH, Ding Y, Ding H. Prediction of chatter stability in high speed milling using the numerical differentiation method. Int J Adv Manuf Technol 2016;89:2535–2544 42. Dong XF, Zhang WM, Deng S. The reconstruction of a semi-discretization method for milling stability prediction based on Shannon standard orthogonal basis. Int J Adv Manuf Technol 2016; 85(5–8):1501–1511 43. Dong XF, Zhang WM. Chatter suppression analysis in milling process with variable spindle speed based on the reconstructed semi-discretization method. Int J Adv Manuf Technol 2019; Doi: 10.1007/s00170-019-04363-0 44. Hariharan G, Kannan K. Review of wavelet methods for the solution of reaction-diffusion problems in science and engineering. Applied Mathematical Modelling 2016;38(3):799–813 45. Hariharan G, Rajaraman R, Sathiyaseelan D. Wavelet based spectral algorithm for nonlinear dynamical systems arising in ship dynamics. Ocean Engineering 2016;126:321–328 46. Sahu K, Saha S. Chebyshev wavelet method for numerical solutions of integro-differential form of Lane-Emden type differential equations. International Journal of Wavelets, Multiresolution and Information Processing 2017;15 (2):1750015 47. Yi M, Ma Y, Wang L. An efficient method based on the second kind Chebyshev wavelets for solving variable-order fractional convection diffusion equations. International Journal of Computer Mathematics 2017;95(10):1973–1991

22

Highlights This work focuses on the accurate and fast stability prediction in milling operations by developing a Chebyshev-wavelet-based method. After divide the period of the coefficient matrix into two subintervals, the forced vibration time interval is mapped equivalently to the definition interval of the second kind Chebyshev wavelets. It approximates the state term by the second kind Chebyshev wavelets, and obtains the corresponding derivative term by a simple and explicit operational matrix of derivative. It is shown that the proposed method is spectral convergent and its convergence rate is much higher than the 1st SDM, the 3rd UFDM, the GRKM, the ASM and the EMHPM. Moreover, it demonstrates that the proposed approach achieves higher speed than the other methods for stability prediction of the benchmark milling models.

Conflict of interest The authors declare that they have no conflict of interest.