# Stability Analysis for Milling Process with Variable Pitch and Variable Helix Tools by High-Order Full-Discretization Methods

Authors: Yang Zhang; Kenan Liu; Wuyun Zhao; Wei Zhang; Fei Dai
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Chatter is one of the significant limitations in the milling process, which may cause poor surface quality, reduced productivity, and accelerated tool wear. Variable pitch and variable helix tools can be used to suppress regenerative chatter. This study extends the high-order full-discretization methods (FDMs) to predict the stability of milling with variable pitch and variable helix tools. The time-periodic delay-differential equation (DDE) with multiple delays is used to model the milling process using variable pitch and variable helix tools. Then, the DDE with multiple delays is reexpressed by the state-space equation. Meanwhile, the spindle rotational period is divided into many small-time intervals, and the state space equation is integrated on the small-time interval. Then, the high-order interpolation polynomials are used to approximate the state term, and the weights related to the time delay are employed to approximate the time-delay term. The second-order, third-order, and fourth-order extended FDMs (2nd EFDM, 3rd EFDM, and 4th EFDM) are compared with the benchmark in terms of the rate of convergence. It is found that the 2nd EFDM, 3rd EFDM, and 4th EFDM converge faster than the benchmark method. The difference between the curves obtained by different EFDMs and the reference curve is very small. There is no need to extend hypersecond FDMs to analyze the stability of milling with variable pitch and variable helix tools.


---

## Body

## 1. Introduction

High precision products machined by using high-speed milling technology have an increasing demand in various industries. Chatter is one of the significant limitations in the milling process, which may cause poor surface quality, reduced productivity, and accelerated tool wear. Therefore, a great deal of effort has been made to develop chatter avoidance methods. The stability lobe diagram (SLD) is an effective way to predict the stability in the milling process.In the milling process, regenerative chatter can be modeled by time-periodic DDE. The mapping relations of the dynamic responses between two consecutive periods can be obtained by solving the DDE. Then, the stability lobe diagrams can be determined according to the obtained mapping relations. Until now, different methods for calculating stability lobe diagrams have been reported. Altintas and Budak [1] developed a famous and efficient zeroth-order approximation (ZOA) method that determines the stability lobe diagram in the frequency domain. In this method, the time-varying cutting force coefficients are approximated by their Fourier series components. However, this method cannot accurately predict the stability lobes at low radial immersion milling operations. Hence, Merdol and Altintas [2] presented a so-called multifrequency method. This method can be used to accurately predict the milling stability under low radial immersion conditions in the frequency domain. After that, many time-domain methods have been developed to predict the stability in milling. Bayly et al. [3] proposed a temporal finite element analysis method, Butcher et al. [4] presented a Chebyshev collocation method, and Insperger and Stépán reported the updated semidiscretization method (SDM) [5] and first-order semidiscretization method (1st SDM) [6]. The abovementioned time-domain methods can be used to predict the stability of milling accurately; however, their computational efficiency should be further improved. To improve computational efficiency, Ding et al. [7] presented a full-discretization method (FDM) based on the direct integration scheme. Then, on the framework of FDM, Ding et al. [8] introduced the second-order FDM (2nd FDM), Guo et al. [9] proposed the third-order FDM (3rd FDM), and Ozoegwu reported the least squares approximation methods [10] and hyper-third-order full-discretization methods [11]. Besides, the numerical integration method [12] and differential quadrature method [13] developed by Ding and his coworkers, the Runge–Kutta-based methods proposed by Niu et al. [14], the improved precise integration method proposed by Li et al. [15], the Simpson-based method presented by Zhang et al. [16], numerical differentiation method reported by Zhang et al. [17], and so on are proposed for the stability prediction in the milling process. More recently, Olvera et al. [18] presented the stability analysis for a single degree of freedom down-milling operation in a thin-walled workpiece by using the enhanced multistage homotopy perturbation (EMHP) method.The above works mainly studied the chatter avoidance methods by using stability lobe diagrams. However, chatter suppression is also a critical issue for the study of stability in milling. It is well known that variable pitch and variable helix (nonuniform) milling tools can be used to suppress regenerative chatter by disturbing the regeneration mechanism. Many studies on the stability analysis in milling with variable pitch and variable helix tools have been reported. Based on the ZOA method, Altintas et al. [19] presented an analytical model for stability prediction in milling with variable pitch cutters. Based on this model, a method for selecting optimal pitch angles is also presented. Sims et al. [20] utilized the semidiscretization method, time-averaged semidiscretization method, and temporal finite element method to predict the chatter stability for variable pitch and variable helix milling tools. Sellmeier and Denkena [21] proposed a method based on Ackermann’s approach for investigating the stability of an unequally pitched end mill. Compeán et al. [22] suggested the enhanced multistage homotopy perturbation method (EMHPM) to predict the stability of a multivariable milling tool. Wan et al. [23] presented a unified method for predicting the stability lobes of the milling process with multiple delays. Jin et al. [24] proposed a frequency-domain solution for efficient stability prediction of variable helix cutters milling. After that, Jin et al. [25] extended the semidiscretization method to predict the stability of the variable pitch and variable helix milling cutter. Guo et al. [26] proposed a new variable interpolation method to predict the stability of milling with multidelays combining the effect of the cutter’s variable-helix angle’s effect. Sims [27] described a new formulation for predicting the chatter stability of variable helix tools based on Laplace formulation. Wang et al. [28] proposed an improved semidiscretization algorithm (ISDA) for modeling and simulation with variable pitch and variable helix milling tools. Niu et al. [29] studied the mechanics and multiregenerative stability of variable pitch and variable helix milling tools considering runout. Huang et al. [30] developed the robust active chatter control for milling processes with variable pitch cutters whose dynamics are governed by multidelay nonlinear differential equations. Stepan et al. [31] presented the achievable upper and lower capability bounds by introducing so-called stabilizability diagrams of a hypothetical variable pitch milling cutter that is tuned continuously along the stability boundaries. Huang et al. [32] presented a stability analysis method for milling system with variable pitch cutters under variable speed. Regarding the mechanics model for nonuniform milling tools, Chen et al. [33] proposed a unified analytical cutting force model based on a predictive machining theory for variable helix end mill considering cutter runout. Otto et al. [34] derived a dynamic process model for milling with nonuniform pitch and variable helix tools. In this model, the nonlinear cutting force behavior and the effect of runout are included. Guo et al. [35] presented a compensated-chip thickness-based cutting force model for nonuniform helix tools in the five-axis milling process. Cai et al. [36] proposed an integrated process-machine model based on the computer graphics method for simulating the milling process of a variable pitch cutter. Additionally, the stability lobe diagrams are usually considered in the design process of the nonuniform milling tools. Regarding the milling tool design with consideration of stability, Budak [37, 38] presented a design method for optimizing the pitch angles of the milling tool. Yusoff and Sims [39] proposed a semidiscretization method combined with differential evolution to predict milling stability as well as optimize variable helix end milling tools. Comak and Budak [40] presented a practical but accurate design method for the selection of the best variation combination to maximize chatter-free material removal rate without using time-consuming computer simulations. Mei et al. [41] proposed an analytical method for designing milling cutters with alternating variable pitches. Guo et al. [42] proposed an effective optimization method for variable helix cutter by introducing an index called “suppression factor” to measure stability quantitatively. Moreover, serrated milling tools can also be employed to suppress chatter in the milling process. Merdol and Altintas [43, 44] analyzed the mechanics and dynamics of milling with serrated cutters, where the periodic serrations ground on the helical flutes are modeled by fitting a cubic spline on the design profile. Dombovari et al. [45, 46] presented the stability model of serrated cutters in analytical semidiscrete time domain. Hosseini et al. [47] proposed mechanistic modeling for cutting with serrated end mills using a parametric representation approach. Tehranizadeh and Budak [48] presented the methods for designing optimum serration shapes. In their methods, the effects of serration shapes on the mechanics of milling are investigated. Guo et al. [49] presented a mechanistic cutting force model of the serrated end mill to predict cutting forces. Bari et al. [50] presented an improved chip thickness model for serrated end mills that account for the actual trochoidal path traced by the tooth. Tehranizadeh et al. [51] investigated the effects of different waveforms on the mechanics of process. In their work, the effects of phase shift direction and local cutting angles on milling forces are also experimentally verified. Pelayo et al. [52] studied the mechanics of the milling system with serrated end mills using force and surface topography models. In their study, a stationary milling force model is developed including the main geometric parameters of the serrated profile, and a surface topography model is developed to predict the resulting machined surfaces.In the literature [11], Ozoegwu et al. pointed out that the accuracy of stability result rises with the order of full-discretization method peaking at fourth order and then declined at the fifth order. That is, in the hyper-fourth-order full-discretization methods, the accuracy of the stability results may be affected by the Runge phenomenon. Although it is indicated from the literature [11] that the fourth-order FDM is the best for a uniform pitch and uniform helix tool, whether it is the best for a variable pitch and variable helix cutter is also worthy of study. In this study, the high-order FDMs which are proposed for the analysis of milling with uniform tools are extended to analyze the stability of milling with variable pitch and variable helix tools, and the performances of extended FDMs are also evaluated.This paper is organized as follows: in Section2, the mathematical model of milling dynamics is introduced; in Section 3, the high-order FDMs are extended to predict the stability of milling with variable pitch and variable helix tools; in Section 4, the rates of convergence of the extended FDMs are analyzed, and the stability lobe diagrams for different conditions are obtained; and conclusions are drawn in Section 5.


## 2. Mathematical Model of Milling Dynamics

In the process of modeling the dynamics of the milling process with variable pitch and variable helix tools, the axial depth of cutap is equally discretized into L axial layers with the thickness of Δz; thus, Δz=ap/L.


### 2.1. Angular Position

The angular position of thejth tooth on the lth (1 ≤ l ≤ L) axial layer φl,jt is represented as(1)φl,jt=2πΩ60t−apl−1/2tanβjRL,ifj=1,2πΩ60t+∑j=2jPj−apl−1/2tanβjRL,if1<j≤N,where Ω is the spindle speed in rpm, ap is the axial depth of cut, βj is the helix angle of the jth tooth, R is the radius of the tool shank envelope, Pj is the pitch angle between the jth tooth and (j-1)th tooth, and N is the number of cutter tooth.


### 2.2. The Dynamic Chip Thickness

The dynamic chip thickness caused by regenerative effect,dl,jt, can be given by(2)dl,jt=xt−xt−τl,jsinφl,jt+yt−yt−τl,jcosφl,jt,where xt and yt are the displacements of the current tooth in the X and Y directions, respectively; xt−τl,j and yt−τl,j are the displacements of the previous tooth in the X and Y directions, respectively; and the time delay τl,j depends on the separation angle between the jth tooth and (j−1)th tooth on the lth layer, Δφl,j. Thus, the time delay τl,j can be given by(3)τl,j=TΔφl,j2π,where the period T is equal to the spindle rotational period, namely, T=60/Ω. The separation angle Δφl,j can be expressed as follows [24, 29]:(4)Δφl,j=Pj+l−1/2tanβj−tanβNLR,ifj=1,Pj+l−1/2tanβj−tanβj−1LR,if1<j≤N.


### 2.3. Milling Dynamics

The governing equation of milling dynamics for a two-degree-of-freedom milling system can be written as follows [5]:(5)x¨ty¨t+2ζxωnx002ζyωnyx˙ty˙t+ωnx200ωny2xtyt=∑j=1N∑l=1L−hxxtmx−hxytmx−hyxtmy−hyytmyxt−xt−τl,jyt−yt−τl,j,where ζx, ωnx, and mx are the damping ratio, angular natural frequency, and modal mass of the milling system in the X direction, respectively; and ζy, ωny, and my are the damping ratio, angular natural frequency, and modal mass of the milling system in the Y direction, respectively. The specific cutting force coefficients hxxt, hxyt, hyxt, and hyyt are given as follows:(6)hxxt=gφl,jtsinφl,jtKtcosφl,jt+Knsinφl,jtΔz,hxyt=gφl,jtcosφl,jtKtcosφl,jt+Knsinφl,jtΔz,hyxt=gφl,jtsinφl,jt−Ktsinφl,jt+Kncosφl,jtΔz,hyyt=gφl,jtcosφl,jt−Ktsinφl,jt+Kncosφl,jtΔz,where Kt and Kn are the tangential and normal cutting coefficients, respectively. The window function gφl,jt is used to determine whether the tooth is in or out of cut [5].


## 2.1. Angular Position

The angular position of thejth tooth on the lth (1 ≤ l ≤ L) axial layer φl,jt is represented as(1)φl,jt=2πΩ60t−apl−1/2tanβjRL,ifj=1,2πΩ60t+∑j=2jPj−apl−1/2tanβjRL,if1<j≤N,where Ω is the spindle speed in rpm, ap is the axial depth of cut, βj is the helix angle of the jth tooth, R is the radius of the tool shank envelope, Pj is the pitch angle between the jth tooth and (j-1)th tooth, and N is the number of cutter tooth.


## 2.2. The Dynamic Chip Thickness

The dynamic chip thickness caused by regenerative effect,dl,jt, can be given by(2)dl,jt=xt−xt−τl,jsinφl,jt+yt−yt−τl,jcosφl,jt,where xt and yt are the displacements of the current tooth in the X and Y directions, respectively; xt−τl,j and yt−τl,j are the displacements of the previous tooth in the X and Y directions, respectively; and the time delay τl,j depends on the separation angle between the jth tooth and (j−1)th tooth on the lth layer, Δφl,j. Thus, the time delay τl,j can be given by(3)τl,j=TΔφl,j2π,where the period T is equal to the spindle rotational period, namely, T=60/Ω. The separation angle Δφl,j can be expressed as follows [24, 29]:(4)Δφl,j=Pj+l−1/2tanβj−tanβNLR,ifj=1,Pj+l−1/2tanβj−tanβj−1LR,if1<j≤N.


## 2.3. Milling Dynamics

The governing equation of milling dynamics for a two-degree-of-freedom milling system can be written as follows [5]:(5)x¨ty¨t+2ζxωnx002ζyωnyx˙ty˙t+ωnx200ωny2xtyt=∑j=1N∑l=1L−hxxtmx−hxytmx−hyxtmy−hyytmyxt−xt−τl,jyt−yt−τl,j,where ζx, ωnx, and mx are the damping ratio, angular natural frequency, and modal mass of the milling system in the X direction, respectively; and ζy, ωny, and my are the damping ratio, angular natural frequency, and modal mass of the milling system in the Y direction, respectively. The specific cutting force coefficients hxxt, hxyt, hyxt, and hyyt are given as follows:(6)hxxt=gφl,jtsinφl,jtKtcosφl,jt+Knsinφl,jtΔz,hxyt=gφl,jtcosφl,jtKtcosφl,jt+Knsinφl,jtΔz,hyxt=gφl,jtsinφl,jt−Ktsinφl,jt+Kncosφl,jtΔz,hyyt=gφl,jtcosφl,jt−Ktsinφl,jt+Kncosφl,jtΔz,where Kt and Kn are the tangential and normal cutting coefficients, respectively. The window function gφl,jt is used to determine whether the tooth is in or out of cut [5].


## 3. The Proposed Methods

Letut=xtytx˙ty˙tT, and equation (5) can be represented in state-space form as follows:(7)u˙t=Aut+∑j=1N∑l=1LBl,jtut−ut−τl,j,where(8)A=00100001−ωnx20−2ζxωnx00−ωny20−2ζyωny,Bl,jt=00000000−hxxmx−hxymx00−hyxmy−hyymy00.To solve equation (7) numerically, the period T is equally divided into k small-time intervals firstly; thus, T=kΔt, where k is an integer and Δt is the length of the small-time interval. The integer kl,j related to the time delay τl,j can be expressed as(9)kl,j=intτl,j+Δt/2Δt,where int (∗) is the function that rounds positive numbers towards zero.Then, equation (7) is integrated on the ith small-time interval ti,ti+1 resulting in(10)ut=eAt−tiuti+∫titeAt−ξ∑j=1N∑l=1LBl,jξuξ−Bl,jξuξ−τl,jdξ.Equation (10) can be equivalently written as follows [7]:(11)ui+1=eAΔtui+∫titi+1eAti+1-ξ∑j=1N∑l=1LBl,jξuξ−Bl,jξuξ−τl,jdξ.In equation (11), ui+1 and ui denote uti+1 and uti, respectively.In this work, the existing 2nd FDM, 3rd FDM, and 4th FDM are all extended to analyze the stability of milling with variable pitch and variable helix tools, and the corresponding extended methods are denoted as 2nd EFDM, 3rd EFDM, and 4th EFDM. In the 2nd EFDM, 3rd EFDM, and 4th EFDM, the second-, third-, and fourth-order interpolation polynomials are used to approximate the state termuξ, respectively. The periodic coefficient term Bl,jξ is approximated by first-order interpolation polynomial, and the time-delay term uξ−τl,j is approximated by using the weights related to time delay. Since the derivation processes of the 2nd EFDM, 3rd EFDM, and 4th EFDM are similar, to avoid repetitive work, only the detailed derivation process of the 4th EFDM is given in this work.In the 4th EFDM, the state termuξ and the periodic coefficient term Bl,jξ can be, respectively, approximated by the fourth-order and first-order interpolation polynomials using the Lagrange interpolation scheme. In the calculation process, the interpolation nodes ti−3,ti−2,ti−1,ti,ti+1 and the corresponding node values ui−3, ui−2, ui−1, ui, ui+1 are employed to approximate the state term. Meanwhile, the interpolation nodes ti, ti+1 and the node values Bi, Bi+1 are employed to approximate the periodic coefficient term. Then, the state term uξ and the periodic coefficient term Bl,jξ can be expressed as follows:(12)uξ=ξ424Δt4+ξ312Δt3−ξ224Δt2−ξ12Δtui−3+−ξ46Δt4−ξ32Δt3+ξ26Δt2+ξ2Δtui−2+ξ44Δt4+ξ3Δt3+ξ24Δt2−3ξ2Δtui−1+−ξ46Δt4−5ξ36Δt3−5ξ26Δt2+5ξ6Δt+1ui+ξ424Δt4+ξ34Δt3+11ξ224Δt2+ξ4Δtui+1,(13)Bl,jξ=Δt−ξΔtBil,j+ξΔtBi+1l,j.The time-delay termuξ−τl,j is approximated by using the weights ωal,j and ωbl,j, which can be written as(14)uξ−τl,j=ωbl,j−ξΔtui−kl,j+1+ωal,j+ξΔtui−kl,j,where the weights ωal,j and ωbl,j are calculated as follows:(15)ωal,j=kl,jΔt+Δt/2−τl,jΔt,(16)ωbl,j=τl,j+Δt/2−kl,jΔtΔt=1−ωal,j.Substituting equations (12)–(16) into equation (11) yields(17)ui+1=PiF0+∑j=1N∑l=1LG17Bil,j+G18Bi+1l,jui+∑j=1N∑l=1LG15Bil,j+G16Bi+1l,jui−1+∑j=1N∑l=1LG13Bil,j+G14Bi+1l,jui−2+∑j=1N∑l=1LG11Bil,j+G12Bi+1l,jui−3−∑j=1N∑l=1LG23Bil,j+G24Bi+1l,jui−kl,j+1−∑j=1N∑l=1LG21Bil,j+G22Bi+1l,jui−kl,j,where(18)G11=−F624Δt5−F524Δt4+F48Δt3+F324Δt2−F212Δt,G12=F624Δt5+F512Δt4−F424Δt3−F312Δt2,G13=F66Δt5+F53Δt4−2F43Δt3−F33Δt2+F22Δt,G14=−F66Δt5−F52Δt4+F46Δt3+F32Δt2,G15=−F64Δt5−3F54Δt4+3F44Δt3+7F34Δt2−3F22Δt,G16=F64Δt5+F5Δt4+F44Δt3−3F32Δt2,G17=F66Δt5+2F53Δt4−5F33Δt3−F26Δt2+F1,G18=−F66Δt5−5F56Δt4−5F46Δt3+5F36Δt2+F2Δt,G19=−F624Δt5−5F524Δt4−5F424Δt3+5F324Δt2+F24Δt,G110=F624Δt5+F54Δt4+11F424Δt3+F34Δt2,G21=F3Δt2+ωal,j−2ΔtF2+1−ωal,jF1,G22=−F3Δt2+1−ωal,jΔtF2,G23=−F3Δt2+1−ωal,jΔtF2+ωal,jF1,G24=F3Δt2+ωal,jΔtF2,Pi=I−∑j=1N∑l=1LG19Bil,j+G20Bi+1l,j−1,where I denotes the identity matrix and the matrixes F1–F6 can be expressed as follows:(19)F0=eAΔt,F1=F0−IA−1,F2=F1−ΔtIA−1,F3=2F2−Δt2IA−1,F4=3F3−Δt3IA−1,F5=4F4−Δt4IA−1,F6=5F5−Δt5IA−1.According to equation (17), the following local discrete map can be obtained:(20)ui+1uiui−1⋮ui+1−kmax=Diuiui−1ui−2⋮ui−kmax,where kmax denotes the maximum value of kl,j, and matrix Di can be expressed as follows:(21)Di=M11iM12iM13iM14i⋯000I000⋯0000I00⋯000⋮⋮⋮⋮⋮⋮⋮⋮000000I0+∑j=1N∑l=1L0⋯M1,kl,jiM1,kl,j+1i⋯00⋯00000⋯0000⋮⋮⋮⋮⋮⋮0⋯00⋯0,where(22)M11i=PiF0+∑j=1N∑l=1LG17Bil,j+G18Bi+1l,j,M12i=Pi∑j=1N∑l=1LG15Bil,j+G16Bi+1l,j,M13i=Pi∑j=1N∑l=1LG13Bil,j+G14Bi+1l,j,M14i=Pi∑j=1N∑l=1LG11Bil,j+G12Bi+1l,j,M1,kl,ji=−PiG23Bil,j+G24Bi+1l,j,M1,kl,j+1i=−PiG21Bil,j+G22Bi+1l,j.The positions of the matricesM1,kl,ji and M1,kl,j+1i in equation (21) depend on the value of kl,j related to the time delay τl,j. As for the two-degree-of-freedom milling system, the matrix M1,kl,ji is located in columns (4kl,j–3) to (4kl,j) of the matrix Di, and the matrix M1,kl,j+1i is located in columns (4kl,j + 1) to (4kl,j + 4) of the matrix Di. Additionally, the dimension of the matrix Di depends on the value of kmax.Then, the Floquet state transition matrixψ is obtained as(23)ψ=Dk−1Dk−2…D0.According to Floquet theory, the stability of the milling system can be determined using the criterion expressed by equation (24). The Matlab code for the stability analysis by using the 4th EFDM is in Appendix:(24)maxeigψ<1,stable,=1,critical stable,>1,unstable.


## 4. Stability Analysis

The two-degree-of-freedom milling system is employed for verifying the accuracy and efficiency of the extended FDMs. The rate of convergence and computational time are analyzed for the study of accuracy and efficiency, respectively. The parameters are chosen as the same as those used in the literature [28]. The detailed parameters are as follows: the number of cutter teeth is N = 4, the diameter of the cutter is 19.05 mm, the modal masses are mx = 1.4986 kg and my = 1.199 kg, the damping ratios are ζx = 0.055801 and ζy = 0.025004, the angular natural frequencies are ωnx = 3540.89 rad/s and ωny = 3243.44 rad/s, and the tangential and the normal cutting force coefficients are Kt = 697 MPa and Kn = 256 MPa, down-milling.


### 4.1. Rate of Convergence

The rate of convergence is used to estimate the local error between the absolute value of the maximum-magnitude eigenvalue of the system |μ(k)| and the reference eigenvalue μ0. In this section, the parameter combinations (Ω = 6000 rpm, ap = 0.5 mm) and (Ω = 7000 rpm, ap = 2 mm) are employed to analyze the rates of convergence of the extended FDMs. The radial immersion ratio is set as a/D = 1, the pitch angle combination of the cutter is P = [70°, 110°, 70°, 110°], and the corresponding helix angle combination is β = [30°, 35°, 30°, 35°]. The axial depth of cut is discretized into 20 axial layers, that is, L = 20. To determine the reference eigenvalues for the adopted parameter combinations, the convergences of the eigenvalues obtained by different methods are studied, as shown in Figure 1.Figure 1
The convergences of the eigenvalues obtained by ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM: (a) Ω = 6000 rpm andap = 0.5 mm; (b) Ω = 7000 rpm and ap = 2.0 mm.
(a)(b)It is seen from Figure1 that the ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM all converge to the same reference value when the parameter k reaches 800. In this study, the reference eigenvalue μ0 is determined by using the 2nd EFDM with k = 800. It should be noted that the reference eigenvalue μ0 can also be determined by using the ISDA, 3rd EFDM, and 4th EFDM since all methods converge to the same reference value. After calculation, the reference eigenvalues for the parameter combinations (Ω = 6000 rpm and ap = 0.5 mm) and (Ω = 7000 rpm and ap = 2 mm) are 0.483935 and 0.955073, respectively.In practice, the calculation burden will increase severely when the parameterk is chosen as a very large value. Therefore, the parameter k is usually smaller than that used to determine the value of μ0. With the aim of estimating the rate of convergence of the extended FDMs, the rate of convergence of the improved semidiscretization algorithm (ISDA) [28] is taken as the benchmark. Figure 2 illustrates the rates of convergence of the ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM with different values of parameter k.Figure 2
Rates of convergence of the ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM: (a)ap = 0.5 mm, Ω = 6000 rpm, and |μ0| = 0.483935; (b) ap = 2 mm, Ω = 7000 rpm, and |μ0| = 0.955073.
(a)(b)From Figure2, it is found that the difference between the approximated eigenvalues and the ideals ones calculated by the ISDA has more significant fluctuation than that calculated by the 2nd EFDM, 3rd EFDM, and 4th EFDM. Besides, the eigenvalues obtained by the 2nd EFDM, 3rd EFDM, and 4th EFDM converge to a constant when the parameter k is equal to 300, while the eigenvalues obtained by the ISDA do not converge to a constant under the same condition. It is also seen from Figure 2 that the 2nd EFDM, 3rd EFDM, and 4th EFDM converge faster than the ISDA. The analysis results indicate that the 2nd EFDM, 3rd EFDM, and 4th EFDM are more robust and accurate than the benchmark method.According to Figure2, it is also indicated that there is no significant difference between the 2nd EFDM, 3rd EFDM, and 4th EFDM in terms of the rate of convergence. It is difficult to evaluate the accuracy of different extended FDMs by using the rates of convergence of limited parameter combinations. In the following section, the accuracy of different extended FDMs is further studied by taking the stability lobe diagrams into consideration.


### 4.2. Stability Lobe Diagrams

To verify the effectiveness of the 2nd EFDM, 3rd EFDM, and 4th EFDM, the time-averaged SDM suggested by Sims et al. [20] and the extended SDM suggested by Jin et al. [25] are taken for comparison. The parameters derived from the literature [20] are used for analysis. The stability lobe diagrams obtained by the time-averaged SDM, the extended SDM, the 2nd EFDM, 3rd EFDM, and 4th EFDM are described in Figure 3.Figure 3
Stability lobe diagrams obtained by the time-averaged SDM [20], the extended SDM [25], and the 2nd EFDM, 3rd EFDM, and 4th EFDM.As shown in Figure3, the parameter combinations A (Ω = 8500 rpm and ap = 5 mm) and B (Ω = 8500 rpm and ap = 7 mm) are considered. In the literature [21], the parameter combinations A and B are proved to be a stable case and an unstable case through time-domain simulations, respectively. However, in the stability lobe diagram obtained by the time-averaged SDM, the parameter combination A corresponds to an unstable situation. That is, the time-averaged SDM cannot predict the stability correctly when the spindle speed ranges from 8000 rpm to 9000 rpm. It is seen from Figure 3 that the stabilities of parameter combinations A and B are correctly predicted in the stability lobe diagrams obtained by the 2nd EFDM, 3rd EFDM, and 4th EFDM. Additionally, the stability lobe diagrams obtained by the extended FDMs agree well with that obtained by the extended SDM [25]. Therefore, the 2nd EFDM, 3rd EFDM, and 4th EFDM are proved to be effective methods.The computational time consumed by the 2nd EFDM, 3rd EFDM, and 4th EFDM is compared with that consumed by enhanced multistage homotopy perturbation method (EMHPM) [22] and to illustrate the efficiencies of the extended FDMs. In the calculation process, the pitch angle combination of the cutter is P = [70°, 110°, 70°, 110°], and the corresponding helix angle combination is β = [30°, 35°, 30°, 35°]. The parameter plane, which includes the parameter combinations of spindle speed and axial depth of cut, is divided into a 100 × 100 size grid. The spindle speed ranges from 2 × 103 rpm to 10 × 103 rpm, and the axial depth of cut ranges from 0 mm to 20 mm. The parameter k is chosen as 80, and the axial depth of cut is discretized into 10 axial layers, that is, L = 10. The referenced stability lobe diagrams denoted by the red line are obtained by the 2nd EFDM with k = 300. The stability lobe diagrams obtained by the EMHPM, 2nd EFDM, 3rd EFDM, and 4th EFDM as well as the corresponding computational time under full immersion condition (a/D = 1) are shown in Table 1.Table 1
Stability lobe diagrams obtained by the EMHPM, 2nd EFDM, 3rd EFDM, and 4th EFDM under full immersion condition (a/D = 1).
2nd EFDM, 1232 seconds3rd EFDM, 1255 seconds4th EFDM, 1276 secondsEMHPM, 1298 secondsIt is seen from Table1 that the stability lobe diagrams calculated by the EMHPM, 2nd EFDM, 3rd EFDM, and 4th EFDM are consistent with the referenced ones. The results indicate that all these methods are effective in predicting the stability of milling with variable pitch and variable helix cutters. As for the computational time, it is seen from Table 1 that the 4th EFDM takes more time than the 2nd EFDM and 3rd EFDM to obtain stability lobe diagrams. The calculation time depends on the number of matrices. More computing time is needed for calculating more matrices. The higher the order of the method is, the more the matrixes there are. Therefore, 4th EFDM takes more time than the other extended methods.The computational time consumed by extended FDMs is also compared with that consumed by the enhanced multistage homotopy perturbation method (EMHPM). Although more matrixes need to be calculated in the 4th EFDM, one more calculation cycle related to expansion order needs to be calculated in the EMHPM. Generally, the EMHPM takes a little more time to obtain stability lobe diagrams.To further evaluate the different EFDMs, the stability lobe diagrams under low immersion condition (a/D = 0.1) and half immersion (a/D = 0.5) are also obtained, as shown in Table 2.Table 2
Stability lobe diagrams obtained by the 2nd EFDM, 3rd EFDM, and 4th EFDM under half immersion (a/D = 0.5) and low immersion (a/D = 0.1) conditions.
a/D = 0.5a/D = 0.12nd EFDM3rd EFDM4th EFDMFor comparison, the difference between the curves obtained by different EFDMs and the reference curve is also presented, as shown in Figure4.Figure 4
Difference between the curves obtained by different EFDMs and the reference curve: (a)a/D = 0.1; (b) a/D = 0.5; (c) a/D = 1.0.
(a)(b)(c)As shown in Figure4, it can be seen that the difference between the curves obtained by different EFDMs and the reference curve is very small. This result is consistent with the conclusion mentioned above that there is no significant difference between the 2nd EFDM, 3rd EFDM, and 4th EFDM in terms of convergence rate. Therefore, there is no need to extend hypersecond FDMs to analyze the stability in milling with variable pitch and variable helix tools. Besides, higher-order methods require a more computational cost.


## 4.1. Rate of Convergence

The rate of convergence is used to estimate the local error between the absolute value of the maximum-magnitude eigenvalue of the system |μ(k)| and the reference eigenvalue μ0. In this section, the parameter combinations (Ω = 6000 rpm, ap = 0.5 mm) and (Ω = 7000 rpm, ap = 2 mm) are employed to analyze the rates of convergence of the extended FDMs. The radial immersion ratio is set as a/D = 1, the pitch angle combination of the cutter is P = [70°, 110°, 70°, 110°], and the corresponding helix angle combination is β = [30°, 35°, 30°, 35°]. The axial depth of cut is discretized into 20 axial layers, that is, L = 20. To determine the reference eigenvalues for the adopted parameter combinations, the convergences of the eigenvalues obtained by different methods are studied, as shown in Figure 1.Figure 1
The convergences of the eigenvalues obtained by ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM: (a) Ω = 6000 rpm andap = 0.5 mm; (b) Ω = 7000 rpm and ap = 2.0 mm.
(a)(b)It is seen from Figure1 that the ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM all converge to the same reference value when the parameter k reaches 800. In this study, the reference eigenvalue μ0 is determined by using the 2nd EFDM with k = 800. It should be noted that the reference eigenvalue μ0 can also be determined by using the ISDA, 3rd EFDM, and 4th EFDM since all methods converge to the same reference value. After calculation, the reference eigenvalues for the parameter combinations (Ω = 6000 rpm and ap = 0.5 mm) and (Ω = 7000 rpm and ap = 2 mm) are 0.483935 and 0.955073, respectively.In practice, the calculation burden will increase severely when the parameterk is chosen as a very large value. Therefore, the parameter k is usually smaller than that used to determine the value of μ0. With the aim of estimating the rate of convergence of the extended FDMs, the rate of convergence of the improved semidiscretization algorithm (ISDA) [28] is taken as the benchmark. Figure 2 illustrates the rates of convergence of the ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM with different values of parameter k.Figure 2
Rates of convergence of the ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM: (a)ap = 0.5 mm, Ω = 6000 rpm, and |μ0| = 0.483935; (b) ap = 2 mm, Ω = 7000 rpm, and |μ0| = 0.955073.
(a)(b)From Figure2, it is found that the difference between the approximated eigenvalues and the ideals ones calculated by the ISDA has more significant fluctuation than that calculated by the 2nd EFDM, 3rd EFDM, and 4th EFDM. Besides, the eigenvalues obtained by the 2nd EFDM, 3rd EFDM, and 4th EFDM converge to a constant when the parameter k is equal to 300, while the eigenvalues obtained by the ISDA do not converge to a constant under the same condition. It is also seen from Figure 2 that the 2nd EFDM, 3rd EFDM, and 4th EFDM converge faster than the ISDA. The analysis results indicate that the 2nd EFDM, 3rd EFDM, and 4th EFDM are more robust and accurate than the benchmark method.According to Figure2, it is also indicated that there is no significant difference between the 2nd EFDM, 3rd EFDM, and 4th EFDM in terms of the rate of convergence. It is difficult to evaluate the accuracy of different extended FDMs by using the rates of convergence of limited parameter combinations. In the following section, the accuracy of different extended FDMs is further studied by taking the stability lobe diagrams into consideration.


## 4.2. Stability Lobe Diagrams

To verify the effectiveness of the 2nd EFDM, 3rd EFDM, and 4th EFDM, the time-averaged SDM suggested by Sims et al. [20] and the extended SDM suggested by Jin et al. [25] are taken for comparison. The parameters derived from the literature [20] are used for analysis. The stability lobe diagrams obtained by the time-averaged SDM, the extended SDM, the 2nd EFDM, 3rd EFDM, and 4th EFDM are described in Figure 3.Figure 3
Stability lobe diagrams obtained by the time-averaged SDM [20], the extended SDM [25], and the 2nd EFDM, 3rd EFDM, and 4th EFDM.As shown in Figure3, the parameter combinations A (Ω = 8500 rpm and ap = 5 mm) and B (Ω = 8500 rpm and ap = 7 mm) are considered. In the literature [21], the parameter combinations A and B are proved to be a stable case and an unstable case through time-domain simulations, respectively. However, in the stability lobe diagram obtained by the time-averaged SDM, the parameter combination A corresponds to an unstable situation. That is, the time-averaged SDM cannot predict the stability correctly when the spindle speed ranges from 8000 rpm to 9000 rpm. It is seen from Figure 3 that the stabilities of parameter combinations A and B are correctly predicted in the stability lobe diagrams obtained by the 2nd EFDM, 3rd EFDM, and 4th EFDM. Additionally, the stability lobe diagrams obtained by the extended FDMs agree well with that obtained by the extended SDM [25]. Therefore, the 2nd EFDM, 3rd EFDM, and 4th EFDM are proved to be effective methods.The computational time consumed by the 2nd EFDM, 3rd EFDM, and 4th EFDM is compared with that consumed by enhanced multistage homotopy perturbation method (EMHPM) [22] and to illustrate the efficiencies of the extended FDMs. In the calculation process, the pitch angle combination of the cutter is P = [70°, 110°, 70°, 110°], and the corresponding helix angle combination is β = [30°, 35°, 30°, 35°]. The parameter plane, which includes the parameter combinations of spindle speed and axial depth of cut, is divided into a 100 × 100 size grid. The spindle speed ranges from 2 × 103 rpm to 10 × 103 rpm, and the axial depth of cut ranges from 0 mm to 20 mm. The parameter k is chosen as 80, and the axial depth of cut is discretized into 10 axial layers, that is, L = 10. The referenced stability lobe diagrams denoted by the red line are obtained by the 2nd EFDM with k = 300. The stability lobe diagrams obtained by the EMHPM, 2nd EFDM, 3rd EFDM, and 4th EFDM as well as the corresponding computational time under full immersion condition (a/D = 1) are shown in Table 1.Table 1
Stability lobe diagrams obtained by the EMHPM, 2nd EFDM, 3rd EFDM, and 4th EFDM under full immersion condition (a/D = 1).
2nd EFDM, 1232 seconds3rd EFDM, 1255 seconds4th EFDM, 1276 secondsEMHPM, 1298 secondsIt is seen from Table1 that the stability lobe diagrams calculated by the EMHPM, 2nd EFDM, 3rd EFDM, and 4th EFDM are consistent with the referenced ones. The results indicate that all these methods are effective in predicting the stability of milling with variable pitch and variable helix cutters. As for the computational time, it is seen from Table 1 that the 4th EFDM takes more time than the 2nd EFDM and 3rd EFDM to obtain stability lobe diagrams. The calculation time depends on the number of matrices. More computing time is needed for calculating more matrices. The higher the order of the method is, the more the matrixes there are. Therefore, 4th EFDM takes more time than the other extended methods.The computational time consumed by extended FDMs is also compared with that consumed by the enhanced multistage homotopy perturbation method (EMHPM). Although more matrixes need to be calculated in the 4th EFDM, one more calculation cycle related to expansion order needs to be calculated in the EMHPM. Generally, the EMHPM takes a little more time to obtain stability lobe diagrams.To further evaluate the different EFDMs, the stability lobe diagrams under low immersion condition (a/D = 0.1) and half immersion (a/D = 0.5) are also obtained, as shown in Table 2.Table 2
Stability lobe diagrams obtained by the 2nd EFDM, 3rd EFDM, and 4th EFDM under half immersion (a/D = 0.5) and low immersion (a/D = 0.1) conditions.
a/D = 0.5a/D = 0.12nd EFDM3rd EFDM4th EFDMFor comparison, the difference between the curves obtained by different EFDMs and the reference curve is also presented, as shown in Figure4.Figure 4
Difference between the curves obtained by different EFDMs and the reference curve: (a)a/D = 0.1; (b) a/D = 0.5; (c) a/D = 1.0.
(a)(b)(c)As shown in Figure4, it can be seen that the difference between the curves obtained by different EFDMs and the reference curve is very small. This result is consistent with the conclusion mentioned above that there is no significant difference between the 2nd EFDM, 3rd EFDM, and 4th EFDM in terms of convergence rate. Therefore, there is no need to extend hypersecond FDMs to analyze the stability in milling with variable pitch and variable helix tools. Besides, higher-order methods require a more computational cost.


## 5. Conclusions

In this work, the high-order full-discretization methods are extended for stability analysis of milling with variable pitch and variable helix tools. The two-degree-of-freedom milling system is employed to validate the extended methods. The following conclusions can be drawn:(1)
The 2nd EFDM, 3rd EFDM, and 4th EFDM are compared with the benchmark method in terms of the rate of convergence. It is found that the 2nd EFDM, 3rd EFDM, and 4th EFDM converge faster than the benchmark method.(2)
The 4th EFDM takes more time than the 2nd EFDM and 3rd EFDM to obtain stability lobe diagrams. Besides, the EMHPM takes a little more time than the 4th EFDM to obtain stability lobe diagrams.(3)
The difference between the curves obtained by different EFDMs and the reference curve is very small. There is no need to extend hypersecond FDMs to analyze the stability of milling with variable pitch and variable helix tools.



---

*Source: 4517969-2020-07-26.xml*