An Updated Method for Stability Analysis of Milling Process with Multiple and Distributed Time Delays and Its Application

: Predicting and avoiding the onset of milling chatter are desirable to reduce its harm to machine tools, workpieces, and cutters. This paper presents an updated method to complete the stability prediction for the milling process with multiple and distributed time delays. After the dynamic of the combination milling process with variable helix cutter (VHC) and variable spindle speed (VSS) is modeled as linear delay differential equations with multiple and distributed time delays, the presented method is applied to carrying out its stability prediction for the ﬁrst time. By comparing with the existing researches and time-domain simulations, the effectiveness of the presented method has been validated. The inﬂuence and feasibility of the combination process on chatter suppression are explored and investigated for the associated one- and two-degree-of-freedom systems. Results show that the application of the combination process can realize a further suppression of milling chatter in practice. It can result in nearly 2-fold as high as the minimum depth of cut for the traditional milling or VSS milling and about 1.3-fold for VHC milling for some special domain, and can respectively lead to the average increase of stable area by 30.4%, 23.5%, and 1.5% for the adopted simulations. However, consider the contribution, the combination process is actually one process in which VHC plays an absolutely leading role but VSS plays an auxiliary role, in terms of milling stability.


Introduction
Chatter in the machining process may cause problems such as poor workpiece quality, upset noise, increased tool wear, and continuous damage to the cutting spindle. Thus, predicting and avoiding the onset of chatter are important and desirable in practice. In the past, a large number of scholars carried out the research associated with cutting parameter optimization in machining to ensure the maximum material removal rate without chatter.
It is generally believed that one of the greatest sources of cutting chatter is the regeneration effect, where the displacement of the current and previous teeth of a cutter will change dynamically, resulting in dynamic cutting thickness and dynamic milling force with phase difference. Tobias and Fishwick [1] and Tlusty and Polacek [2] were the pioneers in researching chatter mechanisms and they determined the presence of the regeneration mechanism firstly in the 1960s. The main method to predict and control regeneration chatter is to obtain a stability lobe diagram (SLD), from which one can choose the process parameters corresponding to chatter-free cutting. Consequently, following the regeneration mechanism, great deal of investigations have been done on the stability prediction methods for cutting process, trying to obtain SLD accurately and efficiently. Some typical methods are: zeroth-order frequency-domain method proposed by Altintas and Budak [3], temporal finite element method by Bayly et al. [4], multi-frequency .

y(t) = Ay(t) + B(t)y(t) −
where A is a constant matrix, N is the number of the distributed delays, and B(t), C j (t, φ) and θ j (t, φ) are periodic coefficient matrices that respectively satisfy B(t) = B(t + T), C j (t) = C j (t + T) and θ j (t, φ) = θ j (t + T, φ) with the system period T and j = 1,2, . . . , N.
Starting with the idea in Ref. [7], Equation (1) can be written in .

y(t) = Ay(t) + B(t) y(t) −
where B(t) and C j,k (t) respectively are the linear interpolation approximations of B(t) and C j,k (t) on boundary t = t i + ξ ( ξ ∈ [0, ∆t]), and they can be defined by C j,k (t) = C j,k (ξ + t i ) = C j,k,i + C j,k,i+1 − C j,k,i ∆t ξ. (5) and the delay item θ j,k,i is equal to θ j,k,i = t i+1 t i θ j,k (t)dt/∆t. When the number of intervals corresponding to θ j,k,i is defined by m j,k,i = ceil( θ j,k,i /∆t), the time-delayed item y(t − θ j,k,i ) in Equation (2) can be obtained by Lagrange quadratic interpolation [15,32] at node points y i−m j,k.i , y i−m j,k,i +1 and y i−m j,k,i +2 , i.e., y(t − θ j,k,i ) = y(ξ + t i − θ j,k,i ) = L j,k,i (ξ)y i−m j,k,i + M j,k,i (ξ)y i−m j,k,i +1 + N j,k,i (ξ)y i−m j,k,i +2 , (6) where L j,k,i (ξ) = (ξ−v j,k,i ∆t) * (ξ−v j,k,i ∆t−∆t)I 2∆t 2 , M j,k,i (ξ) = (ξ+w j,k,i ∆t) * (ξ−v j,k,i ∆t−∆t)I −∆t 2 , N j,k,i (ξ) = (ξ−v j,k,i ∆t) * (ξ+w j,k,i ∆t)I 2∆t 2 , w j,k,i = m j,k,i − τ j,k,i /∆t, v j,k,i = 1 − w j,k,i . (7) Considering y i , y i+1 and the derivatives of y(t) at t = t i and t = t i+1 , i.e., . y(t i ) = . (4) can be obtained by Cubic Hermite interpolation [32]. By collating the formula, one can get the expression Taking Equations (5), (6) and (9) into Equation (4), the solution of Equation (4) with the initial condition y i = y(t i+1 ) can be rearranged as the matrix d, d, H 1 , H 2 , e j,k,i , e j,k,i , f j,k,i , f j,k,i , g j,k,i and g j,k,i in Equation (9) following the idea in [15] can be expressed as (8) can be recast into a discrete map form [7][8][9] as where Z i is the coefficient matrix in the form of Through coupling the coefficient matrix Z i in Equation (15) for l continuous time intervals in the system period, the Floquet transition matrix [8,9] can be obtained approximately as Then, the eigenvalues of the transition matrix Φ in Equation (16) can be calculated. If all the eigenvalues are less than unity in modulus, the investigated system is stable. Otherwise, it is unstable [7].

Stability Analysis and Discussion
Based on the method presented in the Section 2, the cutting process stability about the process with VHCVSS will be analyzed. For the convenience of analysis, similar to the idea in [49], four kinds of cutting processes are defined at first: Process I (i.e., conventional cutting, see the combination of Figure 1a,c), Process II (i.e., the cutting with VSS, see the combination of Figure 1b Based on the method presented in the Section 2, the cutting process stability about the process with VHCVSS will be analyzed. For the convenience of analysis, similar to the idea in [49], four kinds of cutting processes are defined at first: Process I (i.e., conventional cutting, see the combination of Figure 1a,c), Process II (i.e., the cutting with VSS, see the combination of Figure 1b  The main purpose, here, is as follows: (1) Apply the method proposed in this paper to above four processes, carry out the stability prediction about the corresponding 1-DOF (see Figure 2a) and 2-DOF (see Figure 2b) systems, and effectively expand the application scope of the original methods; (2) The effectiveness and feasibility of chatter suppression by the combination process with VHCVSS are explored, and its influence on cutting stability is investigated. The main purpose, here, is as follows: (1) Apply the method proposed in this paper to above four processes, carry out the stability prediction about the corresponding 1-DOF (see Figure 2a) and 2-DOF (see Figure 2b) systems, and effectively expand the application scope of the original methods; (2) The effectiveness and feasibility of chatter suppression by the combination process with VHCVSS are explored, and its influence on cutting stability is investigated.

Single DOF Model of VHCVSS
For a milling process, when the workpiece is regarded to be significantly more flexible in comparison with the cutter [5], the vibration of the workpiece need to be considered but that of the cutter. Then, the associated system equations can be simplified and be expressed by the 1-DOF DDEs. This is more suitable for some specially situations, such as thin-wall cutting system or flexible support system.

Modeling
Compared with a VPC, the pitches of a VHC vary with the axial height z. Considering regeneration effect, the related time delay will be a function of z. When the effect of VSS is also taken into account, the multiple and time-varying delays ( , ) j t τ φ associated with VPCVSS (for more details, see [49]) will be changed to the multiple and distributed ones ( , ) j t z τ for VHCVSS.

Single DOF Model of VHCVSS
For a milling process, when the workpiece is regarded to be significantly more flexible in comparison with the cutter [5], the vibration of the workpiece need to be considered but that of the cutter. Then, the associated system equations can be simplified and be expressed by the 1-DOF DDEs. This is more suitable for some specially situations, such as thin-wall cutting system or flexible support system. Compared with a VPC, the pitches of a VHC vary with the axial height z. Considering regeneration effect, the related time delay will be a function of z. When the effect of VSS is also taken into account, the multiple and time-varying delays τ j (t, φ) associated with VPCVSS (for more details, see [49]) will be changed to the multiple and distributed ones τ j (t, z) for VHCVSS.
Starting with this idea, considering the nonlinear force model in [7], the 1-DOF dynamic equation of a VHCVSS milling shown in Figure 2a,c, can be expressed as m x ..
where m x ,c x = m x ω 2 x and k x = 2m x ζ x ω x respectively are the mode mass, the damping coefficient and the spring stiffness where ζ x is the damping ratios and ω x is the natural frequencies, N is the number of the cutter teeth, h j (t, z) are the dynamic chip thickness for the j-th tooth defined as where f z is the feed per tooth, q is the nonlinear factor related to the milling force, K t and K r are cutting force coefficients in tangential and radial direction, the angular position for the j tooth is where Ω(s) is the variable spindle speed and ρ is equal to tan β j /R. For a sinusoidal modulation [7] as shown in Figure 1b, the time period is T = 60/Ω 0 /RVF, the amplitude is Ω 1 = RVA × Ω 0 and its function can be expressed by where RVA = Ω 1 /Ω 0 and RVF = 60/(Ω 0 T) are respectively used to characterize the ratio of the speed variation amplitude to nominal spindle speed and the ratio of the speed variation frequency to the nominal spindle speed in a VSS milling. Considering Equations (2) and (3), and slicing the axial depth of cut into f disks axially (see Figure 3), then Equation (17) can be changed into where Appl. Sci. 2021, 11, 4203 In order to select the appropriate value of f to reduce the calculation cost, the values of h in Equation (22) can be calculated by Under the circumstances of VHCVSS, the time delay , ( ) j k t τ in Equation (20) can be written as In order to select the appropriate value of f to reduce the calculation cost, the values of h in Equation (22) can be calculated by 2πR/ tan β 1 /l. Under the circumstances of VHCVSS, the time delay τ j,k (t)in Equation (20) can be written as Clearly, one cannot obtain a closed form of τ j,k (t) in Equation (23), thus, only numerical calculation can be used [7]. It should be noted that because compared to Ω 0 , the values of Ω 1 in Equation (20) can be small enough in practice, one can expand Equation (23) in series, and obtain the approximation of τ j,k (t) (for more details, see [7,49]) as Let (21) can be represented as . where

Simulation and Analysis
Based on the 1-DOF model in Section 3.1.1 and the method in Section 2, the computer programs are written in MATLAB(R) 2009a (MathWorks, Natick, USA) and implemented on a personal computer (Lenovo, Beijing, China). The multi-parameters calculation corresponding to milling stability is carried out for above defined processes. Considering that the difference of each process is mainly reflected in the change of the variable speed parameters and the geometric parameters of variable helix cutter, therefore except for the core parameters, i.e., RVF and RVA parameters associated with VSS and the pitch at the end of the tool (ψ) and helix angle (β) parameters associated with VHC, other parameters are the same as that in [7] and remain unchanged in the subsequent calculation.
For the convenience of analysis, two related variations are defined to analyze the effect of the cutting parameters on chatter quantitatively firstly. One is the stable area ratio [49], which can be expressed by and is applied to characterize the ratio between the stable area A 0 and the total area A for a certain SLD. Another is the contribution factor κ that is to measure the contribution of a certain process to Process IV. For Processes II and III, their contributions to process IV can be respectively expressed as κ I I and κ I I I , and be defined by  [35,30,35,30] • . The material of workpieces is Al7050. The other calculation parameters are: up milling, the radial immersion ratio a D = 0.5, the feed per tooth f z = 0.1mm, the cutting-force coefficients K t = 1.07 × 10 8 N/m 1+q and K r = 0.4 × 10 8 N/m 1+q , the nonlinear parameter in cutting force q = 0.75, the modal mass m x = 3.1663 kg, the natural frequency ω x = 400 Hz and the damping ratio ξ x = 0.02. The black thin lines in Figure 4 are associated with Process I, they are the reference stability boundaries of the traditional milling with Cutter-1; the red dotted lines correspond to process III with Cutter-2; the cyan thin lines in the subgraphs (a)-(d) of Figure 4 respectively are the stability boundary of Process II with different RVF (RVF = 0.05, 0.1, 0.2, and 0.5) but the same RVA (RVA = 0.1) when Cutter-1 is used; the SLDs shown by the green thick lines are related to Process IV, who uses Cutter-2 on the basis of above Processes II. Obviously, it can be found from Figure 4 that compared to Process I, Processes II and III do lead to better milling stability gains, especially for low-speed region. Also, their corresponding values of χ are always larger (see Figure 5). The phenomena above imply the active influence of VHC or VSS on milling stability by the disruption of regeneration chatter, and they are also coincident with that in [35][36][37][38][40][41][42]. Then, Process IV is focused on. As seen from the green lines in Figure 4, Process IV almost always leads to the highest axial depth of cut for every adopted RVF. For instance, the values of absolute (minimum) depth of cut for Ω = 3600 rpm in Figure 4b-d, can reach up to 5.4, 5.4, 5.5, and 6.5 mm, which are nearly 2-fold as high as these for processes I and II and are about 1.3-fold for Process III. The values of χ associated with SLDs in Figure 4 are compared and shown in Figure 5. It can be seen that Process IV could also gain the biggest χ. From Process I to II to III to IV, the obtained values of χ keep increasing gradually. Compared with Processes I, II, and III, the average increases of stability region for Process IV at different RVF values are about 30.4%, 23.5%, and 1.5%, respectively. This illustrates that the application of Process IV can realize a further suppression of milling chatter.
However, it should be noted that compared with the stable region difference between Processes IV and III, the difference between Processes IV and II is more obvious. This seems indicate that although Processes II and III both play a contribution for realizing the combination effect, the contribution degree is obviously different. Figure 6 shows the contribution factors of Processes II and III to the combination effect for different RVF values, based on Equation (27). It can be seen that although the values of κ for Process II is almost larger with the increase of RVF, Process III obviously has a greater contribution (the values of κ I I I are all greater than 85%). This seems indicate that Process III plays a more important role in affecting the cutting stability of Process IV. In other words, VHCVSS is essentially a process in which VHC plays an absolutely leading role but VSS plays an auxiliary role, in terms of cutting stability.

Model Verification
Here, the presented method will be verified by comparing with the existing research and time-domain simulation. In order to evaluate the effectiveness of this method, the researches associated with VSS process [7] and with the combination process of VPC and VSS [50] are taken into accout at first. If the comparisons of the stability boundaries plotted by the black lines in Figure 4 with those in Figure 5.17 in [7], and the comparisons of the stability boundaries plotted by the cyan lines in Figure 4 with those in Figure 3 in [50] are respectively carried out, one can see that the results gained by the presented method coincide well with those by the methods in [7,50]. This illustrates the applicability of the presented method for the stability prediction about Processes I and II.
For Process IV, to the knowledge of the authors, there is no reference literature at present. Thus, the time-domain simulation is mainly used to carry out the verification. In [53], authors constructed the contour plot of peak-to-peak (PTP) vibration by investigating the vibration amplitude of a cutter. In view of this, based on the method in [53], an improved time-domain program suitable for VHCVSS is proposed. It should be noted that compared with the program in [53], the one in this paper can not only consider the infulence of regeneration effect, cutter runout, multimodal of the system, loss of contact, and variable pitch/helix angles in the original method, but also that of variable spindle speed.
The top diagram of Figure 7 shows the SLDs in Figure 4d between the rotation speed Ω = 1500-5000 rpm and the axial depth of cut p a = 0-14mm. The Subgraphs (a)-(d) in Figure 7 respectively show the corresponding contour plot of PTP by the time-domain method for above four processes under the same parameters. It can be seen that the stability boundaries obtained through the presented method are basically consistent with the that through time-domain method, not only for their curve shapes but also for the stability evolution trends caused by the processes differences. This implies the effectiveness of this method in predicting stability for VHCVSS milling. There are little differences between the predicting results and simulation ones. This may attribute to the fact that the presented method can not take into account the influence of loss-of-contact character on cutting stability, but the time domain method can.

Model Verification
Here, the presented method will be verified by comparing with the existing research and time-domain simulation. In order to evaluate the effectiveness of this method, the researches associated with VSS process [7] and with the combination process of VPC and VSS [50] are taken into accout at first. If the comparisons of the stability boundaries plotted by the black lines in Figure 4 with those in Figure 5.17 in [7], and the comparisons of the stability boundaries plotted by the cyan lines in Figure 4 with those in Figure 3 in [50] are respectively carried out, one can see that the results gained by the presented method coincide well with those by the methods in [7,50]. This illustrates the applicability of the presented method for the stability prediction about Processes I and II.
For Process IV, to the knowledge of the authors, there is no reference literature at present. Thus, the time-domain simulation is mainly used to carry out the verification. In [53], authors constructed the contour plot of peak-to-peak (PTP) vibration by investigating the vibration amplitude of a cutter. In view of this, based on the method in [53], an improved time-domain program suitable for VHCVSS is proposed. It should be noted that compared with the program in [53], the one in this paper can not only consider the infulence of regeneration effect, cutter runout, multimodal of the system, loss of contact, and variable pitch/helix angles in the original method, but also that of variable spindle speed.
The top diagram of Figure 7 shows the SLDs in Figure 4d between the rotation speed Ω = 1500-5000 rpm and the axial depth of cut a p = 0-14mm. The Subgraphs (a)-(d) in Figure 7 respectively show the corresponding contour plot of PTP by the time-domain method for above four processes under the same parameters. It can be seen that the stability boundaries obtained through the presented method are basically consistent with the that through time-domain method, not only for their curve shapes but also for the stability evolution trends caused by the processes differences. This implies the effectiveness of this method in predicting stability for VHCVSS milling. There are little differences between the predicting results and simulation ones. This may attribute to the fact that the presented method can not take into account the influence of loss-of-contact character on cutting stability, but the time domain method can.

2-DOF Model of VHCVSS
For a milling process, when the contained cutter is regarded as a cantilever beam but the workpiece is assumed to be absolute rigidity, the vibration of the tool has to be considered but that of the workpiece. Then, the associated system equation can be simplified and defined by the DDE with 2-DOF. Compared with the 1-DOF system in Section 3.1, the 2-DOF one is a more universal and general expression, so more attentions are paid on it for the existing research.

Modeling
Mathematically, the 2-DOF model associated with VHCVSS milling process, shown in subgraphs (b) and (c) of Figure 2 can be written as , , where the cutting force coefficients , ( , )

2-DOF Model of VHCVSS
For a milling process, when the contained cutter is regarded as a cantilever beam but the workpiece is assumed to be absolute rigidity, the vibration of the tool has to be considered but that of the workpiece. Then, the associated system equation can be simplified and defined by the DDE with 2-DOF. Compared with the 1-DOF system in Section 3.1, the 2-DOF one is a more universal and general expression, so more attentions are paid on it for the existing research.

Modeling
Mathematically, the 2-DOF model associated with VHCVSS milling process, shown in subgraphs (b) and (c) of Figure 2 can be written as ..
if the axial depth of cut is sliced, considering Equations (2), (3) and (28) will be changed into M .. where Denoting p(t) = M . q + cq/2 and x(t) = q(t) p(t) T , one can transform the 2-DOF model in Equation (28) into where

Simulation and Analysis
Here, multi-parameters calculation for the milling stability of Processes I, II, III, and IV, is carried out based on the two DOF model in Section 3.2.1 and the method in Section 2. Similar to the situation in Section 3.1.2, the other parameters in the subsequent calculation remain unchanged except for the core parameters RVF, RVA, ψ, and β. The mainly parameters are selected from the representative papers [54], and are: down milling, half immersion, q = 1, K r = 2.56 × 10 8 N/m, K t = 6.79 × 10 8 N/m, N = 4, f z = 0.1 mm, ω x = 563.6 Hz, ω y = 516.21 Hz, ζ x = 0.0558, ζ y = 0.025, m x = 1.4986 kg, and m y = 1.199 kg. Figure 8 shows some SLDs corresponding to the 2-DOF systems for Processes I, II, III, and IV. The black thin lines are the stability boundaries corresponding to Process I for reference, whose core parameters are:  Figure 8, it can be found that the combined process (i.e., Process IV) still shows a better cutting stability suppression advantage than that of Processes I, II, and III. For different combinations of variable helix angles β, the values of χ associated with Process IV are the biggest (see Figure 9). For Process IV, the values of χ will increase for a bigger helix angle increment ∆β (defined by β 1 − β 2 ). This shows that the increasing of ∆β plays a positive role in enhancing the stability area for VHCVSS. Thus, a bigger ∆β should be suggested in practice. [30,30,30,30]°. The cyan thin lines are related to the stability boundaries of Process II with RVF = 0.5, RVA = 0.1 and the same ψ and β in Process I. The boundaries plotted by red dotted lines in Figure 8a-d are associated with Process III, where four VHCs with the different β ( β = [30,30,30,30]°, β = [32,30,32,30]°, β = [34,30,34,30]°, and β = [36,30,36,30]°) but the same ψ (ψ = [75, 85, 95, 105]°) are utilized. The green thick lines exhibit the SLDs related to Process IV for the relevant parameters combinations of Process II and Process III. Comparing with the curves in the Figure 8, it can be found that the combined process (i.e., Process IV) still shows a better cutting stability suppression advantage than that of Processes I, II, and III. For different combinations of variable helix angles β , the values of χ associated with Process IV are the biggest (see Figure 9). For Process IV, the values of χ will increase for a bigger helix angle increment β Δ (defined by 1 2 β β − ). This shows that the increasing of β Δ plays a positive role in enhancing the stability area for VHCVSS. Thus, a bigger β Δ should be suggested in practice.  Figure 10 shows the contribution factor diagram of Processes II and III to Process IV under different β . Obviously, the corresponding values of III κ are all greater than 61%, Process III still has a significant advantage in the contribution of the combination effect, compared to Process II. However, it should be noted that the contribution factor changed most significantly when the variable helix angles change from β = [30,30,30,30]° to β compared to Process II. However, it should be noted that the contribution factor changed most significantly when the variable helix angles change from β = [30,30,30,30]° to β = [32,30,32,30]°. This is mainly due to the fact that the former actually constructs a VPC milling process that is consistent with the model described in [49], while the latter is a VHC one. In comparison, variable helix characteristics of a VHC aggravate the degree of regeneration disturbance, and thus relatively weaken the stability influence ability of VSS in the combination process.   Figure 10 shows the contribution factor diagram of Processes II and III to Process IV under different β. Obviously, the corresponding values of κ I I I are all greater than 61%, Process III still has a significant advantage in the contribution of the combination effect, compared to Process II. However, it should be noted that the contribution factor changed most significantly when the variable helix angles change from β = [30,30,30,30] • to β = [32,30,32,30] • . This is mainly due to the fact that the former actually constructs a VPC milling process that is consistent with the model described in [49], while the latter is a VHC one. In comparison, variable helix characteristics of a VHC aggravate the degree of regeneration disturbance, and thus relatively weaken the stability influence ability of VSS in the combination process. regeneration disturbance, and thus relatively weaken the stability influence ability of VSS in the combination process.

Conclusions
This paper presents an improved method to gain the SLDs for milling process with multiple and distributed time delays. After the milling dynamic of the VHCVSS process, taking into account regeneration effect, is modeled by linear DDEs with multiple and distributed time delays, the presented method is applied to deal with the problems about the stability prediction for VHCVSS process for the first time. By comparing with the existing research and time-domain simulations, the presented method has been validated. The effectiveness and feasibility of chatter suppression by the combination of VHC and VSS milling are explored and investigated for associated 1-and 2-DOF systems. Moreover, the following conclusion can be drawn: (1). The presented method is effective for obtaining the SLDs for milling processes with multiple and distributed time delays. Thus, the theory in [15,32] is developed and its application scope is widened.

Conclusions
This paper presents an improved method to gain the SLDs for milling process with multiple and distributed time delays. After the milling dynamic of the VHCVSS process, taking into account regeneration effect, is modeled by linear DDEs with multiple and distributed time delays, the presented method is applied to deal with the problems about the stability prediction for VHCVSS process for the first time. By comparing with the existing research and time-domain simulations, the presented method has been validated. The effectiveness and feasibility of chatter suppression by the combination of VHC and VSS milling are explored and investigated for associated 1-and 2-DOF systems. Moreover, the following conclusion can be drawn: (1) The presented method is effective for obtaining the SLDs for milling processes with multiple and distributed time delays. Thus, the theory in [15,32] is developed and its application scope is widened. (2) The presented method is applied to dealing with the stability prediction problem with regard to the new process, i.e., VHCVSS for the first time. Obviously, the method is not only applicable to the cases associated with VPC, VHC, or VSS, but also to their combination cases, e.g., VPCVSS and VHCVSS. (3) Compared with Processes I, II, and III, application of Process IV can result in nearly 2-fold as high as the minimum depth of cut for Processes I and II and are about 1.3-fold for Process III for some special domain, and can respectively lead to the average increases of stable area by 30.4%, 23.5%, and 1.5% at different RVF values. This indicates that one can improve the chatter suppression performance of the cutting system through VHCVSS in practice. The proposed method can also provide a good help for the subsequent process optimization to achieve the chatter-free maximum removal rate.
(4) For Process IV, the increasing of the helix angle increment ∆β plays a positive role in enhancing the stability area for VHCVSS. Thus, bigger values of ∆β should be suggested in practice. (5) The contribution of Process III to the cutting stability of Process IV is more than 85% and 61% for the associated 1-and 2-DOF systems, respectively. Thus, essentially, VHCVSS is a process in which VHC plays an absolutely leading role for milling stability but VSS plays an auxiliary role.

Conflicts of Interest:
The authors declare no conflict of interest.
Nomenclature a p axial depth of cut (m) a D radial immersion ratio c x ,c y damping coefficient in x and y directions CCM Chebyshev collocation method f number of slicing the axial depth of cut f z feed per tooth (mm/tooth) FDM full-discretization method g functions which determine whether the tooth is in contact with workpiece h j dynamic chip thickness for the j-th tooth HPM Homotopy perturbation method l number of discrete time intervals k x ,k y spring stiffness in x and y directions K t ,K r cutting force coefficients in tangential and radial direction (N/m) m j,k,i approximate number of discrete parts interval corresponding to τ j,k,i m x ,m y modal mass in x and y directions (kg) N number of teeth PTP peak-to-peak q nonlinear parameter in cutting force R radius of a cutter (m) RVA ratio of the speed variation amplitude to the nominal spindle speed RVF ratio of the speed variation frequency to the nominal spindle speed SDM semi-discretization method SLD stability lobe diagram VHC variable helix cutter VHCVSS combination of VHC and VSS VPC variable pith cutter VSS variable spindle speed β helix angles ( • ) ∆β helix angle increment ( • ) κ contribution factor ζ x ,ζ y damping ratios in x and y directions τ j,k,i approximate delay time φ j angular position for the j-th tooth ( • ) χ stable area ratio ψ pitches at the end of the tool ( • ) ω x ,ω y the natural frequencies in x and y directions (Hz) Ω(s) spindle speed of VSS (Rev/min) Ω 0 foundation speed of VSS (Rev/min) Ω 1 amplitude variation of spindle speed (Rev/min)