# Stability Analysis of Milling Process with Variable Spindle Speed and Pitch Angle considering Helix Angle and Process Phase Difference

Authors: Gang Jin; Haotian Jiang; Jianxin Han; Zhanjie Li; Hua Li; Bing Yan
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Suppression of milling chatter by disrupting regenerative effect is a well-known method to obtain higher cutting stability domain. In this paper, a dynamic model of the milling process with variable spindle speed and pitch angle considering helix angle and process phase difference is presented. Then, an updated semidiscretization method is applied to obtain the stability chart. After the effectiveness of the proposed method is confirmed by comparisons with the previously published works and the time-domain simulations, lots of analyses are conducted to deeply evaluate the influence of the helix angle, the process phase difference, and feed per tooth on milling stability. Results show that the change of helix angle can result in significant stability discrepancies for both high-speed and low-speed regions. Though the process phase difference has the randomness and immeasurability in the practical application, it has an important influence on the stability and will result in a periodic evolution of the stability with a periodπ. Also, its recommended values are given for the practical milling process.


---

## Body

## 1. Introduction

Chatter in metal cutting processes is a form of self-excited vibration and always causes lower material removal rates, machining quality, poor accuracy and surface finish, unpleasant noise and sound, accelerated tool wear, and so on in practice. The most powerful source of chatter is regeneration, which is associated with the surface waviness variation during successive cuts. Thus, in order to avoid chatter and achieve high productivity in machining, lots of investigations had been done based on the regeneration mechanism for optimizing the machining conditions in the past few decades.Generally, the regenerative chatter is modeled and analyzed in the frequency domain or time domain. The most established method for predicting and preventing regeneration chatter is of cutting parameters selection from a stability lobe diagram (SLD), which can predict the stability boundary and distinguish the stable and unstable cutting region. Around this theme, various typical methods, such as zeroth-order approximation method [1], temporal finite element analysis [2], multifrequency solution [3], Chebyshev polynomial approximation [4], semidiscretization method (SDM) [5–7], full-discretization method (FDM) [8, 9], spectral element method [10], have been proposed for more rapidly and accurately predicting the stability chart successively.Alternatively, some methods focused on enlarging the stable cutting region of the stability chart by disrupting the regenerative effect. One typical example is the utilization of some special cutters, such as variable pitch cutter (VPC) suggested by Slavicek [11]. In this case, the modulation in chip thickness caused by a nonuniform pitch can reduce the vibration under disturbance. Aiming at predicting the stability of milling process with VPC, frequency-domain model [12, 13], the Cluster Treatment of Characteristic Roots [14], unified SDM [15], improved FDM [16, 17], and improved SDM [18] have been proposed successively. Furthermore, Sims et al. [19], Sims [20], Dombovari and Stepan [21], and Jin et al. [22, 23] proposed efficient methods to predict the stability for VPC and variable helix cutter, while Yusoff and Sims [24] optimized the geometry of variable helix tool to suppress regenerative chatter. Besides, the use of serrated cutters [25] also exhibited good capacity and potential for a higher stable depth of cuts.Another typical example to achieve regenerative chatter disruption is the use of variable spindle speed (VSS), which has a similar mechanism to the use of cutter with variable pitch/helix cutter. Insperger and Stepan [26] and Zatarain et al. [27], respectively, used the SDM and frequency-domain method to construct stability diagrams in the VSS process. Seguy et al. [28] studied VSS in high-speed milling and showed its effective suppression on the unstable domain of the first flip lobe. Long and Balachandran [29] presented a stability treatment for upmilling and downmilling processes with a VSS and discussed the benefits of VSS milling operations. Later, improved SDM [30] and Chebyshev collocation method [31] are presented to efficiently evaluate the process stability in milling with VSS, respectively.Recently, Jin et al. [32] constructed a new dynamic model to investigate the potential of chatter suppression for the milling process in which the variable pitch cutter and variable spindle speed (VPCVSS) are considered simultaneously. Results show that the combined effect of VPC and VSS on improving the stability region exists and it is more remarkable in certain cases. Obviously, VPCVSS is not a simple combination of VPC and VSS. For modeling, the special characteristics caused by the processes combination, such as time-varying delay, process phase difference, and time-varying feed, significantly increase the difficulty of dynamic modeling, solving, and simulation. For stability impact, the form of the perturbation of the regeneration has realized upgrading to some extent. This upgrade is not a simple superposition of the positive effects of VPC and VSS on system stability but has its own mechanism and phenomenon. However, the research still exists following deficiencies:(1)
The actual feature of a tool, helix angle, is not considered during modeling. However, for a general cutting process, the effect of helix angle is crucially important [33, 34]. Thus, the related results cannot reflect the influence of the helix angle on the stability, and its real influence is still worthy of exploration and investigation.(2)
The analysis is based on an assumption that the value of the process phase difference is equal to 0; however, it is nonzero usually (note that this will be discussed in detail in SectionSection 2.1.1). Thus, its influence on system stability needs to be explored further.(3)
The analysis focuses on the influence of the process parameters, that is, the parameters associated with variable tooth and variable speed, on the system stability. However, other important cutting conditions (such as feed rate, feed direction, and helix angle) are not studied yet.With the above questions in mind, in this paper, a dynamical model of milling process of VPCVSS considering the process phase difference and the helix angle is established for the first time, and then the influence analyses associated with process phase angle, helix angle, and feed per tooth are conducted based on an updated SDM; some relevant engineering suggestions are given. The structure of the paper is as follows: inSection 2, the mathematical model and the method for milling stability are introduced. In Section 3, the method verification and the stability analysis are conducted. In Section 4, conclusions are presented in brief.


## 2. Stability Prediction of Milling with VPCVSS

Generally, the dynamics of the tool-workpiece system in the milling process is described by the models in the form of linear delay-differential equations considering regenerative effect. Based on the number of system modes, these equations usually are single DOF or two DOFs ones. For the former, the workpiece is assumed to be more flexible compared to the rigid tool [5, 35, 36], but the latter models the tool as a cantilever beam because of the consideration of rigid workpiece [1–3].Here, for the convenience of the following analysis, the single DOF model introduced in [5] is chosen and its schematic diagram is shown in Figure 1(a), where a helix tool with radius R, helix angle β, and N unequally spaced teeth is rotating at a variable velocity Ωt in rpm, but the workpiece is excited by the milling force at the tooltip. The mathematical representation of this system is given by (1)mxx¨t+cxx˙t+kxxt=−Fxt,where xt and Fxt are the position and cutting force in the x direction. mx, cx, and kx are the modal mass, damping coefficient, and spring stiffness, respectively.Figure 1
Schematic mechanical model of a single DOF milling system (a), regenerative effect (b), and different regions of cutting caused by the helix angle ((c) and (d)).


### 2.1. Cutting ForceFxt

Because of the existence of the cutter helix angle, the height of the cutter teeth in the cutting changes all the time. Thus, the cutting force acting on the whole axial cutting height is the integration of differential force with respect to the differential axial depthdz along the compliant direction of the workpiece, as shown in Figure 1(c). If the nonlinear force model suggested in [32] is considered, the cutting force here can be expressed as(2)Fxt=∑j=1N∫0apdFjxt,z=∑j=1N∫0apgϕjsinϕjkt+cosϕjkrhjqtdz,where q is the nonlinear force coefficient. Influenced by the judgment function gϕj in equation (2), Fxt is piecewise continuous [37]. As shown in Figure 1(d), each helix tooth will go through three regions, that is, entry (labeled by ①), middle of the cut (labeled by ②), and exit region (labeled by ③) in a milling process. For one case with the same helix angle and axial depth of cut, there may be two types of such regions due to the different radial depth of cut, as shown in the left and middle graphics of Figure 1(d). Obviously, each region corresponds to different calculation methods of the upper and lower limits of integration in equation (2). Here, function Gϕj is introduced instead of the original gϕj, and the unified calculation models for the upper and lower limits of integration; that is, ϕejt=minϕjt,0,ϕe and ϕsjt=maxϕjt,ap,ϕs are utilized (for more details, see [17, 18]). Then, equation (2) can be transformed into(3)Fxt=γ∑j=1NGϕj∫ϕsjtϕejtsinϕjkt+cosϕjkrhjqtdϕ,where γ=R/tanβ and Gϕjt is the new judgment function which determines whether the tooth is in contact with the workpiece, defined by(4)Gϕjt=1,ifϕs<ϕjt<ϕe+ϕhex,0,otherwise,where ϕs and ϕe, respectively, are the exit and start angles which depend on radial immersion ratios aD and milling style (up- or downmilling), ϕhex is the exit angle difference caused by the helix angle and is equal to ϕhex=aptanβ/R.


#### 2.1.1. Angular Positionϕjt

Generally, angular positionϕjt is the function of helix angle and spindle speed. Here, it is also relevant to the tooth pitch angle. If a linear tooth pitch variation is shown as(5)ψ=ψ1,ψ2,…,ψN−1,ψN=ψ0,ψ0+Δψ,…,ψ0+N−1Δψ,where ψ0=360/N−N−1Δψ/2, and a sinusoidal spindle speed modulation modeled as(6)Ωt=Ω0+Ω1sin2πTmt=Ω01+RVA⋅sinRVF⋅2π60Ω0t,where RVA=Ω1/Ω0 is the ratio of the speed variation amplitude to the nominal spindle speed and RVF=60/Ω0Tm is the ratio of the speed variation frequency to the nominal spindle speed, are considered the same as [32], the angular position ϕjt can be defined by(7)ϕjt=2π60∫0tΩsds−zγ,ifj=1,2π60∫0tΩsds+∑i=1j−12πψi360−zγ,if1<j≤N.As a spindle rotates in the machining process, the pitch angle of VPC changes periodically at the rotation periodτ0=60/Ω0, whereas the spindle speed varies with the modulation period Tm=60/Ω0/RVF. Therefore, when VSS and VPC are combined, there must be some problems associated with the phase. As shown in Figure 2, one VPC is rotating with Ωt, when its 1-th tooth passes through the positive direction of Y-axis, that is, the red lines; the value of its corresponding instantaneous spindle speed is not equal to Ω0, but the value is between intervals [Ω0−Ω1,Ω0+Ω1], which is depending on the phase of Ωt at this moment. In order to consider this issue, equation (6) is rewritten in a new form including the initial phase as(8)Ωt,η=Ω0+Ω1sin2πTmt+η,where η is defined as the process phase difference between VPC and VSS. Obviously, there are kinds of possible η in VPCVSS milling process, such as π/4, π, and 5π/4. In other words, the value of η is a random one between 0 and 2π in the practical application. As a result, the times to pass through any pitch angle ψj are different when the teeth is rotating at Ωt,η, for example, Cases (1) and (2) shown in Figure 2. Considering the mechanism of regenerative chatter, it means that η will affect the time-varying characteristics of time delays in the milling system. Further, η does make sense for milling stability.Figure 2
Schematic diagram of process phase difference for VPCVSS system.


#### 2.1.2. Instantaneous Chip Thicknesshjt

For a milling process, the instantaneous chip thicknesshjt generally consists of two parts, that is, the static chip thickness hjstt mainly determined by the feed per tooth and angular position, and the dynamic displacements vjt and vj−1t, respectively, corresponding to the present tooth j and previous tooth j-1 in the radial direction [1]. Thus, hjt can be expressed as(9)hjt=hjstt,η+vj−1t−vjt=ft,jt,ηsinϕjt−yt−τjt,ηsinϕjt+ytsinϕjt,whereτjt,η represents the time delay between the teeth j and j−1, which is mainly affected by the parameters of VPCVSS, such as spindle speed and pitch angles, and ft,jt,η represents the feed per tooth, which can be described by(10)ft,jt,η=vfτjt,η=fmϑjt,η,where ϑjt,η=Nψj1−RVAsinωmt+η/360, which is applied to represent the relationship between nominal feed per tooth fm and ft,jt,η.


#### 2.1.3. Time Delaysτjt,η

Generally, time delay is determined by calculating the rotation time of the tool over the pitch between two successive teeth except for some special cases, such as big cutter runout [15, 37]. For VPCVSS milling, the time delays can be calculated, based on the idea in [5]. Thus, an implicit form including τjt,η is conducted as(11)∫t−τjt,ηtΩ0+Ω1sinωms+η60ds=ψj360.Substituting equation (6) into equation (11) and rearranging both sides of equation (11), then, τjt,η can be approximately expressed by the following relationship:(12)τjt,η≈ψj/κ6,where κ=Ω0+Ω1sinωmt+η. If Ω1 is small enough compared with Ω0, that is, RVA is small, equation (12) can be expanded and then approximated in the form of(13)τjt,η≈τ0,j−τ1,jsinωmt+η,where τ0,j=ψj/Ω0/6 and τ1,j=τ0,jRVA. Here, it should be noted that the approximation errors of equation (13) increase with the increase of RVA. For small RVA, for example, RVA = 0.2, there are also errors although the effects of errors on the stability lob may be ignored.Substituting equations (7) and (9) into equation (3) and considering equation (1), then linearizing the consequent equation on the basis of the ideal in [5], a linear time-periodic DDE can be given by(14)mxx¨t+cxx˙t+kxxt=−ε∑j=1NKjtxt−xt−τjt,η,with(15)ε=γq,Kjt=ft,jq−1t,ηGϕj∫ϕsjtϕejtktsinqϕjtcosϕjt+krsinq+1ϕjtdϕ,


### 2.2. Updated SDM

SDM is a common method often applied to the stability analysis of various DDEs where the time delayed terms are approximated by a piecewise constant function while the current time terms are left unchanged. Here, the method is presented for the case of VPCVSS.However, to ensure that equation (14) is periodic and applicable to SDM, the ratio of the spindle modulation period Tm and the rotation period τ0 in the current system must be a rational number, that is, q1Tm=q2τ0 where q1 and q2 are relative prime numbers [5]. Thus, the system period Tc can be expressed as(16)Tc=q1Tmorq2τ0.By introducing the system vectorUt=xt,x˙tT and using Cauchy transformation, the second-order DDE, that is, equation (14) can be written as a system of the first-order DDE:(17)U˙t=ΑtUt+∑j=1NBjtUt−τjt,η,At=At+Tc,Bjt=Bjt+Tc,with(18)At=01−ωnx2−εltmx−2ζxωnx,Bjt=00ε∑j=1Nkjtmx0,lt=∑j=1Nkjt,where ωnx=kx/mx and ζx=cx/2ωnxmx are the angular natural frequency and the damping ratio.The first step of SDM is the construction of the discrete time intervalti,ti+1 by dividing its period Tc into k number with length Δt, that is, Δt=Tc/k. Then, the time-varying delay τjt,η in equation (17) for the interval ti,ti+1 can be written in an average form based on the ideal in [5] as(19)τi,j=1Δt∫titi+1τjt,ηdt=τ0,j−τ1,jci,where(20)ci=kν2π∫i2π/k/ν+ηi+12π/k/ν+ηsintdt,i=0,1,…,k−1ν=TmTc.Then, substituting equation (19) into equation (17), equation (17) can be approximated in [ti, ti+1] as new form yields as(21)U˙t=AiUt+∑j=1NBi,jαi,jUi−mi,j+1+χi,jUi−mi,j=AiUt+Wi,with(22)Ai=1Δt∫titi+1Atdt,Bi,j=1Δt∫titi+1Bjtdt,αi,j=mi,jΔt+Δt/2−τi,jΔt,χi,j=1−αi,j,Ui−mi,j+1=Uti−mi,j+1,Ui−mi,j=Uti−mi,j,mi,j=intτi,j/Δt+1/2,where int (∗) indicates the operation that rounds positive number towards zero. For the initial value U(ti) = Ui, the solution of equation (21) is known: (23)Ut=eAit−tiUi+Ai−1Wi−Ai−1Wi.Substitutingt = ti+1 and Uti+1=Ui+1 into the solution leads to(24)Ui+1=QiUi+∑j=1NMi,j,αUi−mi,j+1+Mi,j,χUi−mi,j,where Qi=eAiΔt,Mi,j,α=eAiΔt−IΑi−1Bi,jαi,j,Mi,j,χ=eAiΔt−IAi−1Bi,jχi,j.Defining the maximum value ofmi,j by n=maxmi,j and state augmentation of equation (24) results in the 2(n+1) dimensional discrete map: (25)Vi+1=ZiVi,with(26)Vi=colUi,Ui−1,…,Ui−n+1,Ui−n,(27)Zi=Qi0⋯000I0⋯0000I⋯000⋮⋮⋱⋮⋮⋮00⋯I0000⋯0I0+∑j=1N0⋯Mi,j,αMi,j,χ⋯00⋯00⋯00⋯00⋯0⋮⋮⋱⋮⋮⋮0⋯00⋯00⋯00⋯0.The horizontal position of the discrete input matricesMi,j,α and Mi,j,χ in the right side of equation (27) begins from the column of 2mi,j−1 and 2mi,j+1, respectively. Here, it should be noted that if one reduction about the matrix Vi is carried out by removing the delayed values of the velocities, like [7], the size of vector Vi can be decreased to n+2, the horizontal position of Mi,j,αand Mi,j,χ will change from the column of mi,j+1 and mi,j+2, respectively.The Floquet transition matrix over the principal periodTc is approximated by coupling the solutions of k successive time intervals Δt:(28)Φ=Zk−1Zk−2,…,Z1Z0.Stability of the investigated system is determined by the eigenvalues of the transition matrixΦ. The system is stable if all eigenvalues of Φ are in modulus less than 1.The main advantage associated with the method in this paper can be summarized as follows:(1)
Compared with the original algorithm [5, 6, 36], the method in this paper can be applied not only to variable pitch milling cutter or variable speed milling process but also to their composite systems. Mathematically, it is suitable for the stability solution of DDEs with multiple variable delays.(2)
Due to the addition of equations (16) and (20), the method in this paper is also applicable to some special cases, where number 1 cannot be divisible by the values of RVF, such as RVF = 0.22, 0.3, or 0.4. For the above cases, this method is applicable, but the method in [32] is not.In addition, it should be noted that the 1-DOF VPCVSS model is proposed in this paper; however, the milling system belongs to the 1-DOF model in general. Thus, to realize a conversion from 1-DOF one to 2-DOF one [38], some modifications need to be carried out as follows: (a) the system equation (equation (1)) should consider the dynamics in both X and Y directions instead of the original X. (b) Correspondingly, the cutting force model in X and Y directions (equations (2) and (3)) should be constructed. (c) The dynamic cutting force thickness (equation (9)) is changed to a coupled model considering X and Y directions. (d) The system equation (equation (17)) and the corresponding state matrix are changed from two dimensions to four ones. (e) Equations (25) and (26) change from 2(n+1) to 4(n+1) or 2n+4 if the delayed values of the velocities are removed, thus, the horizontal position of the discrete input matrices Mi,j,α and Mi,j,χ in equation (27), respectively, begins from the column of 4mj,k,i−3 and 4mj,k,i+1or 2mj,k,i+1 and 2mj,k,i+3.


## 2.1. Cutting ForceFxt

Because of the existence of the cutter helix angle, the height of the cutter teeth in the cutting changes all the time. Thus, the cutting force acting on the whole axial cutting height is the integration of differential force with respect to the differential axial depthdz along the compliant direction of the workpiece, as shown in Figure 1(c). If the nonlinear force model suggested in [32] is considered, the cutting force here can be expressed as(2)Fxt=∑j=1N∫0apdFjxt,z=∑j=1N∫0apgϕjsinϕjkt+cosϕjkrhjqtdz,where q is the nonlinear force coefficient. Influenced by the judgment function gϕj in equation (2), Fxt is piecewise continuous [37]. As shown in Figure 1(d), each helix tooth will go through three regions, that is, entry (labeled by ①), middle of the cut (labeled by ②), and exit region (labeled by ③) in a milling process. For one case with the same helix angle and axial depth of cut, there may be two types of such regions due to the different radial depth of cut, as shown in the left and middle graphics of Figure 1(d). Obviously, each region corresponds to different calculation methods of the upper and lower limits of integration in equation (2). Here, function Gϕj is introduced instead of the original gϕj, and the unified calculation models for the upper and lower limits of integration; that is, ϕejt=minϕjt,0,ϕe and ϕsjt=maxϕjt,ap,ϕs are utilized (for more details, see [17, 18]). Then, equation (2) can be transformed into(3)Fxt=γ∑j=1NGϕj∫ϕsjtϕejtsinϕjkt+cosϕjkrhjqtdϕ,where γ=R/tanβ and Gϕjt is the new judgment function which determines whether the tooth is in contact with the workpiece, defined by(4)Gϕjt=1,ifϕs<ϕjt<ϕe+ϕhex,0,otherwise,where ϕs and ϕe, respectively, are the exit and start angles which depend on radial immersion ratios aD and milling style (up- or downmilling), ϕhex is the exit angle difference caused by the helix angle and is equal to ϕhex=aptanβ/R.


### 2.1.1. Angular Positionϕjt

Generally, angular positionϕjt is the function of helix angle and spindle speed. Here, it is also relevant to the tooth pitch angle. If a linear tooth pitch variation is shown as(5)ψ=ψ1,ψ2,…,ψN−1,ψN=ψ0,ψ0+Δψ,…,ψ0+N−1Δψ,where ψ0=360/N−N−1Δψ/2, and a sinusoidal spindle speed modulation modeled as(6)Ωt=Ω0+Ω1sin2πTmt=Ω01+RVA⋅sinRVF⋅2π60Ω0t,where RVA=Ω1/Ω0 is the ratio of the speed variation amplitude to the nominal spindle speed and RVF=60/Ω0Tm is the ratio of the speed variation frequency to the nominal spindle speed, are considered the same as [32], the angular position ϕjt can be defined by(7)ϕjt=2π60∫0tΩsds−zγ,ifj=1,2π60∫0tΩsds+∑i=1j−12πψi360−zγ,if1<j≤N.As a spindle rotates in the machining process, the pitch angle of VPC changes periodically at the rotation periodτ0=60/Ω0, whereas the spindle speed varies with the modulation period Tm=60/Ω0/RVF. Therefore, when VSS and VPC are combined, there must be some problems associated with the phase. As shown in Figure 2, one VPC is rotating with Ωt, when its 1-th tooth passes through the positive direction of Y-axis, that is, the red lines; the value of its corresponding instantaneous spindle speed is not equal to Ω0, but the value is between intervals [Ω0−Ω1,Ω0+Ω1], which is depending on the phase of Ωt at this moment. In order to consider this issue, equation (6) is rewritten in a new form including the initial phase as(8)Ωt,η=Ω0+Ω1sin2πTmt+η,where η is defined as the process phase difference between VPC and VSS. Obviously, there are kinds of possible η in VPCVSS milling process, such as π/4, π, and 5π/4. In other words, the value of η is a random one between 0 and 2π in the practical application. As a result, the times to pass through any pitch angle ψj are different when the teeth is rotating at Ωt,η, for example, Cases (1) and (2) shown in Figure 2. Considering the mechanism of regenerative chatter, it means that η will affect the time-varying characteristics of time delays in the milling system. Further, η does make sense for milling stability.Figure 2
Schematic diagram of process phase difference for VPCVSS system.


### 2.1.2. Instantaneous Chip Thicknesshjt

For a milling process, the instantaneous chip thicknesshjt generally consists of two parts, that is, the static chip thickness hjstt mainly determined by the feed per tooth and angular position, and the dynamic displacements vjt and vj−1t, respectively, corresponding to the present tooth j and previous tooth j-1 in the radial direction [1]. Thus, hjt can be expressed as(9)hjt=hjstt,η+vj−1t−vjt=ft,jt,ηsinϕjt−yt−τjt,ηsinϕjt+ytsinϕjt,whereτjt,η represents the time delay between the teeth j and j−1, which is mainly affected by the parameters of VPCVSS, such as spindle speed and pitch angles, and ft,jt,η represents the feed per tooth, which can be described by(10)ft,jt,η=vfτjt,η=fmϑjt,η,where ϑjt,η=Nψj1−RVAsinωmt+η/360, which is applied to represent the relationship between nominal feed per tooth fm and ft,jt,η.


### 2.1.3. Time Delaysτjt,η

Generally, time delay is determined by calculating the rotation time of the tool over the pitch between two successive teeth except for some special cases, such as big cutter runout [15, 37]. For VPCVSS milling, the time delays can be calculated, based on the idea in [5]. Thus, an implicit form including τjt,η is conducted as(11)∫t−τjt,ηtΩ0+Ω1sinωms+η60ds=ψj360.Substituting equation (6) into equation (11) and rearranging both sides of equation (11), then, τjt,η can be approximately expressed by the following relationship:(12)τjt,η≈ψj/κ6,where κ=Ω0+Ω1sinωmt+η. If Ω1 is small enough compared with Ω0, that is, RVA is small, equation (12) can be expanded and then approximated in the form of(13)τjt,η≈τ0,j−τ1,jsinωmt+η,where τ0,j=ψj/Ω0/6 and τ1,j=τ0,jRVA. Here, it should be noted that the approximation errors of equation (13) increase with the increase of RVA. For small RVA, for example, RVA = 0.2, there are also errors although the effects of errors on the stability lob may be ignored.Substituting equations (7) and (9) into equation (3) and considering equation (1), then linearizing the consequent equation on the basis of the ideal in [5], a linear time-periodic DDE can be given by(14)mxx¨t+cxx˙t+kxxt=−ε∑j=1NKjtxt−xt−τjt,η,with(15)ε=γq,Kjt=ft,jq−1t,ηGϕj∫ϕsjtϕejtktsinqϕjtcosϕjt+krsinq+1ϕjtdϕ,


## 2.1.1. Angular Positionϕjt

Generally, angular positionϕjt is the function of helix angle and spindle speed. Here, it is also relevant to the tooth pitch angle. If a linear tooth pitch variation is shown as(5)ψ=ψ1,ψ2,…,ψN−1,ψN=ψ0,ψ0+Δψ,…,ψ0+N−1Δψ,where ψ0=360/N−N−1Δψ/2, and a sinusoidal spindle speed modulation modeled as(6)Ωt=Ω0+Ω1sin2πTmt=Ω01+RVA⋅sinRVF⋅2π60Ω0t,where RVA=Ω1/Ω0 is the ratio of the speed variation amplitude to the nominal spindle speed and RVF=60/Ω0Tm is the ratio of the speed variation frequency to the nominal spindle speed, are considered the same as [32], the angular position ϕjt can be defined by(7)ϕjt=2π60∫0tΩsds−zγ,ifj=1,2π60∫0tΩsds+∑i=1j−12πψi360−zγ,if1<j≤N.As a spindle rotates in the machining process, the pitch angle of VPC changes periodically at the rotation periodτ0=60/Ω0, whereas the spindle speed varies with the modulation period Tm=60/Ω0/RVF. Therefore, when VSS and VPC are combined, there must be some problems associated with the phase. As shown in Figure 2, one VPC is rotating with Ωt, when its 1-th tooth passes through the positive direction of Y-axis, that is, the red lines; the value of its corresponding instantaneous spindle speed is not equal to Ω0, but the value is between intervals [Ω0−Ω1,Ω0+Ω1], which is depending on the phase of Ωt at this moment. In order to consider this issue, equation (6) is rewritten in a new form including the initial phase as(8)Ωt,η=Ω0+Ω1sin2πTmt+η,where η is defined as the process phase difference between VPC and VSS. Obviously, there are kinds of possible η in VPCVSS milling process, such as π/4, π, and 5π/4. In other words, the value of η is a random one between 0 and 2π in the practical application. As a result, the times to pass through any pitch angle ψj are different when the teeth is rotating at Ωt,η, for example, Cases (1) and (2) shown in Figure 2. Considering the mechanism of regenerative chatter, it means that η will affect the time-varying characteristics of time delays in the milling system. Further, η does make sense for milling stability.Figure 2
Schematic diagram of process phase difference for VPCVSS system.


## 2.1.2. Instantaneous Chip Thicknesshjt

For a milling process, the instantaneous chip thicknesshjt generally consists of two parts, that is, the static chip thickness hjstt mainly determined by the feed per tooth and angular position, and the dynamic displacements vjt and vj−1t, respectively, corresponding to the present tooth j and previous tooth j-1 in the radial direction [1]. Thus, hjt can be expressed as(9)hjt=hjstt,η+vj−1t−vjt=ft,jt,ηsinϕjt−yt−τjt,ηsinϕjt+ytsinϕjt,whereτjt,η represents the time delay between the teeth j and j−1, which is mainly affected by the parameters of VPCVSS, such as spindle speed and pitch angles, and ft,jt,η represents the feed per tooth, which can be described by(10)ft,jt,η=vfτjt,η=fmϑjt,η,where ϑjt,η=Nψj1−RVAsinωmt+η/360, which is applied to represent the relationship between nominal feed per tooth fm and ft,jt,η.


## 2.1.3. Time Delaysτjt,η

Generally, time delay is determined by calculating the rotation time of the tool over the pitch between two successive teeth except for some special cases, such as big cutter runout [15, 37]. For VPCVSS milling, the time delays can be calculated, based on the idea in [5]. Thus, an implicit form including τjt,η is conducted as(11)∫t−τjt,ηtΩ0+Ω1sinωms+η60ds=ψj360.Substituting equation (6) into equation (11) and rearranging both sides of equation (11), then, τjt,η can be approximately expressed by the following relationship:(12)τjt,η≈ψj/κ6,where κ=Ω0+Ω1sinωmt+η. If Ω1 is small enough compared with Ω0, that is, RVA is small, equation (12) can be expanded and then approximated in the form of(13)τjt,η≈τ0,j−τ1,jsinωmt+η,where τ0,j=ψj/Ω0/6 and τ1,j=τ0,jRVA. Here, it should be noted that the approximation errors of equation (13) increase with the increase of RVA. For small RVA, for example, RVA = 0.2, there are also errors although the effects of errors on the stability lob may be ignored.Substituting equations (7) and (9) into equation (3) and considering equation (1), then linearizing the consequent equation on the basis of the ideal in [5], a linear time-periodic DDE can be given by(14)mxx¨t+cxx˙t+kxxt=−ε∑j=1NKjtxt−xt−τjt,η,with(15)ε=γq,Kjt=ft,jq−1t,ηGϕj∫ϕsjtϕejtktsinqϕjtcosϕjt+krsinq+1ϕjtdϕ,


## 2.2. Updated SDM

SDM is a common method often applied to the stability analysis of various DDEs where the time delayed terms are approximated by a piecewise constant function while the current time terms are left unchanged. Here, the method is presented for the case of VPCVSS.However, to ensure that equation (14) is periodic and applicable to SDM, the ratio of the spindle modulation period Tm and the rotation period τ0 in the current system must be a rational number, that is, q1Tm=q2τ0 where q1 and q2 are relative prime numbers [5]. Thus, the system period Tc can be expressed as(16)Tc=q1Tmorq2τ0.By introducing the system vectorUt=xt,x˙tT and using Cauchy transformation, the second-order DDE, that is, equation (14) can be written as a system of the first-order DDE:(17)U˙t=ΑtUt+∑j=1NBjtUt−τjt,η,At=At+Tc,Bjt=Bjt+Tc,with(18)At=01−ωnx2−εltmx−2ζxωnx,Bjt=00ε∑j=1Nkjtmx0,lt=∑j=1Nkjt,where ωnx=kx/mx and ζx=cx/2ωnxmx are the angular natural frequency and the damping ratio.The first step of SDM is the construction of the discrete time intervalti,ti+1 by dividing its period Tc into k number with length Δt, that is, Δt=Tc/k. Then, the time-varying delay τjt,η in equation (17) for the interval ti,ti+1 can be written in an average form based on the ideal in [5] as(19)τi,j=1Δt∫titi+1τjt,ηdt=τ0,j−τ1,jci,where(20)ci=kν2π∫i2π/k/ν+ηi+12π/k/ν+ηsintdt,i=0,1,…,k−1ν=TmTc.Then, substituting equation (19) into equation (17), equation (17) can be approximated in [ti, ti+1] as new form yields as(21)U˙t=AiUt+∑j=1NBi,jαi,jUi−mi,j+1+χi,jUi−mi,j=AiUt+Wi,with(22)Ai=1Δt∫titi+1Atdt,Bi,j=1Δt∫titi+1Bjtdt,αi,j=mi,jΔt+Δt/2−τi,jΔt,χi,j=1−αi,j,Ui−mi,j+1=Uti−mi,j+1,Ui−mi,j=Uti−mi,j,mi,j=intτi,j/Δt+1/2,where int (∗) indicates the operation that rounds positive number towards zero. For the initial value U(ti) = Ui, the solution of equation (21) is known: (23)Ut=eAit−tiUi+Ai−1Wi−Ai−1Wi.Substitutingt = ti+1 and Uti+1=Ui+1 into the solution leads to(24)Ui+1=QiUi+∑j=1NMi,j,αUi−mi,j+1+Mi,j,χUi−mi,j,where Qi=eAiΔt,Mi,j,α=eAiΔt−IΑi−1Bi,jαi,j,Mi,j,χ=eAiΔt−IAi−1Bi,jχi,j.Defining the maximum value ofmi,j by n=maxmi,j and state augmentation of equation (24) results in the 2(n+1) dimensional discrete map: (25)Vi+1=ZiVi,with(26)Vi=colUi,Ui−1,…,Ui−n+1,Ui−n,(27)Zi=Qi0⋯000I0⋯0000I⋯000⋮⋮⋱⋮⋮⋮00⋯I0000⋯0I0+∑j=1N0⋯Mi,j,αMi,j,χ⋯00⋯00⋯00⋯00⋯0⋮⋮⋱⋮⋮⋮0⋯00⋯00⋯00⋯0.The horizontal position of the discrete input matricesMi,j,α and Mi,j,χ in the right side of equation (27) begins from the column of 2mi,j−1 and 2mi,j+1, respectively. Here, it should be noted that if one reduction about the matrix Vi is carried out by removing the delayed values of the velocities, like [7], the size of vector Vi can be decreased to n+2, the horizontal position of Mi,j,αand Mi,j,χ will change from the column of mi,j+1 and mi,j+2, respectively.The Floquet transition matrix over the principal periodTc is approximated by coupling the solutions of k successive time intervals Δt:(28)Φ=Zk−1Zk−2,…,Z1Z0.Stability of the investigated system is determined by the eigenvalues of the transition matrixΦ. The system is stable if all eigenvalues of Φ are in modulus less than 1.The main advantage associated with the method in this paper can be summarized as follows:(1)
Compared with the original algorithm [5, 6, 36], the method in this paper can be applied not only to variable pitch milling cutter or variable speed milling process but also to their composite systems. Mathematically, it is suitable for the stability solution of DDEs with multiple variable delays.(2)
Due to the addition of equations (16) and (20), the method in this paper is also applicable to some special cases, where number 1 cannot be divisible by the values of RVF, such as RVF = 0.22, 0.3, or 0.4. For the above cases, this method is applicable, but the method in [32] is not.In addition, it should be noted that the 1-DOF VPCVSS model is proposed in this paper; however, the milling system belongs to the 1-DOF model in general. Thus, to realize a conversion from 1-DOF one to 2-DOF one [38], some modifications need to be carried out as follows: (a) the system equation (equation (1)) should consider the dynamics in both X and Y directions instead of the original X. (b) Correspondingly, the cutting force model in X and Y directions (equations (2) and (3)) should be constructed. (c) The dynamic cutting force thickness (equation (9)) is changed to a coupled model considering X and Y directions. (d) The system equation (equation (17)) and the corresponding state matrix are changed from two dimensions to four ones. (e) Equations (25) and (26) change from 2(n+1) to 4(n+1) or 2n+4 if the delayed values of the velocities are removed, thus, the horizontal position of the discrete input matrices Mi,j,α and Mi,j,χ in equation (27), respectively, begins from the column of 4mj,k,i−3 and 4mj,k,i+1or 2mj,k,i+1 and 2mj,k,i+3.


## 3. Stability Analysis and Discuss

In this section, the verification of the proposed method is carried out firstly. Then, the influence of the cutting conditions, that is, the process phase angle, the helix angle, and the feed per tooth on the stability of the VPCVSS process is investigated in detail. Note that besides the influence parameters shown in respective figure names, other parameters in the following calculation example are as follows: upmilling, the cutting-force exponent isq=0.75, the cutting-force coefficients are kt=107×106N/m1+q and kr=40×106N/m1+q, the mode mass is mx=3.1663kg, the natural frequency is ωn=400Hz, damping ratios are ζx=0.02, the number of the cutter teeth is N=4, the tool radius is R=9.525×10−3m, and pitch angles are ψ=75,85,95,105°). Here, it must be emphasized and noted that the VPC used in this paper belongs to an asymmetric cutter. Correspondingly, its dynamic balancing ability is poor and this will be more significant in high-speed cutting. However, this problem cannot be considered in the following calculation examples, so there may be some prediction errors associated with cutting stability in the high-speed region.


### 3.1. Model Verification

In this section, the validity of the proposed method will be carried out by comparing with previous works and the time-domain simulations.To illustrate the performance of the proposed approach, the related works on VSS in [5] and on VPCVSS in [32] are considered first. Based on the main calculation parameters in [5, 32], the proposed method is utilized to calculate the SLDs for both VSS and VPCVSS cases with two different helix angles under two different radial depths of cut. The related results are shown in Figures 3(a)–3(h), where (a)–(d) and (e)–(h), respectively, correspond to VSS and VPCVSS, whilst the stability boundary labeled by blue and red dotted lines are, respectively, associated with the cases of helix angle β=0° and 30°.Figure 3
Stability diagram of VSS and VPCVSS systems under different cutting parameters. (a)–(h): obtained by the analytic method; (i)-(j): obtained by time-domain simulation.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)If comparing the SLDs marked by blue lines in Figures3(a) and 3(b) with those in Figures 5.17(a) and (b) in [5], and comparing the SLDs marked by the red dotted lines in Figures 3(e)–3(g) with those in Figures 4(a) and 4(b) and 11(b) in [32], one can see that their calculation results are in good agreement obviously. This indicates that the proposed method has good reliability and accuracy in predicting the stability for both VSS and VPCVSS processes when the helix angle is β=0°.Figure 4
Effect of feed per tooth on milling stability under different situations. (a)aD = 50% , q = 0.75; (b) aD = 10%, q = 0.75; (c) aD = 50%, q = 0.94.
(a)(b)(c)Comparing all the blue stability boundaries with the red ones in Figures3(a)–3(h), it can be found that, no matter for VSS (Figures 3(a)–3(d)) or VPCVSS (Figures 3(e)–3(h)) milling processes, the change of helix angle can result in significant stability difference.In order to effectively verify the above predicting trend of stability caused by helix angle, one updated time-domain program (or called time marching solution) is utilized. It should be noted that this program is mainly crafted on the basis of [39], and its application range is expanded from the stability prediction for a single case of VPC to that of VPCVSS. The most critical operation is that the rotation angle is discrete uniformly for one system cycle but time period. Based on the program, detailed information about the amplitude of the vibrations can be obtained and the contour plot of peak-to-peak (PTP) [39] vibration can be generated for a range of spindle speed and axial depth of cut combinations. Therefore, the appearance of nonlinear vibration (chatter) can be well grasped by examining the abrupt change of vibration displacement with the axial depth of cut increasing uniformly. Remarkably, this program can take into account teeth jumping out of contact, workpiece dynamics, regeneration effect, pitch angles, helix angles, and spindle speed variations simultaneously.Based on the above program, time-domain simulations for the plot of PTP are carried out under conditions the same as those of the red stability boundaries in Figures3(g) and 3(h), and the related results are shown in Figures 3(i) and 3(j). Obviously, the prediction results based on the proposed method are in accordance with those by the time-domain simulations. This proves the validity and reliability of the proposed method when the helix angles are considered. On the other hand, there are also some little differences between their results. The reason is that the above time-domain simulations can consider the effect of loss-of-contact character on milling dynamics but the proposed method cannot.In order to detect the stability difference caused by the helix angle more intuitively, two special cutting points are selected, that is, point “A” in Figure3(g) (Ω0=1600rpm, ap=48mm) and point “B” in Figure 3(h) (Ω0=13000rpm, ap=27mm). One can see from the associated SLDs that points “A” and “B” are unstable in the cases of β=0°; however, they will turn to the stable ones when the helix angle changes to β=30°.Their vibration information is simulated in the time domain and the results associated with the cutting vibration time history and periodic points, Poincare section, and nonuniform fast Fourier transform (NuFFT) are shown in Figure5. Note that Figures 5(a) and 5(c), respectively, correspond to points “A” and “B” under β=0°, but Figures 5(b) and 5(d) correspond to β=30°. Consequently, the associated stability states are as follows:(i)
Figures5(a) and 5(c) exhibit two unstable processes due to the fact that the period-doubling (characterized by two points in the Poincare section shown as the middle graphic in Figure 5(a)) and quasiperiodic (characterized by the appearance of a circle map for Poincare section shown as the middle graphic in Figure 5(c)) bifurcation occur, respectively. Meanwhile, the large amplitude vibrations occur at chatter frequency (CF).(ii)
Figures5(b) and 5(d) exhibit two stable processes because they are both characterized by the single fixed point for Poincare mapping (see the middle graphics in Figures 5(b) and 5(d)). There are only modulation frequency (MF) and its harmonics in NuFFTs, and the clearly smaller vibration peaks compared with their counterparts.Figure 5
The information associated with the continuous and 1/Tc-sampled time histories of vibration displacement (x), the Poincare section, and the Fast Fourier Transform for points “A” and “B” in Figures 3(g) and 3(h) under different helix angles β. (a, c) β=0° = 0°; (b, d) β=90°.
(a)(b)(c)(d)Obviously, not to mention CFs, there are so many other frequencies in Figure5, especially for Figures 5(a) and 5(b). The reason is that the modulated frequency is equal to Ω0RVF/60q1 and its harmonics equal to nΩ0RVF/60q1, with n=1,2,… appearing now. As a result, their values are about 5.33n for Figures 5(a) and 5(b) and 108.33 n for Figures 5(c) and 5(d), respectively.All in all, the results based on the theoretical calculation and the simulations are consistent.


### 3.2. Influence of the Process Phase Difference

Because the VPCVSS system is a combination of VPC and VSS, there must be the process phase differences between VPC and VSS as shown in Figure2. In fact, the value of η is almost impossible to be zero in practice; therefore, whether the conclusion based on the assumption of η=0 in [32] is universal or not needs further discussion and research.Starting with above problem, to further evaluate the effect ofη on cutting stability, the SLDs corresponding to four special process phase differences (i.e., η=0, π/2, π, and 3π/2) under the large and small radial depth of cut, are calculated and shown in Figures 6(a) and 6(b). It can be seen that there are only two stability curves with obvious difference for every SLD: one is the overlapping stability curves corresponding to η=0 and π (see the blue line and red dotted line) and the other corresponding to η=π/2 and 3π/2 (see the blue line and green dotted line). From a more detailed observation, it can be found that when η changes from 0 to π/2, the system stability is evolving gradually. At η=π/2, the evolution difference reaches the peak. Then, η is increasing from π/2 to π, its stability continues to evolve gradually; however, the difference decreases on the whole. Finally, when η=π, its stability is the same as the case of η=0. Obviously, there are two kinds of limit shapes of the SLDs, that is, ones at η=0 or π and the others at η=π/2 or 3π/2.Figure 6
Influence of the process phase difference on cutting stability and time delays. (a–d) the SLDs under different process phase difference; (e, f) the variation of the associated time delays.
(a)(b)(c)(d)(e)(f)In order to evaluate the impact ofη in more detail, four different η are chosen in the interval [0, π/2] corresponding to Figure 6(a) and the interval [π/2, π] corresponding to Figure 6(b), respectively, and the related SLDs are calculated and shown in Figures 6(c) and 6(d). It can be seen that, for the above intervals, η always plays an important role and results in a gradual change of SLDs from one limit to another. If this trend is extended to the whole value domain of η, one can easily know that the limit shapes of the SLDs under the influence of η, respectively, occur at ηa=nπ−π and ηb=nπ−π/2 (n = 1, 2, ...), and η obviously leads to the periodic evolution of system stability whose period is π. This is attributed to the following reason. For the VPCVSS processes described in this article, their time delays exist in the form of sine waves as shown in Figures 6(c) and 6(d) (note that they are corresponding to the SLDs in Figures 6(e) and 6(f), resp.). When the values of η are varying, the waves will move along the time axis (shown as the red lines). As a result, the combinations of the multiple time delays in the milling system will be altered, the stability nature will also be changed consequently.In view of the significant influence ofη on cutting stability, it is necessary to consider it for a VPCVSS process. However, as shown in Figure 2, it is obvious that the value of η is random and unmeasured in the actual machining process. Considering its periodic influence on the cutting stability region and the variation limit of the SLDs, one can take the median on the interval [ηa, ηb] as the value of η for obtaining a relatively compromised SLD; thus, the difference between the predicting SLDs and the actual ones can be minimized. With this in mind, η can be defined as ηa+ηa−ηb/2, that is, η=3nπ/4−π/2. Thus, some η, such as π/4 and 5π/4, are suggested and used in the practical application.


### 3.3. Influence of the Feed per Tooth

Feed per tooth is an important basis for the selection of cutting parameters [40, 41]. However, in most previous studies concerning milling chatter, it is ignored usually because it is considered to have no contribution to the dynamic cutting thickness and system stability. Also, one cannot get the effects of feed rate on the stability of milling by a linear cutting force model through the SMD in general. Here, the influence of feed per tooth on stability is investigated deeply, based on the method in this paper. Note that the following will be based on the nominal feed per tooth fm in equation (10) instead of the time-varying one ft,jt,η for convenience.To illustrate the performance of feed per tooth on the dynamic behaviour of the VPCVSS system, SLDs are here provided for four different feed per tooth (i.e.,fm = 0.01, 0.05, 0.1, and 0.2 mm/tooth) under two different radial immersion ratios (10% and 50%), as shown in Figures 4(a) and 4(b). It can be seen that, with the increasing of fm, the corresponding stability boundaries exhibit a trend of upward movement for both big and small radial depth of cut; thus, the limit depth of cut and stability area increases significantly. This is mainly due to the effect of the force-feed nonlinearity, specifically the nonlinear force coefficient q in equation (2), on chatter in the milling process. Because the value of q is equal to 0.75, that is, a number less than 1, a bigger fm will lead to a weaker nonlinear relationship between the cutting force and the dynamic displacements.In order to further investigate the above effects, the SLDs corresponding to the four examples in Figure4(a) are calculated under q = 0.94 and shown in Figure 4(c). It can be seen from Figure 4(c) that the influence of feed per tooth is the same as that of q = 0.75. However, compared with the latter, the former results in a more significant increase in stability. This is mainly due to the fact that since the values of q are significantly less than unity, the nonlinear feed effect is more evident. The conclusions are consistent with those in [41], where the machining stability of the nonlinear feed is investigated for general turning and milling process through experiments and theoretical analysis. Thus, the combination nature in the VPCVSS system does not change the general influence mechanism of the feed per tooth on milling stability.Here, it must be noted that the above analysis conclusions are based on the premise that the value ofq is less than one, which is suitable for some materials, such as aluminum alloy and steel [39]. However, if q = 1, the force-feed relationship in equations (2) and (15) will become linear and the effect of the force-feed nonlinearity on the system dynamic mentioned above will vanish. At this time, only if some other conditions are added, the influence of the feed effect can be considered again. For example, when the multiple time delays caused by the tool runout is considered, a smaller feed will result in a higher stability boundary [15]. For the case of the variable time delay caused by the tool path, [41] exhibits a shift of SLDs with the change of feed. Further, if q > 1, the associated conclusion in this article will be the opposite, that is, the increase of fm results in the smaller limit depth of cut and stability area. Obviously, different considerations and analysis angles will lead to changes in conclusions; this is understandable.Next, the simplification of the time-varying characteristic of the feed per tooth is carried out and its influence on the stability is investigated. For equation (10), if the VSS effect is not considered, that is, RVA = 0, it will degenerate to(29)ft,jt,η=ftmNψj360,if neither the VSS nor VPC effect is considered, it turns to(30)ft,jt,η=ftm.Figure7 shows the SLDs corresponding to ft,jt,η respectively defined by equations (10), (29), and (30). It can be seen that there is a small difference among their stability boundaries for both the radial immersion ratios of 10% and 50%. This indicates that no matter the multi-time-varying characteristics of feed per tooth reflected by equation (10) or the multiconstant ones by equation (29), their influence on the stability is very limited. Thus, for the related general calculations, they both can be simplified to a single constant feed per tooth as shown in equation (30). However, this does not result in too much error for stability prediction.Figure 7
Effect of the different simplified forms of feed per tooth on milling stability.


### 3.4. Influence of the Helix Angle

Although many efforts have been devoted to studying the influence of helix angle on stability in the past, the related works are almost restricted to the regular cutter [33, 34] or variable pitch cutter [19]; there is little work focusing on this topic for VSS milling, let alone VPCVSS. In this section, the above question will be analyzed in detail through combining some new calculations with the ones in Figure 3.As shown in Figures3(a)–3(h), the stability limits increase evidently as the helix angle increases from β=0° to 30°. And this trend is more significant for the small axial depth of cut. This is attributed to the fact that larger helix angle means more obvious system time variation in small immersion ratios, which is characterized as highly interrupted. Further, it can be found that, for VSS milling the effect of helix angle is mainly reflected in the high-speed region (e.g., among 12000–18000 rpm in Figures 3(a)–3(d)), whereas this effect is relatively small for the low-speed region. However, the effect of helix angle is very significant in both high-speed and low-speed regions for VPCVSS milling, for example, 2000–6000 rpm and 12000–18000 rpm in Figures 3(e)–3(h). This may be attributed to the combined effect of VPC and VSS, that is, the positive effect of VPC (especially linear pitches case) on the maximum depth of cut in low speed as shown in [17] and the effective suppression of period double chatter in the high-speed region for VSS milling as shown in [28].In order to evaluate the stability trend caused by helix angle in more detail, the calculations are conducted under four different helix angles (i.e.,β=0°, 15°, 30°, and 45°) and their results are shown in Figure 8. It can be seen that, with the increase of the helix angle, the stable cutting area increases significantly. In Figure 8(a), the so-called “stable islands” [17, 42] gradually appear around 2000 rpm so that the original unstable region turns to stable, whereas the sharp appendages around 16000 rpm gradually turn to “unstable islands” (flip bifurcation region) and disappear finally when β=45°. However, it should be noted that, with the further increase of the helix angle, the influence of the helix angle on the cutting stability becomes weaker and weaker, especially in the high-speed area. This is mainly due to the smaller influence of the larger helix angle on the time-varying characteristics for the VPCVSS system.Figure 8
Effect of helix angle on cutting stability.
(a)(b)


## 3.1. Model Verification

In this section, the validity of the proposed method will be carried out by comparing with previous works and the time-domain simulations.To illustrate the performance of the proposed approach, the related works on VSS in [5] and on VPCVSS in [32] are considered first. Based on the main calculation parameters in [5, 32], the proposed method is utilized to calculate the SLDs for both VSS and VPCVSS cases with two different helix angles under two different radial depths of cut. The related results are shown in Figures 3(a)–3(h), where (a)–(d) and (e)–(h), respectively, correspond to VSS and VPCVSS, whilst the stability boundary labeled by blue and red dotted lines are, respectively, associated with the cases of helix angle β=0° and 30°.Figure 3
Stability diagram of VSS and VPCVSS systems under different cutting parameters. (a)–(h): obtained by the analytic method; (i)-(j): obtained by time-domain simulation.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)If comparing the SLDs marked by blue lines in Figures3(a) and 3(b) with those in Figures 5.17(a) and (b) in [5], and comparing the SLDs marked by the red dotted lines in Figures 3(e)–3(g) with those in Figures 4(a) and 4(b) and 11(b) in [32], one can see that their calculation results are in good agreement obviously. This indicates that the proposed method has good reliability and accuracy in predicting the stability for both VSS and VPCVSS processes when the helix angle is β=0°.Figure 4
Effect of feed per tooth on milling stability under different situations. (a)aD = 50% , q = 0.75; (b) aD = 10%, q = 0.75; (c) aD = 50%, q = 0.94.
(a)(b)(c)Comparing all the blue stability boundaries with the red ones in Figures3(a)–3(h), it can be found that, no matter for VSS (Figures 3(a)–3(d)) or VPCVSS (Figures 3(e)–3(h)) milling processes, the change of helix angle can result in significant stability difference.In order to effectively verify the above predicting trend of stability caused by helix angle, one updated time-domain program (or called time marching solution) is utilized. It should be noted that this program is mainly crafted on the basis of [39], and its application range is expanded from the stability prediction for a single case of VPC to that of VPCVSS. The most critical operation is that the rotation angle is discrete uniformly for one system cycle but time period. Based on the program, detailed information about the amplitude of the vibrations can be obtained and the contour plot of peak-to-peak (PTP) [39] vibration can be generated for a range of spindle speed and axial depth of cut combinations. Therefore, the appearance of nonlinear vibration (chatter) can be well grasped by examining the abrupt change of vibration displacement with the axial depth of cut increasing uniformly. Remarkably, this program can take into account teeth jumping out of contact, workpiece dynamics, regeneration effect, pitch angles, helix angles, and spindle speed variations simultaneously.Based on the above program, time-domain simulations for the plot of PTP are carried out under conditions the same as those of the red stability boundaries in Figures3(g) and 3(h), and the related results are shown in Figures 3(i) and 3(j). Obviously, the prediction results based on the proposed method are in accordance with those by the time-domain simulations. This proves the validity and reliability of the proposed method when the helix angles are considered. On the other hand, there are also some little differences between their results. The reason is that the above time-domain simulations can consider the effect of loss-of-contact character on milling dynamics but the proposed method cannot.In order to detect the stability difference caused by the helix angle more intuitively, two special cutting points are selected, that is, point “A” in Figure3(g) (Ω0=1600rpm, ap=48mm) and point “B” in Figure 3(h) (Ω0=13000rpm, ap=27mm). One can see from the associated SLDs that points “A” and “B” are unstable in the cases of β=0°; however, they will turn to the stable ones when the helix angle changes to β=30°.Their vibration information is simulated in the time domain and the results associated with the cutting vibration time history and periodic points, Poincare section, and nonuniform fast Fourier transform (NuFFT) are shown in Figure5. Note that Figures 5(a) and 5(c), respectively, correspond to points “A” and “B” under β=0°, but Figures 5(b) and 5(d) correspond to β=30°. Consequently, the associated stability states are as follows:(i)
Figures5(a) and 5(c) exhibit two unstable processes due to the fact that the period-doubling (characterized by two points in the Poincare section shown as the middle graphic in Figure 5(a)) and quasiperiodic (characterized by the appearance of a circle map for Poincare section shown as the middle graphic in Figure 5(c)) bifurcation occur, respectively. Meanwhile, the large amplitude vibrations occur at chatter frequency (CF).(ii)
Figures5(b) and 5(d) exhibit two stable processes because they are both characterized by the single fixed point for Poincare mapping (see the middle graphics in Figures 5(b) and 5(d)). There are only modulation frequency (MF) and its harmonics in NuFFTs, and the clearly smaller vibration peaks compared with their counterparts.Figure 5
The information associated with the continuous and 1/Tc-sampled time histories of vibration displacement (x), the Poincare section, and the Fast Fourier Transform for points “A” and “B” in Figures 3(g) and 3(h) under different helix angles β. (a, c) β=0° = 0°; (b, d) β=90°.
(a)(b)(c)(d)Obviously, not to mention CFs, there are so many other frequencies in Figure5, especially for Figures 5(a) and 5(b). The reason is that the modulated frequency is equal to Ω0RVF/60q1 and its harmonics equal to nΩ0RVF/60q1, with n=1,2,… appearing now. As a result, their values are about 5.33n for Figures 5(a) and 5(b) and 108.33 n for Figures 5(c) and 5(d), respectively.All in all, the results based on the theoretical calculation and the simulations are consistent.


## 3.2. Influence of the Process Phase Difference

Because the VPCVSS system is a combination of VPC and VSS, there must be the process phase differences between VPC and VSS as shown in Figure2. In fact, the value of η is almost impossible to be zero in practice; therefore, whether the conclusion based on the assumption of η=0 in [32] is universal or not needs further discussion and research.Starting with above problem, to further evaluate the effect ofη on cutting stability, the SLDs corresponding to four special process phase differences (i.e., η=0, π/2, π, and 3π/2) under the large and small radial depth of cut, are calculated and shown in Figures 6(a) and 6(b). It can be seen that there are only two stability curves with obvious difference for every SLD: one is the overlapping stability curves corresponding to η=0 and π (see the blue line and red dotted line) and the other corresponding to η=π/2 and 3π/2 (see the blue line and green dotted line). From a more detailed observation, it can be found that when η changes from 0 to π/2, the system stability is evolving gradually. At η=π/2, the evolution difference reaches the peak. Then, η is increasing from π/2 to π, its stability continues to evolve gradually; however, the difference decreases on the whole. Finally, when η=π, its stability is the same as the case of η=0. Obviously, there are two kinds of limit shapes of the SLDs, that is, ones at η=0 or π and the others at η=π/2 or 3π/2.Figure 6
Influence of the process phase difference on cutting stability and time delays. (a–d) the SLDs under different process phase difference; (e, f) the variation of the associated time delays.
(a)(b)(c)(d)(e)(f)In order to evaluate the impact ofη in more detail, four different η are chosen in the interval [0, π/2] corresponding to Figure 6(a) and the interval [π/2, π] corresponding to Figure 6(b), respectively, and the related SLDs are calculated and shown in Figures 6(c) and 6(d). It can be seen that, for the above intervals, η always plays an important role and results in a gradual change of SLDs from one limit to another. If this trend is extended to the whole value domain of η, one can easily know that the limit shapes of the SLDs under the influence of η, respectively, occur at ηa=nπ−π and ηb=nπ−π/2 (n = 1, 2, ...), and η obviously leads to the periodic evolution of system stability whose period is π. This is attributed to the following reason. For the VPCVSS processes described in this article, their time delays exist in the form of sine waves as shown in Figures 6(c) and 6(d) (note that they are corresponding to the SLDs in Figures 6(e) and 6(f), resp.). When the values of η are varying, the waves will move along the time axis (shown as the red lines). As a result, the combinations of the multiple time delays in the milling system will be altered, the stability nature will also be changed consequently.In view of the significant influence ofη on cutting stability, it is necessary to consider it for a VPCVSS process. However, as shown in Figure 2, it is obvious that the value of η is random and unmeasured in the actual machining process. Considering its periodic influence on the cutting stability region and the variation limit of the SLDs, one can take the median on the interval [ηa, ηb] as the value of η for obtaining a relatively compromised SLD; thus, the difference between the predicting SLDs and the actual ones can be minimized. With this in mind, η can be defined as ηa+ηa−ηb/2, that is, η=3nπ/4−π/2. Thus, some η, such as π/4 and 5π/4, are suggested and used in the practical application.


## 3.3. Influence of the Feed per Tooth

Feed per tooth is an important basis for the selection of cutting parameters [40, 41]. However, in most previous studies concerning milling chatter, it is ignored usually because it is considered to have no contribution to the dynamic cutting thickness and system stability. Also, one cannot get the effects of feed rate on the stability of milling by a linear cutting force model through the SMD in general. Here, the influence of feed per tooth on stability is investigated deeply, based on the method in this paper. Note that the following will be based on the nominal feed per tooth fm in equation (10) instead of the time-varying one ft,jt,η for convenience.To illustrate the performance of feed per tooth on the dynamic behaviour of the VPCVSS system, SLDs are here provided for four different feed per tooth (i.e.,fm = 0.01, 0.05, 0.1, and 0.2 mm/tooth) under two different radial immersion ratios (10% and 50%), as shown in Figures 4(a) and 4(b). It can be seen that, with the increasing of fm, the corresponding stability boundaries exhibit a trend of upward movement for both big and small radial depth of cut; thus, the limit depth of cut and stability area increases significantly. This is mainly due to the effect of the force-feed nonlinearity, specifically the nonlinear force coefficient q in equation (2), on chatter in the milling process. Because the value of q is equal to 0.75, that is, a number less than 1, a bigger fm will lead to a weaker nonlinear relationship between the cutting force and the dynamic displacements.In order to further investigate the above effects, the SLDs corresponding to the four examples in Figure4(a) are calculated under q = 0.94 and shown in Figure 4(c). It can be seen from Figure 4(c) that the influence of feed per tooth is the same as that of q = 0.75. However, compared with the latter, the former results in a more significant increase in stability. This is mainly due to the fact that since the values of q are significantly less than unity, the nonlinear feed effect is more evident. The conclusions are consistent with those in [41], where the machining stability of the nonlinear feed is investigated for general turning and milling process through experiments and theoretical analysis. Thus, the combination nature in the VPCVSS system does not change the general influence mechanism of the feed per tooth on milling stability.Here, it must be noted that the above analysis conclusions are based on the premise that the value ofq is less than one, which is suitable for some materials, such as aluminum alloy and steel [39]. However, if q = 1, the force-feed relationship in equations (2) and (15) will become linear and the effect of the force-feed nonlinearity on the system dynamic mentioned above will vanish. At this time, only if some other conditions are added, the influence of the feed effect can be considered again. For example, when the multiple time delays caused by the tool runout is considered, a smaller feed will result in a higher stability boundary [15]. For the case of the variable time delay caused by the tool path, [41] exhibits a shift of SLDs with the change of feed. Further, if q > 1, the associated conclusion in this article will be the opposite, that is, the increase of fm results in the smaller limit depth of cut and stability area. Obviously, different considerations and analysis angles will lead to changes in conclusions; this is understandable.Next, the simplification of the time-varying characteristic of the feed per tooth is carried out and its influence on the stability is investigated. For equation (10), if the VSS effect is not considered, that is, RVA = 0, it will degenerate to(29)ft,jt,η=ftmNψj360,if neither the VSS nor VPC effect is considered, it turns to(30)ft,jt,η=ftm.Figure7 shows the SLDs corresponding to ft,jt,η respectively defined by equations (10), (29), and (30). It can be seen that there is a small difference among their stability boundaries for both the radial immersion ratios of 10% and 50%. This indicates that no matter the multi-time-varying characteristics of feed per tooth reflected by equation (10) or the multiconstant ones by equation (29), their influence on the stability is very limited. Thus, for the related general calculations, they both can be simplified to a single constant feed per tooth as shown in equation (30). However, this does not result in too much error for stability prediction.Figure 7
Effect of the different simplified forms of feed per tooth on milling stability.


## 3.4. Influence of the Helix Angle

Although many efforts have been devoted to studying the influence of helix angle on stability in the past, the related works are almost restricted to the regular cutter [33, 34] or variable pitch cutter [19]; there is little work focusing on this topic for VSS milling, let alone VPCVSS. In this section, the above question will be analyzed in detail through combining some new calculations with the ones in Figure 3.As shown in Figures3(a)–3(h), the stability limits increase evidently as the helix angle increases from β=0° to 30°. And this trend is more significant for the small axial depth of cut. This is attributed to the fact that larger helix angle means more obvious system time variation in small immersion ratios, which is characterized as highly interrupted. Further, it can be found that, for VSS milling the effect of helix angle is mainly reflected in the high-speed region (e.g., among 12000–18000 rpm in Figures 3(a)–3(d)), whereas this effect is relatively small for the low-speed region. However, the effect of helix angle is very significant in both high-speed and low-speed regions for VPCVSS milling, for example, 2000–6000 rpm and 12000–18000 rpm in Figures 3(e)–3(h). This may be attributed to the combined effect of VPC and VSS, that is, the positive effect of VPC (especially linear pitches case) on the maximum depth of cut in low speed as shown in [17] and the effective suppression of period double chatter in the high-speed region for VSS milling as shown in [28].In order to evaluate the stability trend caused by helix angle in more detail, the calculations are conducted under four different helix angles (i.e.,β=0°, 15°, 30°, and 45°) and their results are shown in Figure 8. It can be seen that, with the increase of the helix angle, the stable cutting area increases significantly. In Figure 8(a), the so-called “stable islands” [17, 42] gradually appear around 2000 rpm so that the original unstable region turns to stable, whereas the sharp appendages around 16000 rpm gradually turn to “unstable islands” (flip bifurcation region) and disappear finally when β=45°. However, it should be noted that, with the further increase of the helix angle, the influence of the helix angle on the cutting stability becomes weaker and weaker, especially in the high-speed area. This is mainly due to the smaller influence of the larger helix angle on the time-varying characteristics for the VPCVSS system.Figure 8
Effect of helix angle on cutting stability.
(a)(b)


## 4. Conclusion

This paper presents a dynamical model of the milling process of VPC and VSS considering the helix angle and the process phase difference. Then, an updated SDM is proposed to gain the stability chart. After the verification of the proposed method by comparisons with previously published works and the time-domain simulations, lots of influence analyses are conducted deeply. From the stability analysis, the following conclusions can be drawn:(1)
No matter for VSS or VPCVSS milling, the change of helix angle can result in significant stability difference for both small and big radial immersion ratios. However, for VSS milling, the effect of helix angle is mainly reflected in the high-speed region, whereas the effect of helix angle is very significant in both high-speed and low-speed regions for VPCVSS milling process.(2)
Process phase difference has an important influence on the cutting stability, especially for the low-speed domain. Under the influence of the phase difference, there are two kinds of change limits in SLDs, which occur inη=n−1π and nπ/2, n = 1, 2, .... Therefore, it is necessary to consider the phase difference in the stability calculation of VPCVSS or similar system. Considering the randomness and immeasurability of η in the actual process, its recommended value is η=3nπ/4−π/2, such as π/4 and 5π/4.(3)
Due to the effect of the force-feed nonlinearity, as the feed per tooth is increasing, the stability boundaries exhibit a trend of upward movement, both for big and small radial immersion ratios, thus the limit depths of cut and stability area are increasing significantly.(4)
The multi-time-varying characteristic of feed per tooth in VPCVSS milling has almost negligible influence on the stability; thus, it can be simplified to a single constant one for a general calculation.



---

*Source: 6654176-2021-02-09.xml*