Stability Analysis of Milling Process with Multiple Delays

: Cutting chatter is extremely harmful to the machining process, and it is of great signiﬁcance to eliminate chatter through analyzing the stability of the machining process. In this work, the stability of the milling process with multiple delays is investigated. Considering the regeneration e ﬀ ect, the dynamics of the milling process with variable pitch cutter is modeled as periodic coe ﬃ cients delayed di ﬀ erential equations (DDEs) with multiple delays. An adaptive variable-step numerical integration method (AVSNIM) considering the e ﬀ ect of the helix angle is developed ﬁrstly, which can discretize the cutting period accurately, thereby improving the calculation accuracy of the stability limit of the milling process. The accuracy and e ﬃ ciency of the AVSNIM are veriﬁed through a benchmark milling model. Subsequently, a novel spindle speed-dependent discretization algorithm is proposed, which is combined with the AVSNIM to further reduce the calculation time of the stability lobes diagram (SLD). The simulation experiment results demonstrate that the proposed algorithm can e ﬀ ectively reduce the calculation time.


Introduction
Chatter is a self-excited vibration induced in the machining process which reduces the machining efficiency and surface quality, accelerates the tool wear and shortens the tool durability. There have been a lot of research on the mechanism of chatter [1][2][3][4][5][6], among which regenerative effect was usually considered as a main inducement. Since the regenerative chatter is caused by the phase differences between wavy surfaces left by the adjacent teeth, it can be suppressed by adjusting the pitch angle between the adjacent teeth. The use of variable pitch cutters to improve the stability of the milling process is built on this idea. Different from the uniform pitch cutter, when the variable pitch cutter is used, the dynamics model of cutting vibration changes from DDEs with single delay to DDEs with multiple delays.
The stability analysis of the milling process with multiple delays has been studied through different methods. Altintas et al. [7] adopted the frequency domain method to analyze the milling stability of the variable pitch cutter and presented a method to select the optimal pitch angles. Budak [8,9] proposed an analytical design method for nonconstant pitch milling cutters to improve the stability of the milling process and verified it through cutting experiments. The results showed that chatter stability can be improved significantly even at slow cutting speeds by properly designing the pitch angles. The modified semi-discretization method (SDM) was used by Sims et al. [10] to investigate the stability of variable pitch and variable helix milling tools. Considering the tool runout and the

Milling Model with Multiple Delays
Milling vibration considering the regeneration effect is usually modeled by DDEs which can be written as M ..

y(t)+C
. y(t)+Ky(t) = f(t) (1) where M, C and K are the mass, damping and stiffness matrices respectively, y(t) is the vibration displacements vector and f(t) is the time-varying cutting force vector. When using the variable pitch cutter (as described in Figure 1), the dynamics model of the milling process becomes DDEs with multiple delays. The cutting force f(t) in Equation (1) can therefore be expressed as H j (t) · y(t) − y t − T j (2) where N is the number of teeth, T j is the time delay of the j-th tooth, y(t) and y t − T j are the vibration displacements at the current moment and the previous tooth-passing period of the j-th tooth separately. H j (t) is the time-varying cutting force coefficient matrix of the j-th tooth, which satisfies H j (t) = H j (t + T), where T is the spindle period. H j (t) can be written as where where K t and K n are the tangential and radial cutting force coefficients, respectively. z u, j (t) and z l,j (t) are the upper and lower bounds of the cutting edge participating in cutting on the j-th tooth as shown in Figure 1c, respectively. ϕ j (t, z) is the angular position of the differential element of cutting edge with a height z on the j-th tooth, which can be expressed as ϕ j (t, z) = ϕ j (0, 0) + 2πΩ 60 · t − φ lag (z) (5) where ϕ j (0, 0) is the angular position of the lowest differential element on the j-th tooth at the initial time, Ω is the spindle speed, φ lag (z) is the lag angle of the differential element at height z as shown in Figure 1d, and φ lag (z) = (2 tan β/D) · z, where D is the diameter of the cutter, β is the helix angle. Equation (1) can be converted to the state-space form as .
where A is a constant matrix, representing the time-invariant characteristics of the milling system, while B j is a periodic matrix satisfying B j (t) = B j (t + T), which describes the periodic characteristics of the system. x is known as the state vector

y(t)
where I is an identity matrix with the same size as M.
Obviously, Equation (6) describes a linear periodic system. According to the Floquet theory, the stability of a linear periodic system depends on the spectral radius of the Floquet transition matrix: if the spectral radius is less than one, then the system is asymptotically stable; if it is larger than one, the system is unstable; if the spectral radius equals to one, then the system is in a critical stable state.  Equation (1) can be converted to the state-space form as where A is a constant matrix, representing the time-invariant characteristics of the milling system, while is a periodic matrix satisfying , which describes the periodic characteristics of the system. x is known as the state vector where I is an identity matrix with the same size as M .
Obviously, Equation (6) describes a linear periodic system. According to the Floquet theory, the stability of a linear periodic system depends on the spectral radius of the Floquet transition matrix: if the spectral radius is less than one, then the system is asymptotically stable; if it is larger than one, the system is unstable; if the spectral radius equals to one, then the system is in a critical stable state.

AVSNIM Considering the Helix Angle
The variable-step technique can discrete the milling period accurately. Considering the effect of helix angle, an AVSNIM for the stability analysis of milling process with multiple delays is proposed in this section.

Algorithm Description
For the milling process with multiple delays, the cutters may experience several free and forced vibrations in one spindle period. When the hysteresis effect caused by the helix angle is considered, the free vibration and forced vibration intervals will vary with the axial depth of the cut.
In order to achieve accurate discretization, the first thing is to determine the free vibration and forced vibration interval of the cutter in one spindle period T. Let φ j be the pitch angle between the (j − 1)-th and the j-th tooth, w be the axial depth of cut, ϕ st and ϕ ex be the tool start angle and exit angle, respectively. The angle swept by each tooth from cutting into the workpiece to cutting out of the workpiece can be expressed as Then, the free vibration angle interval before the j-th tooth cutting into the workpiece φ f r,j can be expressed as φ f r,j (w) = φ j − φ s (w) (8) In general, the distribution of free vibration and forced vibration in one spindle period can be divided into two categories: one where, at most, one tooth is in cutting during one spindle period; the other where more than one tooth is in cutting simultaneously at a certain time in one spindle period. If φ f r, j (w) > 0 holds for all j ( j = 1, 2, · · · , N), it corresponds to the first case. There are N free vibration intervals and N forced vibration intervals distributed alternately in one spindle period and each forced vibration angle interval φ f o,j (w) equals to φ s (w), as shown in Figure 2a. Else if φ f r, j (w) ≤ 0, it means that the j-th tooth has already cut into the workpiece before the (j − 1)-th tooth cuts out of the workpiece, which belongs to the second case. Under these circumstances, the continuous forced vibration angle intervals are combined as one forced vibration angle interval, as shown in Figure 2b,c. The length of this continuous forced vibration angle interval φ f oc (w) can be calculated as According to the above analysis, one spindle period is finally divided into N d (N d ≤ N) free vibration and N d forced vibration intervals that are alternately distributed. For the convenience of description, we mark one pair of adjacent free vibration and forced vibration intervals as one segment, hence one whole spindle period is divided into N d segments, as described in Figure 2c. More than one tooth is in cutting simultaneously Figure 2. The distribution of free vibration and forced vibration angle interval in one spindle period, (a) at most one tooth is in cutting; (b) more than one tooth is in cutting simultaneously; (c) combine the continuous forced vibration angle interval.
Once the distribution of free vibration and forced vibration intervals of the cutter is obtained, the spindle period can be discretized according to the following rules. In the free vibration interval, the vibration displacement of the cutter can be solved analytically, so there is no need to discretize. For the forced vibration interval, the following strategies are adopted for discretization.
Firstly, the sum of all forced vibration angle intervals of the cutter is calculated as  Once the distribution of free vibration and forced vibration intervals of the cutter is obtained, the spindle period can be discretized according to the following rules. In the free vibration interval, the vibration displacement of the cutter can be solved analytically, so there is no need to discretize. For the forced vibration interval, the following strategies are adopted for discretization.
Firstly, the sum of all forced vibration angle intervals of the cutter is calculated as Then, φ f o (w) is discretized into m p intervals; each interval has a size of ∆φ p (w), and where m p is the pre-set discretization parameter. Next, the discrete number m i (w) and angle step ∆φ p,i (w) of each forced vibration angle interval φ f o,i (w) can be obtained by the following formula where "round (*)" means rounds "*" towards the nearest integer, i = 1, 2, · · · , N d .
The discrete time step can be easily converted from the discrete angle step through the following formula To date, we have got all the discretization information in one spindle period, including the number of discrete intervals and the length of each interval. The number of discrete intervals is and therefore the index of the discrete time nodes in one spindle period are 1, 2, · · · , m(w), m(w) + 1 . Let n(i, k) be the index of the k-th node of the i-th segment in the discrete time nodes of one spindle period, where k = 1, 2, · · · , 1 + m i (w) + 1, n(i, k) ∈ 1, 2, · · · , m(w), m(w) + 1 and can be expressed as According to the theory of differential equation, the state vector at the k-th time node of the i-th segment x t n(i,k) can be deduced as Exchanging the order of integration and summation in Equation (15), and using the trapezoidal rule to approximate the integral term [26], the state vector x t n(i,k) at forced vibration time nodes can be expressed as To get the Floquet transition matrix, the delay terms x t n(i,k−1) − T j and x t n(i,k) − T j should be interpolated on two adjacent spindle periods, as shown in Figure 3.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 23 can be expressed as the following unified form  are in the free vibration interval, the interpolation coefficients in Equation (18) can be expressed as are in the forced vibration intervals, the interpolation coefficients can be obtained through Lagrange interpolation and expressed as Substituting Equation (18) into Equation (16), and rewriting it as follows Let n i, k, T j be the index of the time node on two adjacent spindle periods and satisfy the following relationship x t n(i,k) − T j can be expressed as the following unified form If time intervals t n(i,k−1,T j ) , t n(i,k−1,T j )+1 and t n(i,k,T j ) , t n(i,k,T j )+1 are in the free vibration interval, the interpolation coefficients in Equation (18) can be expressed as Else, if time intervals t n(i,k−1,T j ) , t n(i,k−1,T j )+1 and t n(i,k,T j ) , t n(i,k,T j )+1 are in the forced vibration intervals, the interpolation coefficients can be obtained through Lagrange interpolation and expressed as Substituting Equation (18) into Equation (16), and rewriting it as follows Based on Equations (15) and (21), the discrete dynamical map in matrix form on two adjacent spindle periods can be obtained as follows.
Then, the Floquet transition matrix Φ can be expressed as where the superscript '−1' indicate the operation of matrix inversion. For details of the coefficients of the state vectors in Equation (21) and matrices P, R 1 j , R 2 j , Q 1 j and Q 2 j in Equations (22) and (23), please refer to Appendix A. The stability of milling process can be determined by the spectral radius of the Floquet transition matrix Φ according to the Floquet theory.
It can be observed from the derivation process that the proposed method can adaptively complete the discrete of the spindle period and get the time step of each discrete interval only by simply setting a discretization parameter in advance. It is worth noting that under some machining parameters, the cutter will not experience free vibration, and then the variable-step method degenerates to the equal-step algorithm. Meanwhile, although the trapezoid rule is used to approximate the integral term in the algorithm derivation, other numerical methods, such as the SBM [27], the SSM [28], etc. can also be adopted to deal with Equation (15).

Algorithm Verification and Results Discussion
For convenience of comparison and verification, the milling model used in [7,13,22] is employed for numerical simulation in this section. A 19.05 mm diameter variable pitch cutter with four flutes is adopted. The helix angle is 30 • and the pitch angles are 70 • -110 • -70 • -110 • . The detailed values of the modal parameters and cutting force coefficients are as follows: the modal masses in x and y directions are m tx = 1.4986 kg and m ty = 1.199 kg, respectively; the natural angular frequencies in x and y directions are ω nx = 563.55 × 2π rad/s and ω ny = 516.27 × 2π rad/s, respectively; the damping ratios in x and y directions are ζ x = 0.0558 and ζ y = 0.025, respectively; the cutting force coefficients are K t = 6.79 × 10 8 N/m 2 and K n = 2.56 × 10 8 N/m 2 . The SLDs constructed with the AVSNIM proposed in Section 3.1 are provided to compare with those obtained with the 1st-SDM presented in [29] and equal-step numerical integration method (ESNIM) presented in [30]. The SLD is constructed over a 200 × 100-sized grid, the machining condition is down-milling and the simulation parameters are set as follows: the spindle speed Ω ranges from 2500 to 12,500 rpm and the axial depth of cut w ranges from 0 to 10 mm. The simulation is performed in Matlab ® ( where s n is the total discrete points of the spindle speed in the SLD, * means the absolute value of "*",  The radial immersion in Figures 4-6 is selected as 0.3, 0.2 and 0.1, respectively. It is found that, compared with the other two methods, the differences between the simulation results of the AVSNIM and the reference values are smaller under the same computational parameters. Consequently, the computational accuracy of the proposed method is higher than those of the other two methods. For quantitative comparison, the mean relative error (MRE) of the 1st-SDM, ESNIM and AVSNIM is provided. The MRE is defined as where n s is the total discrete points of the spindle speed in the SLD, | * | means the absolute value of "*", w a,i is the approximate value of the stability limit obtained by the numerical algorithm at the i-th discrete spindle speed, w re f ,i is the reference value of the stability limit at the i-th discrete spindle speed. Figure 7 displays the MRE of stability limits obtained by 1st-SDM, ESNIM and AVSNIM under the same simulation parameters that were used in Figures 4-6. It can be seen from Figure 7 that the stability limit obtained by the AVSNIM has the minimum MRE under each radial immersion. At the same time, it can be seen that the smaller the radial immersion is, the more obvious the advantage of the variable step size algorithm is. This shows that the proposed algorithm is particularly suitable for the stability analysis of milling process under the small radial depth of cut. This is because, under the same discretization parameter, the equal-step methods (the 1st-SDM and the ESNIM) need to discrete the whole spindle period, while the variable-step method only needs to discrete the forced vibration interval of the whole spindle period, so the time step is smaller than that of the other two methods. In numerical computation, a smaller time step usually means higher computational accuracy.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 23 stability limit obtained by the AVSNIM has the minimum MRE under each radial immersion. At the same time, it can be seen that the smaller the radial immersion is, the more obvious the advantage of the variable step size algorithm is. This shows that the proposed algorithm is particularly suitable for the stability analysis of milling process under the small radial depth of cut. This is because, under the same discretization parameter, the equal-step methods (the 1st-SDM and the ESNIM) need to discrete the whole spindle period, while the variable-step method only needs to discrete the forced vibration interval of the whole spindle period, so the time step is smaller than that of the other two methods.
In numerical computation, a smaller time step usually means higher computational accuracy. The computation efficiency of the AVSNIM is also studied and compared with that of the 1st-SDM and ESNIM. Theoretically, when obtaining the Floquet transition matrix the 1st-SDM needs to perform a large number of matrix multiplication operations. The amount of calculation will increase sharply with the increase in the approximate parameter p m , especially when p m is large. In addition, the matrix exponent needs to be calculated at each spindle speed and axial depth of cutting. While for the ESNIM and AVSNIM, the transition matrix is obtained by matrix element assignment operations, the amount of computation is very small. It should be noted that there is a difference between the AVSNIM and the ESNIM: the AVSNIM needs to calculate the variable time step under each spindle speed and axial depth of cut, while the ESNIM does not need to do so. This usually extends the computation time of the AVSNIM. Thus, we can conclude that the computation efficiency of the AVSNIM is slightly lower than that of the ESNIM but much higher than that of the 1st-SDM. The computation efficiency of the AVSNIM is validated with the two-DOF milling model. The calculation time of each method is listed in Table 1. As outlined in Table 1, the AVSNIM can save up to 75% of the calculation time compared with the 1st-SDM. According to the above analysis, it can be concluded that the AVSNIM developed in this Section has high accuracy and good calculation efficiency. The computation efficiency of the AVSNIM is also studied and compared with that of the 1st-SDM and ESNIM. Theoretically, when obtaining the Floquet transition matrix the 1st-SDM needs to perform a large number of matrix multiplication operations. The amount of calculation will increase sharply with the increase in the approximate parameter m p , especially when m p is large. In addition, the matrix exponent needs to be calculated at each spindle speed and axial depth of cutting. While for the ESNIM and AVSNIM, the transition matrix is obtained by matrix element assignment operations, the amount of computation is very small. It should be noted that there is a difference between the AVSNIM and the ESNIM: the AVSNIM needs to calculate the variable time step under each spindle speed and axial depth of cut, while the ESNIM does not need to do so. This usually extends the computation time of the AVSNIM. Thus, we can conclude that the computation efficiency of the AVSNIM is slightly lower than that of the ESNIM but much higher than that of the 1st-SDM. The computation efficiency of the AVSNIM is validated with the two-DOF milling model. The calculation time of each method is listed in Table 1. As outlined in Table 1, the AVSNIM can save up to 75% of the calculation time compared with the 1st-SDM. According to the above analysis, it can be concluded that the AVSNIM developed in this section has high accuracy and good calculation efficiency.

AVSNIM with Spindle Speed-Dependent Discretization Algorithm
The AVSNIM can provide accurate stability limits, the accuracy of which depends on the pre-set discretization parameter m p . A large discretization parameter usually brings a small time step and high calculation accuracy, and at the same time, increases the calculation time. In most of the research work, the SLD is obtained under constant discretization parameter. A constant discretization parameter will lead to different time steps and calculation accuracies under different spindle speeds. Generally, the accuracy of the stability limits in the region with low spindle speed is not as high as that in the region with a high spindle speed in the SLD. This implies that if the discretization parameter can be properly selected at each different spindle speed (usually, a large discretization parameter is used in the low-speed region and a small discretization parameter is used in the high-speed region), it is possible to improve the calculation accuracy of stability limit at a low speed and save the calculation time. Inspired by this, a spindle speed-dependent discretization algorithm is proposed to obtain the appropriate discretization parameters, and used in conjunction with the AVSNIM. The process of the algorithm is listed as follows: 1.
Sample the spindle speed range to get a series of discrete speed nodes Ω k s (Ω k s = Ω start : ∆Ω : Ω end ), where ∆Ω is the sampling interval; 2.
For each discrete speed nodes Ω k s , the stability limit w re f ,k s is obtained under a large reference discretization parameter m p through the dichotomy search algorithm adopted in [24] and used as the reference stability limit; 3.
Starting from a small discretization parameter m d , the approximate stability limit w a,k s (m d ) is obtained through the dichotomy search algorithm. The relative error between w a,k s (m d ) and w re f ,k s is calculated by the following formula

4.
Compare er k s (m d ) with the pre-set relative error limit ε, if er k s (m d ) > ε, set m d = m d + ∆m and return to step 3, until er k s (m d ) ≤ ε is satisfied. Then, record the current m d and let m(k s ) = m d ; 5.
Repeat steps 2 to 4 for all the spindle speed sampling nodes to get the discretization parameter m(k s ) at each speed node Ω k s ; 6.
Interpolate the discretization parameter m(k s ) at the spindle speed sampling points Ω k s to obtain the approximate discretization parameter m a (k) at all the spindle speed nodes on the SLD. Finally, m a (k) are rounded to obtain the final discretization parameter m f (k) using the following formula where "ceil (*)" means rounds "*" towards the nearest integer greater than or equal to "*"; 7.
The stability limit is calculated using the AVSNIM and the corresponding discretization parameter m f (k) on all the spindle speed nodes, then the SLD is obtained.
In the above algorithm, parameters ∆Ω, ε and ∆m are variable input parameters. The value of ∆Ω determines the density of spindle speed sampling points. The smaller the value of ∆Ω, the denser the sampling nodes, and the higher the interpolation accuracy of discretization parameters. The value of ε sets a relative error. The smaller the relative error is, the closer the stability limit is to the reference value. ∆m is the increment step of the discretization parameter, the value of which determines the speed of the discretization parameter approaching the reference discretization parameter. If ∆m is too small, the number of iterations will increase. On the contrary, if ∆m is too large, a discretization parameter larger than the required value may be obtained, thereby increasing the final calculation load. The values of these three parameters can be selected according to the actual engineering requirements.
Numerical simulation is performed to verify the above algorithm; two milling models are studied in this section. The first milling model is the same as used in Section 3.2. The simulation parameters are: down milling, ε = 1 × e −3 , ∆m = 2, m p = 200, the initial m d = 40, the spindle speed Ω ranges from 2500 to 12,500 rpm and the axial depth of cut w ranges from 0 to 15 mm.
Firstly, the time step under constant discretization parameter and the discretization parameters based on the spindle-speed-dependent algorithm are simulated. The results are shown in Figure 8. Figure 8a,b describes the time step under the equal-step method and the variable-step method, respectively. It can be found that the time step in the equal-step algorithm is only related to the spindle speed, while in the variable-step method, due to the effect of the helix angle, the step size depends on both the spindle speed and the axial depth of cut. Figure 8c shows the discretization parameters obtained through the spindle-speed-dependent algorithm. It can be seen that under the given relative error ε, the discretization parameters obtained by the spindle-speed-dependent algorithm are smaller than the reference discretization parameter, which means that the time required to obtain the SLD will be effectively reduced. At the same time, we can see that the value of the discretization parameter is higher in the low-spindle-speed region than that in the high-spindle-speed region in general, which proves the previous conclusion that, in the low-speed region, in order to maintain the accuracy of the numerical calculation, larger discretization parameters are required. However, the discretization parameters do not monotonically decrease with the increase in spindle speed, but there are some fluctuations. There may be two incentives for these fluctuations. One is that they may be caused by numerical calculation, because the bisection method used in this article may randomly converge to the upper or lower boundary of the relative error limit, which usually causes some numerical fluctuations. The second is that they may be caused by the different sensitivity of the eigenvalues of Floquet transition matrix to the processing parameters (mainly the spindle speed here).
Subsequently, the SLD is provided to verify the effectiveness of the proposed spindlespeed-dependent discretization algorithm. Figure 9 shows the SLD obtained through the algorithm proposed above; the SLD is constructed over a 200 × 150-sized grid; the radial immersion in Figure 9a-c is selected as 1.0, 0.6 and 0.1, respectively. The reference stability limits demoed by the red line are obtained by the AVSNIM with a large reference discretization parameter. From Figure 9a-c, we can see that the stability limits obtained by the proposed algorithm are very close to the reference stability limits. This shows that the proposed algorithm is reliable in terms of calculation accuracy. Meanwhile, since the radial immersion in Figure 9 ranges from 1.0 to 0.1, this indicates that the proposed algorithm is applicable to all kinds of radial immersion.
Finally, the calculation time required for the AVSNIM with spindle speed-dependent discretization parameters is listed in Table 2 to compare with the AVSNIM with constant discretization parameter. As outlined in Table 2, the AVSNIM with spindle-speed-dependent discretization parameters can significantly reduce the calculation time. Subsequently, the SLD is provided to verify the effectiveness of the proposed spindle-speeddependent discretization algorithm. Figure 9 shows the SLD obtained through the algorithm proposed above; the SLD is constructed over a 200 × 150-sized grid; the radial immersion in Figure  9a-c is selected as 1.0, 0.6 and 0.1, respectively. The reference stability limits demoed by the red line are obtained by the AVSNIM with a large reference discretization parameter. From Figure 9a-c, we can see that the stability limits obtained by the proposed algorithm are very close to the reference stability limits. This shows that the proposed algorithm is reliable in terms of calculation accuracy. Meanwhile, since the radial immersion in Figure 9 ranges from 1.0 to 0.1, this indicates that the proposed algorithm is applicable to all kinds of radial immersion.   Finally, the calculation time required for the AVSNIM with spindle speed-dependent discretization parameters is listed in Table 2 to compare with the AVSNIM with constant discretization parameter. As outlined in Table 2, the AVSNIM with spindle-speed-dependent discretization parameters can significantly reduce the calculation time. Another milling model to verify the proposed algorithm is from the literature [31]. In this model, a 12.75 mm diameter cutter with one flute is adopted, the helix angle of which is 30 • . The detailed values of the modal parameters and cutting force coefficients are as follows: the modal masses in x and y directions are m x = m y = 0.23 kg; the modal damping coefficients in x and y directions are c x = c y = 15.73 Ns/m; the modal stiffness in x and y directions are k x = k y = 8.4 × 10 5 N/m; the cutting force coefficients are K t = 5.36 × 10 8 N/m 2 and K n = 1.87 × 10 8 N/m 2 . The simulation parameters are: down milling, ε = 1 × e −3 , ∆m = 2, m p = 60, the initial m d = 10, the spindle speed Ω ranges from 5000 to 30,000 rpm and the axial depth of cut w ranges from 0 to 25 mm. The simulation result is shown in Figure 10. The SLD in Figure 10 Figure 10. The SLD in Figure 10 is constructed over a 500 × 250-sized grid; the radial immersion is selected as 0.02. The reference stability limits demoed by the red line are obtained by the AVSNIM with a constant discretization parameter 100 m = . As shown in Figure 10, the stability boundary is very complex, in which there are some islands. The stability limits obtained by the proposed algorithm are consistent with the reference stability limits and are the same as those obtained in the literature [31]. This shows that the algorithm proposed in this work is still effective when the stability boundary is very complex. It is worth noting that there is only a single time delay in this milling model, which is the same as the milling process with a uniform pitch milling cutter. This proves that the proposed algorithm is also suitable for the milling process with a uniform pitch milling cutter. As shown in Figure 10, the stability boundary is very complex, in which there are some islands. The stability limits obtained by the proposed algorithm are consistent with the reference stability limits and are the same as those obtained in the literature [31]. This shows that the algorithm proposed in this work is still effective when the stability boundary is very complex. It is worth noting that there is only a single time delay in this milling model, which is the same as the milling process with a uniform pitch milling cutter. This proves that the proposed algorithm is also suitable for the milling process with a uniform pitch milling cutter.
Based on the above analysis, we can conclude that the proposed algorithm can not only maintain the calculation accuracy but also effectively reduce the calculation time, and hence can meet the needs in the engineering field well.

Conclusions
This work focuses on the stability analysis of milling process with multiple delays. An AVSNIM considering the helix angle is developed. According to the geometric characteristics of the cutter (pitch helix angle, etc.) and machining parameters (spindle speed and axial cutting depth), one spindle period is divided into several segments; each segment contains a free vibration time interval and a forced vibration time interval. Then, the discretization of one spindle period is completed adaptively according to the pre-set discretization parameter. Numerical integration algorithm and Lagrange interpolation method are used to obtain the vibration state vector on each time node, and the Floquet transition matrix between two adjacent spindle periods is constructed accordingly. Finally, the stability of milling process can be analyzed on the basis of the Floquet theory. Numerical verification is conducted with the two-DOF benchmark milling model, the results of which validate that the proposed method has superior accuracy. Subsequently, a novel spindle-speed-dependent discretization algorithm is proposed to further improve the calculation efficiency of the SLD. Numerical simulation shows that this method can effectively shorten the calculation time of the SLD. The advantages of the methods proposed in this work are as follows: 1.
The AVSNIM can adaptively complete the discretization of a cutting period according to the machining parameters and the tool geometry. This is convenient for the subsequent numerical calculation; 2.
The spindle-speed-dependent discretization algorithm can get the appropriate discretization parameters flexibly according to the spindle speed and accuracy requirements, so as to effectively reduce the time required for stability analysis; 3.
The methods presented in this work are not limited by the geometry of the cutter, and can be extended to the stability analysis of the milling process with uniform pitch cutters and variable helix angle cutters.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The detailed expression of the coefficients of each state vector in Equation (21) are shown as follows. The definition of matrices P, R 1 j , R 2 j , Q 1 j and Q 2 j in Equations (22) and (23) are I −e A·∆t f r,1 I F n(1,2) F n (1,3) . . . −e A·∆t f r,2 I F n(2,2) F n (2,3) . . . F n(i,k−1) F n(i,k) . . .

−e
A·∆t f r,N d I F n(N d ,2) F n(N d , 3) . . .