#### 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 *j*th tooth on the *l*th (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 *j*th tooth, *R* is the radius of the tool shank envelope, is the pitch angle between the *j*th 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 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 *j*th tooth and (*j*−1)th tooth on the *l*th 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 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 aswhere int () is the function that rounds positive numbers towards zero.

Then, equation (7) is integrated on the *i*th 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 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 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 *μ*_{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 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.

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

^{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, 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 *a*D = 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–2*a*D); % exit angle elseif up_or_down = = −1% down-milling fist = acos(2*a*D−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 = *i*dtr; 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) *D*1 = zeros(4l_max + 4,4l_max+4); vlow1 = ones(4l_max,1); *D*1 = D1 + diag(vlow1,−4); *D*1(1 : 4,1 : 4) = M11; *D*1(1 : 4,5 : 8) = M12; *D*1(1 : 4,9 : 12) = M13; *D*1(1 : 4,13 : 16) = M14; *M*2 = zeros(4∗l_max + 4,4 l_max + 4); for *j* = 1 : N *D*2 = 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]; *D*3 = 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,(4*s* + 1):4*s* + 4) = −PiFi_lj; *D*2 = D2 + D3; end *M*2 = M2 + D2; end *D* = *D*1 + M2; FFi = *D*FFi; 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).