Research Article  Open Access
Yang Zhang, Kenan Liu, Wuyun Zhao, Wei Zhang, Fei Dai, "Stability Analysis for Milling Process with Variable Pitch and Variable Helix Tools by HighOrder FullDiscretization Methods", Mathematical Problems in Engineering, vol. 2020, Article ID 4517969, 14 pages, 2020. https://doi.org/10.1155/2020/4517969
Stability Analysis for Milling Process with Variable Pitch and Variable Helix Tools by HighOrder FullDiscretization Methods
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 highorder fulldiscretization methods (FDMs) to predict the stability of milling with variable pitch and variable helix tools. The timeperiodic delaydifferential 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 statespace equation. Meanwhile, the spindle rotational period is divided into many smalltime intervals, and the state space equation is integrated on the smalltime interval. Then, the highorder interpolation polynomials are used to approximate the state term, and the weights related to the time delay are employed to approximate the timedelay term. The secondorder, thirdorder, and fourthorder 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.
1. Introduction
High precision products machined by using highspeed 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 timeperiodic 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 zerothorder approximation (ZOA) method that determines the stability lobe diagram in the frequency domain. In this method, the timevarying 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 socalled 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 timedomain 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 firstorder semidiscretization method (1st SDM) [6]. The abovementioned timedomain 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 fulldiscretization method (FDM) based on the direct integration scheme. Then, on the framework of FDM, Ding et al. [8] introduced the secondorder FDM (2nd FDM), Guo et al. [9] proposed the thirdorder FDM (3rd FDM), and Ozoegwu reported the least squares approximation methods [10] and hyperthirdorder fulldiscretization methods [11]. Besides, the numerical integration method [12] and differential quadrature method [13] developed by Ding and his coworkers, the Runge–Kuttabased methods proposed by Niu et al. [14], the improved precise integration method proposed by Li et al. [15], the Simpsonbased 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 downmilling operation in a thinwalled 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, timeaveraged 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 frequencydomain 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 variablehelix 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 socalled 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 compensatedchip thicknessbased cutting force model for nonuniform helix tools in the fiveaxis milling process. Cai et al. [36] proposed an integrated processmachine 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 chatterfree material removal rate without using timeconsuming 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 fulldiscretization method peaking at fourth order and then declined at the fifth order. That is, in the hyperfourthorder fulldiscretization methods, the accuracy of the stability results may be affected by the Runge phenomenon. Although it is indicated from the literature [11] that the fourthorder 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 highorder 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 Section 2, the mathematical model of milling dynamics is introduced; in Section 3, the highorder 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 cut is equally discretized into L axial layers with the thickness of ; thus, .
2.1. Angular Position
The angular position of the jth tooth on the lth (1 ≤ l ≤ L) axial layer is represented aswhere Ω is the spindle speed in rpm, is the axial depth of cut, is the helix angle of the jth tooth, R is the radius of the tool shank envelope, is the pitch angle between the jth tooth and (j1)th tooth, and N is the number of cutter tooth.
2.2. The Dynamic Chip Thickness
The dynamic chip thickness caused by regenerative effect, , can be given bywhere and are the displacements of the current tooth in the X and Y directions, respectively; and are the displacements of the previous tooth in the X and Y directions, respectively; and the time delay depends on the separation angle between the jth tooth and (j−1)th tooth on the lth layer, . Thus, the time delay can be given bywhere the period T is equal to the spindle rotational period, namely, . The separation angle can be expressed as follows [24, 29]:
2.3. Milling Dynamics
The governing equation of milling dynamics for a twodegreeoffreedom milling system can be written as follows [5]:where , , and are the damping ratio, angular natural frequency, and modal mass of the milling system in the X direction, respectively; and , , and are the damping ratio, angular natural frequency, and modal mass of the milling system in the Y direction, respectively. The specific cutting force coefficients , , , and are given as follows:where and are the tangential and normal cutting coefficients, respectively. The window function is used to determine whether the tooth is in or out of cut [5].
3. The Proposed Methods
Let , and equation (5) can be represented in statespace form as follows:where
To solve equation (7) numerically, the period T is equally divided into k smalltime intervals firstly; thus, , where k is an integer and Δt is the length of the smalltime interval. The integer related to the time delay can be expressed aswhere int () is the function that rounds positive numbers towards zero.
Then, equation (7) is integrated on the ith smalltime interval resulting in
Equation (10) can be equivalently written as follows [7]:
In equation (11), and denote and , 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 fourthorder interpolation polynomials are used to approximate the state term , respectively. The periodic coefficient term is approximated by firstorder interpolation polynomial, and the timedelay term 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 term and the periodic coefficient term can be, respectively, approximated by the fourthorder and firstorder interpolation polynomials using the Lagrange interpolation scheme. In the calculation process, the interpolation nodes and the corresponding node values , , , , are employed to approximate the state term. Meanwhile, the interpolation nodes , and the node values , are employed to approximate the periodic coefficient term. Then, the state term and the periodic coefficient term can be expressed as follows:
The timedelay term is approximated by using the weights and , which can be written aswhere the weights and are calculated as follows:
Substituting equations (12)–(16) into equation (11) yieldswherewhere denotes the identity matrix and the matrixes F_{1}–F_{6} can be expressed as follows:
According to equation (17), the following local discrete map can be obtained:where denotes the maximum value of , and matrix can be expressed as follows:where
The positions of the matrices and in equation (21) depend on the value of related to the time delay . As for the twodegreeoffreedom milling system, the matrix is located in columns (4–3) to (4) of the matrix , and the matrix is located in columns (4 + 1) to (4 + 4) of the matrix . Additionally, the dimension of the matrix depends on the value of .
Then, the Floquet state transition matrix ψ is obtained as
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:
4. Stability Analysis
The twodegreeoffreedom 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 = 1.4986 kg and = 1.199 kg, the damping ratios are = 0.055801 and = 0.025004, the angular natural frequencies are = 3540.89 rad/s and = 3243.44 rad/s, and the tangential and the normal cutting force coefficients are = 697 MPa and = 256 MPa, downmilling.
4.1. Rate of Convergence
The rate of convergence is used to estimate the local error between the absolute value of the maximummagnitude eigenvalue of the system μ(k) and the reference eigenvalue μ_{0}. In this section, the parameter combinations (Ω = 6000 rpm, a_{p} = 0.5 mm) and (Ω = 7000 rpm, a_{p} = 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.
(a)
(b)
It is seen from Figure 1 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 a_{p} = 0.5 mm) and (Ω = 7000 rpm and a_{p} = 2 mm) are 0.483935 and 0.955073, respectively.
In practice, the calculation burden will increase severely when the parameter k 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.
(a)
(b)
From Figure 2, 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 Figure 2, 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 timeaveraged 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 timeaveraged SDM, the extended SDM, the 2nd EFDM, 3rd EFDM, and 4th EFDM are described in Figure 3.
As shown in Figure 3, the parameter combinations A (Ω = 8500 rpm and a_{p} = 5 mm) and B (Ω = 8500 rpm and a_{p} = 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 timedomain simulations, respectively. However, in the stability lobe diagram obtained by the timeaveraged SDM, the parameter combination A corresponds to an unstable situation. That is, the timeaveraged 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 × 10^{3} rpm to 10 × 10^{3} 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.

It is seen from Table 1 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.

For comparison, the difference between the curves obtained by different EFDMs and the reference curve is also presented, as shown in Figure 4.
(a)
(b)
(c)
As shown in Figure 4, 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, higherorder methods require a more computational cost.
5. Conclusions
In this work, the highorder fulldiscretization methods are extended for stability analysis of milling with variable pitch and variable helix tools. The twodegreeoffreedom 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.
Appendix
Matlab code for the stability analysis by using the 4th EFDM. close all; clear all; clc N = 4; % number of tool teeth R = 0.009525; % radius (m) P = [70pi/180,110pi/180,70pi/180,110pi/180]; % pitch angle combination beta = [30pi/180,35pi/180,30pi/180,35pi/180]; % helix angle combination aD = 1.0; % radial immersion ratio a/D Kt = 6.97e8; % tangential cutting force coefficient (N/m^2) Kn = 2.56e8; % normal cutting force coefficient (N/m^2) wnx = 563.552pi; % angular natural frequency (rad/s) wny = 516.212pi; % angular natural frequency (rad/s) m_tx = 1.4986; % mass in x direction (kg) m_ty = 1.199; % mass in y direction (kg) zeta_x = 0.055801; % relative damping in x direction zeta_y = 0.025004; % relative damping in y direction up_or_down = −1; % l:upmilling, −l:downmilling if up_or_down = = 1 % upmilling fist = 0; % start angle fiex = acos(1–2aD); % exit angle elseif up_or_down = = −1% downmilling fist = acos(2aD−1); % start angle fiex = pi; % exit angle end stx = 100; % steps of spindle speed sty = 100; % steps of depth of cut w_st = 0e3; % starting depth of cut (m) w_fi = 20e3; % final depth of cut (m) o_st = 2e3; % starting spindle speed (m) o_fi = 10e3; % final spindle speed (rpm) k = 100; % discrete number for spindle speed L = 10; % discrete number for axial depth of cut dtr = 2pi/k; intk = 10; for y = 1:sty + 1 w(y) = _st + (y−1)(_fiw_st)/sty; for i = 1:k + 1 h2xx(y, i) = 0; h2xy(y, i) = 0; h2yx(y,i) = 0; h2yy(y,i) = 0; ti = idtr; pitch = 0; for j = 1 : N h1xx(y,i,j) = 0; h1xy(y,i,j) = 0; h1yx(y,i,j) = 0; h1yy(y,i,j) = 0; if j = = 1 pitch = 0; else pitch = pitch + P(j); end for n = 1 : L hxx(y,i,j,n) = 0; hxy(y,i,j,n) = 0; hyx(y,i,j,n) = 0; hyy(y,i,j,n) = 0; if j = = 1 dt_helix = (n−1/2)((y)/L)(tan(beta(j))−tan(beta(N)))/R; else dt_helix=(n−1/2)((y)/L)(tan(beta(j))−tan(beta(j−1)))/R; end dt_phi(y,j,n) = P(j) + dt_helix; l_max = fix((max(max(max(dt_phi)))+dtr/2)/dtr); for hh = 1:intk fi = tipitch(n−1)((y)/L)tan(beta(j))/R + hhdtr/intk; fi = mod(fi,2pi); if (fi ≥ fist)(fi ≤ fiex) = 1; else = 0; end hxx(y,i,j,n) = hxx(y,i,j,n) + (Ktcos(fi) + Knsin(fi))sin(fi)((y)/L)/intk; hxy(y,i,j,n) = hxy(y,i,j,n) + (Ktcos(fi) + Knsin(fi))cos(fi)((y)/L)/intk; hyx(y,i,j,n) = hyx(y,i,j,n) + (−Ktsin(fi) + Kncos(fi))sin(fi)((y)/L)/intk; hyy(y,i,j,n) = hyy(y,i,j,n) + (−Ktsin(fi) + Kncos(fi))cos(fi)((y)/L)/intk; end h1xx(y,i,j) = h1xx(y,i,j) + hxx(y,i,j,n); h1xy(y,i,j) = h1xy(y,i,j) + hxy(y,i,j,n); h1yx(y,i,j) = h1yx(y,i,j) + hyx(y,i,j,n); h1yy(y,i,j) = h1yy(y,i,j) + hyy(y,i,j,n); end h2xx(y,i) = h2xx(y,i) + h1xx(y,i,j); h2xy(y,i) = h2xy(y,i) + h1xy(y,i,j); h2yx(y,i) = h2yx(y,i) + h1yx(y,i,j); h2yy(y,i) = h2yy(y,i) + h1yy(y,i,j); end end end A = [0 0 1 0; 0 0 0 1; −wnx^2 0–2zeta_xwnx 0; 0 −wny^2 0–2zeta_ywny]; invA = inv(A); I = eye(4,4); for x = 1:stx + 1 o = o_st + (x−1)(o_fio_st)/stx; T = 60/o; % spindle rotational period dt = T/k; % time step F0 = expm(Adt); F1 = (F0–I)invA; F2 = (F1−dtI)invA; F3 = (2F2−dt^2I)invA; F4 = (3F3−dt^3I)invA; F5 = (4F4−dt^4I)invA; F6 = (5F5−dt^5I)invA; G11 = (−1/24)/(dt^5)F6−1/(24dt^4)F5 + 1/(8dt^3)F4 + 1/(24dt^2)F3−1/(12dt)F2; G12 = (1/24)/(dt^5)F6 + 1/(12dt^4)F5−1/(24dt^3)F4−1/(12dt^2)F3; G13 = (1/6)/(dt^5)F6 + 1/(3dt^4)F5−2/(3dt^3)F4−1/(3dt^2)F3 + 1/(2dt)F2; G14 = (−1/6)/(dt^5)F6−1/(2dt^4)F5 + 1/(6dt^3)F4 + 1/(2dt^2)F3; G15 = (−1/4)/(dt^5)F6−3/(4dt^4)F5 + 3/(4dt^3)F4 + 7/(4dt^2)F3−3/(2dt)F2; G16 = (1/4)/(dt^5)F6 + 1/(dt^4)F5 + 1/(4dt^3)∗F4−3/(2dt^2)F3; G17 = (1/6)/(dt^5)F6 + 2/(3dt^4)F5−5/(3dt^2)F3−1/(6dt)F2 + F1; G18 = (−1/6)/(dt^5)F6−5/(6dt^4)F5−5/(6dt^3)F4 + 5/(6dt^2)F3 + 1/(dt)F2; G19 = (−1/24)/(dt^5)F6−5/(24dt^4)F5−5/(24dt^3)F4 + 5/(24dt^2)F3 + 1/(4dt)F2; G110 = (1/24)/(dt^5)F6 + 1/(4dt^4)F5 + 11/(24dt^3)F4 + 1/(4dt^2)F3; for y = 1 : sty + 1 % sweeping depth of cuts FFi = eye(4l_max + 4,4l_max + 4); % construct transition matrix FFi for i = 1 : k Ai = [0,0,0,0; 0,0,0,0; h2xx(y,i)/m_tx,h2xy(y,i)/m_tx,0,0; h2yx(y,i)/m_ty,h2yy(y,i)/m_ty,0,0]; Aii = [0,0,0,0; 0,0,0,0; h2xx(y,i + 1)/m_tx,h2xy(y,i + 1)/m_tx,0,0; h2yx(y,i + 1)/m_ty, h2yy(y,i + 1)/m_ty,0,0]; Pi = inv(I−G110Aii−G19Ai); %x(i + 1) M11 = Pi(F0 + G18Aii + G17Ai); % X(i) M12 = Pi(G16Aii + G15Ai); % X(i−1) M13 = Pi(G14Aii + G13Ai); % X(i−2) M14 = Pi(G12Aii + G11Ai); % X(i−3) D1 = zeros(4l_max + 4,4l_max+4); vlow1 = ones(4l_max,1); D1 = D1 + diag(vlow1,−4); D1(1 : 4,1 : 4) = M11; D1(1 : 4,5 : 8) = M12; D1(1 : 4,9 : 12) = M13; D1(1 : 4,13 : 16) = M14; M2 = zeros(4∗l_max + 4,4 l_max + 4); for j = 1 : N D2 = zeros(4l_max + 4,4l_max + 4); for n = 1 : L Bi = [0,0,0,0; 0,0,0,0; hxx(y,i,j,n)/m_tx,hxy(y,i,j,n)/m_tx,0,0; hyx(y,i,j,n)/m_ty,hyy(y,i,j,n)/m_ty,0,0]; Bii = [0,0,0,0; 0,0,0,0; hxx(y,i + 1,j,n)/m_tx,hxy(y,i + 1,j,n)/m_tx,0,0; hyx(y,i + 1,j,n)/m_ty, hyy(y,i + 1,j,n)/m_ty,0,0]; D3 = zeros(4l_max + 4,4l_max + 4); s = fix((dt_phi(y,j,n) + 0.5dtr)/dtr); w_a = (sdtr + 0.5dtrdt_phi(y,j,n))/dtr; G21 = F3/(dt^2) + (−2/dt + _a/dt)F2 + (−w_a+1)F1; %x(i−k)B(i) G22 = −F3/(dt^2) + (1/dt−w_a/dt)F2; %x(i−k)B(i + 1) G23 = −F3/(dt^2) + (1/dt−w_a/dt)F2 + _aF1; %x(i−k + 1)B(i) G24 = F3/(dt^2) + (_a/dt)F2; %x(i−k + 1) (i + 1) Fi_lj = G21Bi + G22Bii; Pi_lj = G23Bi + G24Bii; D3(1 : 4,(4s3):4s) = −PiPi_lj; D3(1 : 4,(4s + 1):4s + 4) = −PiFi_lj; D2 = D2 + D3; end M2 = M2 + D2; end D = D1 + M2; FFi = DFFi; end ss(x,y) = o; % matrix of spindle speeds dc(x,y) = (y); % matrix of depth of cuts ei(x,y) = max(abs(eig(FFi))); % matrix of eigenvalues end % End of sweeping depth of cuts stx + 1−x end figure; contour (ss,dc,ei, [1,1], “k”), xlabel(“(rmp)”), ylabel(“w(m)”)
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by the Special Fund for Talents of Gansu Agricultural University (Grant nos. GAUKYQD201829, 2017RCZX21, and 2017RCZX35).
References
 Y. Altintas and E. Budak, “Analytical prediction of stability lobes in milling,” CIRP AnnalsManufacturing Technology, vol. 44, no. 1, pp. 357–362, 1995. View at: Publisher Site  Google Scholar
 S. D. Merdol and Y. Altintas, “Multi frequency solution of chatter stability for low immersion milling,” Journal of Manufacturing Science and Engineering, vol. 126, no. 3, pp. 459–466, 2004. View at: Publisher Site  Google Scholar
 P. V. Bayly, J. E. Halley, B. P. Mann, and M. A. Davies, “Stability of interrupted cutting by temporal finite element analysis,” Journal of Manufacturing Science and Engineering, vol. 125, no. 2, pp. 220–225, 2003. View at: Publisher Site  Google Scholar
 E. A. Butcher, O. A. Bobrenkov, E. Bueler, and P. Nindujarla, “Analysis of milling stability by the Chebyshev collocation method: algorithm and optimal stable immersion levels,” Journal of Computational and Nonlinear Dynamics, vol. 4, no. 3, Article ID 031003, 2009. View at: Publisher Site  Google Scholar
 T. Insperger and G. Stépán, “Updated semidiscretization method for periodic delaydifferential equations with discrete delay,” International Journal for Numerical Methods in Engineering, vol. 61, no. 1, pp. 117–141, 2004. View at: Publisher Site  Google Scholar
 T. Insperger, G. Stépán, and J. Turi, “On the higherorder semidiscretizations for periodic delayed systems,” Journal of Sound and Vibration, vol. 313, no. 12, pp. 334–341, 2008. View at: Publisher Site  Google Scholar
 Y. Ding, L. Zhu, X. Zhang, and H. Ding, “A fulldiscretization method for prediction of milling stability,” International Journal of Machine Tools and Manufacture, vol. 50, no. 5, pp. 502–509, 2010. View at: Publisher Site  Google Scholar
 Y. Ding, L. Zhu, X. Zhang, and H. Ding, “Secondorder fulldiscretization method for milling stability prediction,” International Journal of Machine Tools and Manufacture, vol. 50, no. 10, pp. 926–932, 2010. View at: Publisher Site  Google Scholar
 Q. Guo, Y. Sun, and Y. Jiang, “On the accurate calculation of milling stability limits using thirdorder fulldiscretization method,” International Journal of Machine Tools and Manufacture, vol. 62, pp. 61–66, 2012. View at: Publisher Site  Google Scholar
 C. G. Ozoegwu, “Least squares approximated stability boundaries of milling process,” International Journal of Machine Tools and Manufacture, vol. 79, pp. 24–30, 2014. View at: Publisher Site  Google Scholar
 C. G. Ozoegwu, S. N. Omenyi, and S. M. Ofochebe, “Hyperthird order fulldiscretization methods in milling stability prediction,” International Journal of Machine Tools and Manufacture, vol. 92, pp. 1–9, 2015. View at: Publisher Site  Google Scholar
 Y. Ding, L. M. Zhu, X. J. Zhang, and H. Ding, “Numerical integration method for prediction of milling stability,” Journal of Manufacturing Science and Engineering, vol. 133, no. 3, Article ID 031005, 2011. View at: Publisher Site  Google Scholar
 Y. Ding, L. M. Zhu, X. J. Zhang, and H. Ding, “Stability analysis of milling via the differential quadrature method,” Journal of Manufacturing Science and Engineering, vol. 135, no. 4, Article ID 044502, 2013. View at: Publisher Site  Google Scholar
 J. Niu, Y. Ding, L. Zhu, and H. Ding, “RungeKutta methods for a semianalytical prediction of milling stability,” Nonlinear Dynamics, vol. 76, no. 1, pp. 289–304, 2014. View at: Publisher Site  Google Scholar
 H. Li, Y. Dai, and Z. Fan, “Improved precise integration method for chatter stability prediction of twoDOF milling system,” The International Journal of Advanced Manufacturing Technology, vol. 101, no. 58, pp. 1235–1246, 2019. View at: Publisher Site  Google Scholar
 Z. Zhang, H. Li, G. Meng, and C. Liu, “A novel approach for the prediction of the milling stability based on the Simpson method,” International Journal of Machine Tools and Manufacture, vol. 99, pp. 43–47, 2015. View at: Publisher Site  Google Scholar
 X. Zhang, C. Xiong, Y. Ding, and H. Ding, “Prediction of chatter stability in high speed milling using the numerical differentiation method,” The International Journal of Advanced Manufacturing Technology, vol. 89, no. 912, pp. 2535–2544, 2017. View at: Publisher Site  Google Scholar
 D. Olvera, G. Urbikain, A. ElíasZuñiga, and L. López de Lacalle, “Improving stability prediction in peripheral milling of Al7075T6,” Applied Sciences, vol. 8, no. 8, p. 1316, 2018. View at: Publisher Site  Google Scholar
 Y. Altıntas, S. Engin, and E. Budak, “Analytical stability prediction and design of variable pitch cutters,” Journal of Manufacturing Science and Engineering, vol. 121, no. 2, pp. 173–178, 1999. View at: Publisher Site  Google Scholar
 N. D. Sims, B. Mann, and S. Huyanan, “Analytical prediction of chatter stability for variable pitch and variable helix milling tools,” Journal of Sound and Vibration, vol. 317, no. 35, pp. 664–686, 2008. View at: Publisher Site  Google Scholar
 V. Sellmeier and B. Denkena, “Stable islands in the stability chart of milling processes due to unequal tooth pitch,” International Journal of Machine Tools and Manufacture, vol. 51, no. 2, pp. 152–164, 2011. View at: Publisher Site  Google Scholar
 F. I. Compeán, D. Olvera, F. J. Campa, L. N. López de Lacalle, A. ElíasZúñiga, and 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, vol. 57, pp. 27–33, 2012. View at: Publisher Site  Google Scholar
 M. Wan, W.H. Zhang, J.W. Dang, and Y. Yang, “A unified stability prediction method for milling process with multiple delays,” International Journal of Machine Tools and Manufacture, vol. 50, no. 1, pp. 29–41, 2010. View at: Publisher Site  Google Scholar
 G. Jin, Q. Zhang, H. Qi, and B. Yan, “A frequencydomain solution for efficient stability prediction of variable helix cutters milling,” Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 228, no. 15, pp. 2702–2710, 2014. View at: Publisher Site  Google Scholar
 G. Jin, Q. Zhang, S. Hao, and Q. Xie, “Stability prediction of milling process with variable pitch and variable helix cutters,” Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 228, no. 2, pp. 281–293, 2014. View at: Publisher Site  Google Scholar
 Q. Guo, Y. Jiang, B. Zhao, and P. Ming, “Chatter modeling and stability lobes predicting for nonuniform helix tools,” The International Journal of Advanced Manufacturing Technology, vol. 87, no. 14, pp. 251–266, 2016. View at: Publisher Site  Google Scholar
 N. D. Sims, “Fast chatter stability prediction for variable helix milling tools,” Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 230, no. 1, pp. 133–144, 2016. View at: Publisher Site  Google Scholar
 Y. Wang, T. Wang, Z. Yu, Y. Zhang, Y. Wang, and H. Liu, “Chatter prediction for variable pitch and variable helix milling,” Shock and Vibration, vol. 2015, Article ID 419172, 9 pages, 2015. View at: Publisher Site  Google Scholar
 J. Niu, Y. Ding, L. Zhu, and H. Ding, “Mechanics and multiregenerative stability of variable pitch and variable helix milling tools considering runout,” International Journal of Machine Tools and Manufacture, vol. 123, pp. 129–145, 2017. View at: Publisher Site  Google Scholar
 T. Huang, L. Zhu, S. Du, Z. Chen, and H. Ding, “Robust active chatter control in milling processes with variable pitch cutters,” Journal of Manufacturing Science and Engineering, vol. 140, Article ID 101005, 2018. View at: Publisher Site  Google Scholar
 G. Stepan, D. Hajdu, A. Iglesias, D. Takacs, and Z. Dombovari, “Ultimate capability of variable pitch milling cutters,” CIRP Annals, vol. 67, no. 1, pp. 373–376, 2018. View at: Publisher Site  Google Scholar
 J. Huang, P. Deng, H. Li, and B. Wen, “Stability analysis for milling system with variable pitch cutters under variable speed,” Journal of Vibroengineering, vol. 21, no. 2, pp. 331–347, 2019. View at: Publisher Site  Google Scholar
 D. Chen, X. Zhang, Y. Xie, X. Zhang, and H. Ding, “A unified analytical cutting force model for variable helix end mills,” The International Journal of Advanced Manufacturing Technology, vol. 92, no. 912, pp. 3167–3185, 2017. View at: Publisher Site  Google Scholar
 A. Otto, S. Rauh, S. Ihlenfeldt, and G. Radons, “Stability of milling with nonuniform pitch and variable helix tools,” The International Journal of Advanced Manufacturing Technology, vol. 89, no. 912, pp. 2613–2625, 2017. View at: Publisher Site  Google Scholar
 Q. Guo, B. Zhao, Y. Jiang, and W. Zhao, “Cutting force modeling for nonuniform helix tools based on compensated chip thickness in fiveaxis flank milling process,” Precision Engineering, vol. 51, pp. 659–681, 2018. View at: Publisher Site  Google Scholar
 S. Cai, B. Yao, W. Feng, Z. Cai, B. Chen, and Z. He, “Milling process simulation for the variable pitch cutter based on an integrated processmachine model,” The International Journal of Advanced Manufacturing Technology, vol. 106, no. 78, pp. 2779–2791, 2020. View at: Publisher Site  Google Scholar
 E. Budak, “An analytical design method for milling cutters with nonconstant pitch to increase stability, part I: theory,” Journal of Manufacturing Science and Engineering, vol. 125, no. 1, pp. 29–34, 2003. View at: Publisher Site  Google Scholar
 E. Budak, “An analytical design method for milling cutters with nonconstant pitch to increase stability, part 2: application,” Journal of Manufacturing Science and Engineering, vol. 125, no. 1, pp. 35–38, 2003. View at: Publisher Site  Google Scholar
 A. R. Yusoff and N. D. Sims, “Optimisation of variable helix tool geometry for regenerative chatter mitigation,” International Journal of Machine Tools and Manufacture, vol. 51, no. 2, pp. 133–141, 2011. View at: Publisher Site  Google Scholar
 A. Comak and E. Budak, “Modeling dynamics and stability of variable pitch and helix milling tools for development of a design method to maximize chatter stability,” Precision Engineering, vol. 47, pp. 459–468, 2017. View at: Publisher Site  Google Scholar
 J. Mei, M. Luo, J. Guo, H. Li, and D. Zhang, “Analytical modeling, design and performance evaluation of chatterfree milling cutter with alternating pitch variations,” IEEE Access, vol. 6, pp. 32367–32375, 2018. View at: Publisher Site  Google Scholar
 Y. Guo, B. Lin, and W. Wang, “Optimization of variable helix cutter for improving chatter stability,” The International Journal of Advanced Manufacturing Technology, vol. 104, no. 58, pp. 2553–2565, 2019. View at: Publisher Site  Google Scholar
 S. D. Merdol and Y. Altintas, “Mechanics and dynamics of serrated end mills,” in Proceedings of the ASME 2002 International Mechanical Engineering Congress and Exposition, pp. 337–342, New Orleans, LA, USA, November 2002. View at: Publisher Site  Google Scholar
 S. D. Merdol and Y. Altintas, “Mechanics and dynamics of serrated cylindrical and tapered end mills,” Journal of Manufacturing Science and Engineering, vol. 126, no. 2, pp. 317–326, 2004. View at: Publisher Site  Google Scholar
 Z. Dombovari, Y. Altintas, and G. Stepan, “Stability of serrated milling cutters,” in Proceedings of the 12th CIRP Conference on Modelling of Machining Operations, pp. 873–878, Mondragon, Spain, 2009. View at: Google Scholar
 Z. Dombovari, Y. Altintas, and G. Stepan, “The effect of serration on mechanics and stability of milling cutters,” International Journal of Machine Tools and Manufacture, vol. 50, no. 6, pp. 511–520, 2010. View at: Publisher Site  Google Scholar
 A. Hosseini, B. MoetakefImani, and H. A. Kishawy, “Mechanistic modelling for cutting with serrated end mills—a parametric representation approach,” Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, vol. 225, no. 7, pp. 1019–1032, 2011. View at: Publisher Site  Google Scholar
 F. Tehranizadeh and E. Budak, “Design of serrated end mills for improved productivity,” Procedia CIRP, vol. 58, pp. 493–498, 2017. View at: Publisher Site  Google Scholar
 Y. Guo, B. Lin, and W. Wang, “Modeling of cutting forces with a serrated end mill,” Mathematical Problems in Engineering, vol. 2019, 2019. View at: Publisher Site  Google Scholar
 P. Bari, M. Law, and P. Wahi, “Improved chip thickness model for serrated end milling,” CIRP Journal of Manufacturing Science and Technology, vol. 25, pp. 36–49, 2019. View at: Publisher Site  Google Scholar
 F. Tehranizadeh, R. Koca, and E. Budak, “Investigating effects of serration geometry on milling forces and chatter stability for their optimal selection,” International Journal of Machine Tools and Manufacture, vol. 144, p. 103425, 2019. View at: Publisher Site  Google Scholar
 G. U. Pelayo and D. O. Trejo, “Modelbased phase shift optimization of serrated end mills: minimizing forces and surface location error,” Mechanical Systems and Signal Processing, vol. 144, p. 106860, 2020. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2020 Yang Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.