Mathematical Problems in Engineering
Volume 2020 (2020), Article ID 4517969, 14 pages
https://doi.org/10.1155/2020/4517969

Research Article

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

Yang Zhang, Kenan Liu, Wuyun Zhao, Wei Zhang, and Fei Dai

Mechanical and Electrical Engineering College, Gansu Agricultural University, Lanzhou 730070, China

Correspondence should be addressed to Kenan Liu; 85263795@qq.com

Received 12 April 2020; Revised 23 June 2020; Accepted 4 July 2020; Published 26 July 2020

Guest Editor: Cuimei Jiang

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.

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.

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 Section 2, 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 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 as

where Ω 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 (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, , can be given by

where 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 by

where 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 two-degree-of-freedom 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 state-space form as follows:

where

To solve equation (7) numerically, the period T is equally divided into k small-time intervals firstly; thus, , where k is an integer and Δt is the length of the small-time interval. The integer related to the time delay can be expressed as

where int () is the function that rounds positive numbers towards zero.

Then, equation (7) is integrated on the ith small-time 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 fourth-order interpolation polynomials are used to approximate the state term , respectively. The periodic coefficient term is approximated by first-order interpolation polynomial, and the time-delay 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 fourth-order and first-order 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 time-delay term is approximated by using the weights and , which can be written as

where the weights and are calculated as follows:

Substituting equations (12)–(16) into equation (11) yields

where

where denotes the identity matrix and the matrixes F₁–F₆ 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 two-degree-of-freedom 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 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  = 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, 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 μ₀. 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.

(a)

(b)

Figure 1: The convergences of the eigenvalues obtained by ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM: (a) Ω = 6000 rpm and ap = 0.5 mm; (b) Ω = 7000 rpm and ap = 2.0 mm.

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 μ₀ is determined by using the 2nd EFDM with k = 800. It should be noted that the reference eigenvalue μ₀ 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 parameter k is chosen as a very large value. Therefore, the parameter k is usually smaller than that used to determine the value of μ₀. 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)

Figure 2: Rates of convergence of the ISDA, 2nd EFDM, 3rd EFDM, and 4th EFDM: (a) ap = 0.5 mm, Ω = 6000 rpm, and |μ₀| = 0.483935; (b) ap = 2 mm, Ω = 7000 rpm, and |μ₀| = 0.955073.

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 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 Figure 3, 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 × 10³ rpm to 10 × 10³ 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 seconds

3rd EFDM, 1255 seconds

4th EFDM, 1276 seconds

EMHPM, 1298 seconds

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.

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.5
a/D = 0.1

2nd EFDM

3rd EFDM

4th EFDM

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)

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.

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, 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.

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:up-milling, −l:down-milling if up_or_down = = 1 % up-milling     fist = 0; % start angle     fiex = acos(1–2aD); % exit angle elseif up_or_down = = −1% down-milling     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 = 0e-3; % starting depth of cut (m) w_fi = 20e-3; % 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)(_fi-w_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 = ti-pitch-(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_fi-o_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.5dtr-dt_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,(4s-3):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. GAU-KYQD-2018-29, 2017RCZX-21, and 2017RCZX-35).

References

Y. Altintas and E. Budak, “Analytical prediction of stability lobes in milling,” CIRP Annals-Manufacturing Technology, vol. 44, no. 1, pp. 357–362, 1995.

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.

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.

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.

T. Insperger and G. Stépán, “Updated semi-discretization method for periodic delay-differential equations with discrete delay,” International Journal for Numerical Methods in Engineering, vol. 61, no. 1, pp. 117–141, 2004.

T. Insperger, G. Stépán, and J. Turi, “On the higher-order semi-discretizations for periodic delayed systems,” Journal of Sound and Vibration, vol. 313, no. 1-2, pp. 334–341, 2008.

Y. Ding, L. Zhu, X. Zhang, and H. Ding, “A full-discretization method for prediction of milling stability,” International Journal of Machine Tools and Manufacture, vol. 50, no. 5, pp. 502–509, 2010.

Y. Ding, L. Zhu, X. Zhang, and H. Ding, “Second-order full-discretization method for milling stability prediction,” International Journal of Machine Tools and Manufacture, vol. 50, no. 10, pp. 926–932, 2010.

Q. Guo, Y. Sun, and Y. Jiang, “On the accurate calculation of milling stability limits using third-order full-discretization method,” International Journal of Machine Tools and Manufacture, vol. 62, pp. 61–66, 2012.

C. G. Ozoegwu, “Least squares approximated stability boundaries of milling process,” International Journal of Machine Tools and Manufacture, vol. 79, pp. 24–30, 2014.

C. G. Ozoegwu, S. N. Omenyi, and S. M. Ofochebe, “Hyper-third order full-discretization methods in milling stability prediction,” International Journal of Machine Tools and Manufacture, vol. 92, pp. 1–9, 2015.

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.

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.

J. Niu, Y. Ding, L. Zhu, and H. Ding, “Runge-Kutta methods for a semi-analytical prediction of milling stability,” Nonlinear Dynamics, vol. 76, no. 1, pp. 289–304, 2014.

H. Li, Y. Dai, and Z. Fan, “Improved precise integration method for chatter stability prediction of two-DOF milling system,” The International Journal of Advanced Manufacturing Technology, vol. 101, no. 5-8, pp. 1235–1246, 2019.

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.

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. 9-12, pp. 2535–2544, 2017.

D. Olvera, G. Urbikain, A. Elías-Zuñiga, and L. López de Lacalle, “Improving stability prediction in peripheral milling of Al7075T6,” Applied Sciences, vol. 8, no. 8, p. 1316, 2018.

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.

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. 3-5, pp. 664–686, 2008.

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.

F. I. Compeán, D. Olvera, F. J. Campa, L. N. López de Lacalle, A. Elías-Zúñ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.

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.

G. Jin, Q. Zhang, H. Qi, and B. Yan, “A frequency-domain 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.

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.

Q. Guo, Y. Jiang, B. Zhao, and P. Ming, “Chatter modeling and stability lobes predicting for non-uniform helix tools,” The International Journal of Advanced Manufacturing Technology, vol. 87, no. 1-4, pp. 251–266, 2016.

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.

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.

J. Niu, Y. Ding, L. Zhu, and H. Ding, “Mechanics and multi-regenerative stability of variable pitch and variable helix milling tools considering runout,” International Journal of Machine Tools and Manufacture, vol. 123, pp. 129–145, 2017.

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.

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.

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.

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. 9-12, pp. 3167–3185, 2017.

A. Otto, S. Rauh, S. Ihlenfeldt, and G. Radons, “Stability of milling with non-uniform pitch and variable helix tools,” The International Journal of Advanced Manufacturing Technology, vol. 89, no. 9-12, pp. 2613–2625, 2017.

Q. Guo, B. Zhao, Y. Jiang, and W. Zhao, “Cutting force modeling for non-uniform helix tools based on compensated chip thickness in five-axis flank milling process,” Precision Engineering, vol. 51, pp. 659–681, 2018.

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 process-machine model,” The International Journal of Advanced Manufacturing Technology, vol. 106, no. 7-8, pp. 2779–2791, 2020.

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.

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.

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.

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.

J. Mei, M. Luo, J. Guo, H. Li, and D. Zhang, “Analytical modeling, design and performance evaluation of chatter-free milling cutter with alternating pitch variations,” IEEE Access, vol. 6, pp. 32367–32375, 2018.

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. 5-8, pp. 2553–2565, 2019.

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.

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.

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.

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.

A. Hosseini, B. Moetakef-Imani, 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.

F. Tehranizadeh and E. Budak, “Design of serrated end mills for improved productivity,” Procedia CIRP, vol. 58, pp. 493–498, 2017.

Y. Guo, B. Lin, and W. Wang, “Modeling of cutting forces with a serrated end mill,” Mathematical Problems in Engineering, vol. 2019, 2019.

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.

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.

G. U. Pelayo and D. O. Trejo, “Model-based phase shift optimization of serrated end mills: minimizing forces and surface location error,” Mechanical Systems and Signal Processing, vol. 144, p. 106860, 2020.