# Stability Analysis Method for Periodic Delay Differential Equations with Multiple Distributed and Time-Varying Delays

Authors: Gang Jin; Xinyu Zhang; Kaifei Zhang; Hua Li; Zhanjie Li; Jianxin Han; Houjun Qi
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Dynamic stability problems leading to delay differential equations (DDEs) are found in many different fields of science and engineering. In this paper, a method for stability analysis of periodic DDEs with multiple distributed and time-varying delays is proposed, based on the well-known semidiscretization method. In order to verify the correctness of the proposed method, two typical application examples, i.e., milling process with a variable helix cutter and milling process with variable spindle speed, which can be, respectively, described by DDEs with the multidistributed and time-varying delays are considered. Then, comparisons with prior methods for stability prediction are made to verify the accuracy and efficiency of the proposed approach. As far as the milling process is concerned, the proposed method supplies a generalized algorithm to analyze the stability of the single milling systems associated with variable pith cutter, variable helix cutter, or variable spindle speed; it also can be utilized to analyze the combined systems of the aforementioned cases.


---

## Body

## 1. Introduction

Time-delay systems, where rate of change of state is determined by both present and past state variables, are encountered in many different fields of science and engineering, such as machining processes, chemical processes, wheel dynamics, feedback controller dynamics, and population dynamics. However, time delays are frequently a source of system instability, thus may lead to poor performance, unpleasant noise and sound, or other potential damage for the engineering practice. Therefore, the study of dynamical systems with time delays has received considerable attention in the past decades.Stability charts are a useful tool for studying dynamic problems in engineering because they can present the stability of the linearized system in the plane of the system parameters. As a result of the time delays, the study of stability of the systems becomes an infinite-dimensional problem and the governing equations are delay differential equations (DDEs) rather than the traditional ordinary differential equations. For simple time-delay systems, the stability charts can be derived analytically; however, only numerical techniques can be used for complex systems. Around the dynamic analysis for such systems, a lot of research has been carried out by the scientists and mathematicians committed themselves in the past, and many methods including analytical, semianalytical, and numerical ones have been proposed, such as D-subdivision [1], cluster treatment method [2], Galerkin projection [3], harmonic balance [4], Chebyshev collocation method [5], Lambert W function-based method [6], temporal finite element analysis [7], semidiscretization method (SDM) [8–10], full-discretization method [11], and spectral element method [12, 13].Moreover, many scholars are also committed to dealing with DDEs with various types of time delays caused by more complex engineering practice [14–16]. Taking the typical application of time-delay system, stability of the machining process, as an example, the corresponding model of the turning process is an autonomous DDE, while the milling operation can be described by DDEs with time-periodic coefficients. The associated delays originate from the surface regeneration, which is caused by the current and previous positions of the tool and the workpiece [17]. For some special cases, e.g., the utilization of the unconventional tools or spindle speed, the time delays may be in various forms as well. For the milling cases of variable pitch cutter or cutter runout, they can be characterized by multiple time delays. The frequency-domain model [18, 19], cluster treatment of characteristic roots [20], unified SDM [21], improved FDM [22], and improved SDM [23–25] are utilized to analyze their dynamics. Also, Turner et al. [26], Dombovari et al. [27], Dombovari and Stepan [28], Sims et al. [29], Yusoff and Sims [30], and Jin et al. [31, 32] modeled and analyzed the milling stability for special serrated tools or variable helix tools, where the distributed time delays exist. Besides, DDEs with time-varying time delays corresponding to variable spindle speed milling [33, 34] and that with state-dependent ones associated with the trochoidal path of the teeth [35] are also analyzed by SDM effectively.It should be noted that although the methods in [18–33, 35] have been proved to be effective for the cases with time-varying, multiple or distributed delays, they can only be effective in one or two forms of time delay. For the method in [34], it can be utilized to deal with the DDEs with time-varying and multiple delays; however, it is incapable of solving the distributed delay problem. Obviously, the scope of application of the aforementioned methods is limited. From the point of view of the convenience of calculation and the collaborative optimization of the related parameters, it is necessary to propose a generalized stability solution method.With the aforementioned issue in mind, in this paper, a generalized method is proposed to obtain the stability chart for the periodic DDEs with multiple distributed time-varying delays. The focus of the current work is to present a method which has a greater scope of application. As far as the practical application of milling process is concerned, the method is applicable not only to the stability prediction of variable pitch, variable helix, or variable speed systems but also to the stability analysis of combined systems of variable pitch angle and speed or variable helix angle and speed. The structure of the paper is as follows. In Section2, the mathematical model is introduced. In Section 3, two typical examples are used to verify the effectiveness of the proposed method. In Section 4, conclusions with a brief discussion are presented.


## 2. Mathematical Model

Without loss of generality, n-dimensional linear, time-periodic DDEs with multiple distributed and time-varying delays can be expressed by the following state-space form:(1)y˙t=A0yt+Atyt+∑j=1N∫α1α2Bjt,ϕut−τjt,ϕdϕ,ut=Dyt,where yt∈ℝnis the state, ut∈ℝmis the input, A0 is a constant matrix, At and Bjt,ϕ, respectively, are n×n and n×m periodic coefficient matrices that satisfy At=At+Tand Bjt,ϕ=Bjt+T,ϕ, j = 1, 2, …, N, and τjt,ϕ=τjt+T,ϕ>0, T is the time period, and N is the number of time delay. Considering the idea in [10] and defining Cjt,ϕ=Bjt,ϕD, equation (1) can be written in the following form:(2)y˙t=A0yt+Atyt+∑j=1N∑k=1fC˜j,ktyt−τ˜j,kt,where(3)C˜j,kt=∫α1α1+khCjt,ϕdϕ,k=1,2,...,f−1,∫α1+f−1hα2Cjt,ϕdϕ,k=f.τ˜j,kt=τjt,k−12h,k=1,2,...,f−1,τjt,α2+fh−h2,k=f.f=ceilα2−α1/h, with ceil being the ceiling function (i.e., ceil (x) is the smallest integer not less than x) and hdefines the discretization step for integral interval α1,α2. Then, equation (2) is approximated [10, 11] in interval ti,ti+1 by(4)y˙t=A0yt+A˜ty˜t+∑j=1N∑k=1fC˜j,ktyt−τ˜j,k,i,where(5)A˜t=Ai+Ai+1−AiΔtt−ti,y˜t=yi+yi+1−yiΔtt−ti,C˜j,kt=Cj,k,i+Cj,k,i+1−Cj,k,iΔtt−ti,yt−τ˜j,k,i=yi−mj,k,i+yi+1−mj,k,i−yi−mj,k,iΔtt−ti,τ˜j,k,i=1Δt∫titi+1τ˜j,ktdt,mj,k,i=intτ˜j,k,i+Δt/2Δt.Substituting equation (5) into equation (4) leads to(6)yi+1=Hiyi+Ji∑j=1N∑k=1fRj,k,iyi−mj,k,i+Sj,k,iyi+1−mj,k,i,where(7)Hi=JiΦ0+Fi,Ji=I−Pi−1,Fi=Φ1−2ΔtΦ2+1Δt2Φ3Ai+1ΔtΦ2−1Δt2Φ3Ai+1,Pi=1ΔtΦ2−1Δt2Φ3Ai+1Δt2Φ3Ai+1,Rj,k,i=Φ1−2ΔtΦ2+1Δt2Φ3Cj,k,i+1ΔtΦ2−1Δt2Φ3Cj,k,i+1,Sj,k,i=1ΔtΦ2−1Δt2Φ3Cj,k,i+1Δt2Φ3Cj,k,i+1.Clearly,Φ0, Φ1, Φ2, and Φ3 can be expressed as follows:(8)Φ0=eA0Δt,Φ1=∫0ΔteA0Δt−sds=A0−1Φ0−I,Φ2=∫0ΔtseA0Δt−sds=A0−1Φ1−ΔtI,Φ3=∫0Δts2eA0Δt−sds=A0−12Φ2−Δt2I,where I denotes the identity matrix. Let M=maxmj.k,i and(9)Zi=colyi,yi−1,…,yi−M.Then, combining equation (6) with (9), one can recast it into a discrete map as(10)Zi+1=DiZi,where each Di matrix is given by(11)Di=Hi0⋯000I0⋯0000I⋯000⋮⋮⋱⋮⋮⋮00⋯I0000⋯0I0+∑j=1N∑k=1f0⋯JiRj,k,iJiSj,k,i⋯00⋯00⋯00⋯00⋯0⋮⋱⋮⋮⋱⋮0⋯00⋯00⋯00⋯0.It should be noted that the horizontal position of the discrete input matricesJiRj,k,i and JiSj,k,i in equation (11) depends on the value of mj,k,icorresponding to τ˜j,k,i in equation (5) and, respectively, begins from the column of 2mj,k,i − 1 and 2mj,k,i + 1 for a single degree of freedom (DOF) system (4mj,k,i−3 and 4mj,k,i+1 for a two-DOF system).Then, the approximate Floquet transition matrix is obtained asΦ=Dk−1Dk−2⋯D1D0. Now according to Floquet theory, the stability of the system can be determined: if the moduli of all the eigenvalues of the transition matrix Φ are less than unity, the system is stable; otherwise, it is unstable.


## 3. Verification of the Method

In this section, two typical applications of periodic DDEs associated with milling operations are chosen and analyzed. One is the dynamics of milling process with a variable helix cutter, characterized by multiple distributed time delays [29–32], whereas the other is corresponding to variable spindle speed where multiple time-varying delays exist [33, 34].


### 3.1. Milling with Variable Helix Cutter

#### 3.1.1. Dynamic Model

The two-DOF mathematical model for the milling process with a variable helix cutter [32], as shown in Figure 1, can be written as(12)mx00myx¨ty¨t+cx00cyx˙ty˙t+kx00kyxtyt=∑j=1N∫0ap−hxx,jt,z−hxy,jt,z−hyx,jt,z−hyy,jt,zxt−xt−τjzyt−yt−τjzdz,with the modal mass m•, the spring stiffness k•=2m•ζ•ω•, and the damping coefficient c•=m•ω•2, where •=x and y, ζ•are the damping ratios, and ω•are the natural frequencies. hj,xxt,z, hj,xyt,z, hj,yxt,z, and hj,yyt,zare the cutting force coefficients for the jth tooth, which can be defined as(13)hxx,jt,z=gϕjt,zKtcosϕjt,z+Krsinϕjt,zsinϕjt,z,hxy,jt,z=gϕjt,zKtcosϕjt,z+Krsinϕjt,zcosϕjt,z,hyx,jt,z=gϕjt,z−Ktsinϕjt,z+Krcosϕjt,zsinϕjt,z,hyy,jt,z=gϕjt,z−Ktsinϕjt,z+Krcosϕjt,zcosϕjt,z,where g is a unit step function that determines whether the tooth is in or out of cut, Kt and Kr are the tangential and the radial cutting force coefficients, respectively, and ϕjt,z is the angular position of the jth tooth defined by(14)ϕjt,z=Ωt+∑i=1j−12πψi360−zγ,and the delay τjz can be derived from pitch angle ψj and helix angles βj as(15)τjz=30ψjπΩ−30zγjπRΩ,where γj=tanβj−1−tanβj1<j≤NtanβN−tanβ1j=1, Ω is the spindle speed in rpm, and R is the cutter radius.Figure 1
Schematic mechanical model of a two-DOF milling system with variable pitch/helix tool.DenotingM=mx00my, C=cx00cy, K=kx00ky, qt=xtyt, pt=Mq˙+cq/2, and xt=qtptT, the two-DOF milling model in equation (12) can be represented as(16)y˙t=A0yt+Atyt+∑j=1N∫0apCjt,zyt−τjt,zdz,where(17)A0=−M−1C2M−1CM−1C4−K−CM−12,At=00000000−∑j=1N∫0aphxx,jt,zdz−∑j=1N∫0aphxy,jt,zdz00−∑j=1N∫0aphyx,jt,zdz−∑j=1N∫0aphyy,jt,zdz00,Cjt,z=00000000hxx,jt,zhxy,jt,z00hyx,jt,zhyy,jt,z00.Here, it should be noted that the value ofh in equation (3) can be determined by helix angle β, cutter radius R, and discrete number of time period k and is equal to 2πR/tanβ/k.


#### 3.1.2. Verification of the Proposed Method

To illustrate the performance of the proposed approach for the milling process with a variable helix cutter, the related studies [29, 32] are considered firstly. Based on the proposed approach and the methods in [29, 32], the stability charts are, respectively, calculated in cases of four types of milling cutters with different combinations of pitch angles ψ and helix angles β (see Table 1) and are plotted in Figure 2. The main parameters are as follows: down milling, the cutting-force coefficients Kt=679×106N/m2and Kr=256×106N/m2, the mode masses mx = 1.4986 kg and my = 1.199 kg, the natural frequencies ωx = 563.6 Hz and ωy = 516.21 Hz, the damping ratios ζx = 0.0558 and ζy = 0.025, the number of the cutter teeth N = 4, and the tool radius R=9.525×10−3m.Table 1
Tool parameters.
NumberTool typesPitch anglesψHelix anglesβ1Variable pitch(70, 110, 70, 110)°(30, 30, 30, 30)°2Variable helix 1(90, 90, 90, 90)°(30, 35, 30, 35)°3Variable helix 2(70, 110, 70, 110)°(30, 35, 30, 35)°4Variable helix 3(70, 110, 70, 110)°(35, 30, 35, 30)°Figure 2
The stability charts for 50% radial immersion down milling with different combinations of variable helix angle. (a)ψ = (70, 110, 70, 110)° and β = (30, 30, 30, 30)°; (b) ψ = (90, 90, 90, 90)° and β = (30, 35, 30, 35)°; (c) ψ = (70, 110, 70, 110)° and β = (30, 35, 30, 35)°; (d) ψ = (70, 110, 70, 110)° and β = (35, 30, 35, 30)°.
(a)(b)(c)(d)It can be seen from Figure2 that the results obtained via the proposed method in this paper are in basic agreement with those by the methods in [29, 32]. This indicates that the proposed method is effective. On the other hand, there are also some little differences among the three results. This is mainly attributed to the discrepancies in the method themselves or the selection of discrete parameters.


#### 3.1.3. Comparisons of Methods

In Figure2, three typical methods are utilized to analyze the stability of milling with variable pith/helix cutters. In the following, the analysis of the aforementioned methods will be carried out through comparisons.(1)
The proposed method and the method in [32]:(i)
The proposed method is a numerical one as opposed to the frequency-domain one in [32] that is on the basis of the well-known zeroth-order approximation method [14]. Thus, compared with the method in [32], the proposed method has worse computational efficiency to predict stability charts naturally.(ii)
However, for a traditional milling, the method in [32] cannot accurately predict the stability lobes under small axial depth of cut because of the replacement of time-varying system matrices with time-averaged ones [14]. Moreover, if the nature corresponding to variable pith/helix cutter is also considered, the prediction error may be bigger regardless of low and high radial immersion due to the more markedly time-varying character [24, 25]. Fortunately, the aforementioned problems do not exist in the proposed method. In other words, the proposed method has a much better ability to accurately predict milling stability.(2)
The proposed method and the method in [29]:(i)
Both methods are motivated by the work in [10].(ii)
The method in [29] is on the basis of a new semidiscretization formulation that performs spatial and temporal discretization of the tool, but the proposed one is based on the work of Ding et al. [11], where one utilizes a different concept in the discretization scheme compared with SDM.(iii)
The rates of convergence of both methods are different. In [31], authors presented an updated SDM and pointed out that their method converges much faster than that in [29]. Actually, the method in [31] is a zeroth-SDM, which has similar rate of convergence to the basic algorithm of the proposed method, as known from [36]. By a comprehensive comparison, one can say that the proposed method has better rates of convergence than the method in [29].In addition, the biggest advantage of this method is that besides the aforementioned case related to distributed time delay, it is also effective to deal with the problem with regard to varying time delay, such as the variable spindle speed milling, that will be analyzed in the next section.


### 3.2. Milling with Varying Spindle Speed

#### 3.2.1. Dynamic Model

Based on the nonlinear force model in [10], the one-DOF mathematical model for milling processes with spindle speed variation considering helix angles, as shown in Figure 3, can be written as(18)mxx¨t+cxx˙t+kxxt=−∑j=1N∫0aphjz,txt−xt−τjtdz,where hjt,z are the cutting force coefficients for the jth tooth defined as(19)hjt,z=apqfzq−1gjt,zsinqϕjt,zKtcosϕjt,z+Krsinϕjt,z,where q is the nonlinear parameter in cutting force, and the angular position of tooth j is(20)ϕjt,z=2π60∫0tΩsds+j2πN−zγ,where Ωsis the spindle speed and is assumed to change in the form of a sinusoidal wave [10], which is periodic at a time period T=60/Ω0/RVF, with a nominal value, Ω0, and an amplitude, Ω1=RVA×Ω0, as shown in Figure 3. For this sinusoidal modulation, the shape function is modeled as(21)Ωt=Ω0+Ω1sin2πTt=Ω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/Ω0T is the ratio of the speed variation frequency to the nominal spindle speed. Under these circumstances, the time delay τjt can be obtained in the following relationship:(22)τjt=ψjκ/6,where κ=Ω0+Ω1sinωmt.Figure 3
Schematic mechanical model of the one-DOF milling system with variable spindle speed.Letyt=mxx˙t+mxζxωxxt and xt=xtytT, then equation (18) can be represented as(23)y˙t=A0yt+Atyt+∑j=1N∫0apCjt,zyt−τjtdz,where(24)A0=−ζxωx1mxmxζωx2−mxωx2−ζxωx,At=00−1mx∑j=1N∫0aphjt,zdz0,Cjt,z=001mxhjt,z0.


#### 3.2.2. Results and Discussion

To illustrate the effectiveness of the proposed method for milling process with variable spindle speed, the method and results in [10] are taken into consideration. The stability charts corresponding to the milling processes for RVA = 0.1 and for four different RVF values are calculated using the proposed method and the method in [10] and plotted in Figure 4. The main calculation parameters are as follows: up milling, the cutting-force exponent q = 0.75, the cutting-force coefficients kt=107×106N/m1+q and kr=40×106N/m1+q, the mode mass mx = 3.1663 kg, the natural frequency ωx = 400 Hz, damping ratios ζx = 0.02, the number of the cutter teeth N = 4, the tool radius R=9.525×10−3 m, and the feed per tooth fz = 0.1 mm/tooth. It can be seen from Figure 4 that the results obtained via the above two methods are in close agreement with each other. This indicates the effectiveness of the proposed method for the case of variable spindle speed.Figure 4
Comparison of stability charts for milling processes with sinusoidal spindle speed modulation with RVA = 0.1 in the high-speed domain: (a) RVF = 0.5; (b) RVF = 0.2; (c) RVF = 0.1; (d) RVF = 0.05.
(a)(b)(c)(d)Moreover, because the proposed method is based on the ones in [23, 31, 34], but suitable for variable spindle speed with helix angle, the stability charts associated with helix angle β = 30° are also plotted in Figure 4 for comparison. It can be seen that the effect of helix angle is mainly reflected in the high-speed region (e.g., among 12000–18000 rpm in every graphic of Figure 4), whereas this effect is relatively small for the low-speed region. The results are consistent with the conclusion in [33], where the effective suppression of period double chatter is investigated in the high-speed region for variable spindle speed milling.In addition, two points about the above two methods have to be emphasized here. First, they utilize different policies in the process of equation approximations. For the method in [10], the delayed term is approximated by a linear function of time and the periodic coefficient is approximated by a piecewise constant function. However, the delayed term yt−τ˜j,kt, the state term yt, and the periodic terms At and C˜j,kt in equation (2) are all discretized by linear interpolation for the proposed method. Second, the matrices Φ0, Φ1, Φ2, and Φ3 in equation (8) are dependent on spindle speed but on depth of cut; thus, the calculation time in the process of sweeping the range of the depth of cut will be decreased consequently [11].


## 3.1. Milling with Variable Helix Cutter

### 3.1.1. Dynamic Model

The two-DOF mathematical model for the milling process with a variable helix cutter [32], as shown in Figure 1, can be written as(12)mx00myx¨ty¨t+cx00cyx˙ty˙t+kx00kyxtyt=∑j=1N∫0ap−hxx,jt,z−hxy,jt,z−hyx,jt,z−hyy,jt,zxt−xt−τjzyt−yt−τjzdz,with the modal mass m•, the spring stiffness k•=2m•ζ•ω•, and the damping coefficient c•=m•ω•2, where •=x and y, ζ•are the damping ratios, and ω•are the natural frequencies. hj,xxt,z, hj,xyt,z, hj,yxt,z, and hj,yyt,zare the cutting force coefficients for the jth tooth, which can be defined as(13)hxx,jt,z=gϕjt,zKtcosϕjt,z+Krsinϕjt,zsinϕjt,z,hxy,jt,z=gϕjt,zKtcosϕjt,z+Krsinϕjt,zcosϕjt,z,hyx,jt,z=gϕjt,z−Ktsinϕjt,z+Krcosϕjt,zsinϕjt,z,hyy,jt,z=gϕjt,z−Ktsinϕjt,z+Krcosϕjt,zcosϕjt,z,where g is a unit step function that determines whether the tooth is in or out of cut, Kt and Kr are the tangential and the radial cutting force coefficients, respectively, and ϕjt,z is the angular position of the jth tooth defined by(14)ϕjt,z=Ωt+∑i=1j−12πψi360−zγ,and the delay τjz can be derived from pitch angle ψj and helix angles βj as(15)τjz=30ψjπΩ−30zγjπRΩ,where γj=tanβj−1−tanβj1<j≤NtanβN−tanβ1j=1, Ω is the spindle speed in rpm, and R is the cutter radius.Figure 1
Schematic mechanical model of a two-DOF milling system with variable pitch/helix tool.DenotingM=mx00my, C=cx00cy, K=kx00ky, qt=xtyt, pt=Mq˙+cq/2, and xt=qtptT, the two-DOF milling model in equation (12) can be represented as(16)y˙t=A0yt+Atyt+∑j=1N∫0apCjt,zyt−τjt,zdz,where(17)A0=−M−1C2M−1CM−1C4−K−CM−12,At=00000000−∑j=1N∫0aphxx,jt,zdz−∑j=1N∫0aphxy,jt,zdz00−∑j=1N∫0aphyx,jt,zdz−∑j=1N∫0aphyy,jt,zdz00,Cjt,z=00000000hxx,jt,zhxy,jt,z00hyx,jt,zhyy,jt,z00.Here, it should be noted that the value ofh in equation (3) can be determined by helix angle β, cutter radius R, and discrete number of time period k and is equal to 2πR/tanβ/k.


### 3.1.2. Verification of the Proposed Method

To illustrate the performance of the proposed approach for the milling process with a variable helix cutter, the related studies [29, 32] are considered firstly. Based on the proposed approach and the methods in [29, 32], the stability charts are, respectively, calculated in cases of four types of milling cutters with different combinations of pitch angles ψ and helix angles β (see Table 1) and are plotted in Figure 2. The main parameters are as follows: down milling, the cutting-force coefficients Kt=679×106N/m2and Kr=256×106N/m2, the mode masses mx = 1.4986 kg and my = 1.199 kg, the natural frequencies ωx = 563.6 Hz and ωy = 516.21 Hz, the damping ratios ζx = 0.0558 and ζy = 0.025, the number of the cutter teeth N = 4, and the tool radius R=9.525×10−3m.Table 1
Tool parameters.
NumberTool typesPitch anglesψHelix anglesβ1Variable pitch(70, 110, 70, 110)°(30, 30, 30, 30)°2Variable helix 1(90, 90, 90, 90)°(30, 35, 30, 35)°3Variable helix 2(70, 110, 70, 110)°(30, 35, 30, 35)°4Variable helix 3(70, 110, 70, 110)°(35, 30, 35, 30)°Figure 2
The stability charts for 50% radial immersion down milling with different combinations of variable helix angle. (a)ψ = (70, 110, 70, 110)° and β = (30, 30, 30, 30)°; (b) ψ = (90, 90, 90, 90)° and β = (30, 35, 30, 35)°; (c) ψ = (70, 110, 70, 110)° and β = (30, 35, 30, 35)°; (d) ψ = (70, 110, 70, 110)° and β = (35, 30, 35, 30)°.
(a)(b)(c)(d)It can be seen from Figure2 that the results obtained via the proposed method in this paper are in basic agreement with those by the methods in [29, 32]. This indicates that the proposed method is effective. On the other hand, there are also some little differences among the three results. This is mainly attributed to the discrepancies in the method themselves or the selection of discrete parameters.


### 3.1.3. Comparisons of Methods

In Figure2, three typical methods are utilized to analyze the stability of milling with variable pith/helix cutters. In the following, the analysis of the aforementioned methods will be carried out through comparisons.(1)
The proposed method and the method in [32]:(i)
The proposed method is a numerical one as opposed to the frequency-domain one in [32] that is on the basis of the well-known zeroth-order approximation method [14]. Thus, compared with the method in [32], the proposed method has worse computational efficiency to predict stability charts naturally.(ii)
However, for a traditional milling, the method in [32] cannot accurately predict the stability lobes under small axial depth of cut because of the replacement of time-varying system matrices with time-averaged ones [14]. Moreover, if the nature corresponding to variable pith/helix cutter is also considered, the prediction error may be bigger regardless of low and high radial immersion due to the more markedly time-varying character [24, 25]. Fortunately, the aforementioned problems do not exist in the proposed method. In other words, the proposed method has a much better ability to accurately predict milling stability.(2)
The proposed method and the method in [29]:(i)
Both methods are motivated by the work in [10].(ii)
The method in [29] is on the basis of a new semidiscretization formulation that performs spatial and temporal discretization of the tool, but the proposed one is based on the work of Ding et al. [11], where one utilizes a different concept in the discretization scheme compared with SDM.(iii)
The rates of convergence of both methods are different. In [31], authors presented an updated SDM and pointed out that their method converges much faster than that in [29]. Actually, the method in [31] is a zeroth-SDM, which has similar rate of convergence to the basic algorithm of the proposed method, as known from [36]. By a comprehensive comparison, one can say that the proposed method has better rates of convergence than the method in [29].In addition, the biggest advantage of this method is that besides the aforementioned case related to distributed time delay, it is also effective to deal with the problem with regard to varying time delay, such as the variable spindle speed milling, that will be analyzed in the next section.


## 3.1.1. Dynamic Model

The two-DOF mathematical model for the milling process with a variable helix cutter [32], as shown in Figure 1, can be written as(12)mx00myx¨ty¨t+cx00cyx˙ty˙t+kx00kyxtyt=∑j=1N∫0ap−hxx,jt,z−hxy,jt,z−hyx,jt,z−hyy,jt,zxt−xt−τjzyt−yt−τjzdz,with the modal mass m•, the spring stiffness k•=2m•ζ•ω•, and the damping coefficient c•=m•ω•2, where •=x and y, ζ•are the damping ratios, and ω•are the natural frequencies. hj,xxt,z, hj,xyt,z, hj,yxt,z, and hj,yyt,zare the cutting force coefficients for the jth tooth, which can be defined as(13)hxx,jt,z=gϕjt,zKtcosϕjt,z+Krsinϕjt,zsinϕjt,z,hxy,jt,z=gϕjt,zKtcosϕjt,z+Krsinϕjt,zcosϕjt,z,hyx,jt,z=gϕjt,z−Ktsinϕjt,z+Krcosϕjt,zsinϕjt,z,hyy,jt,z=gϕjt,z−Ktsinϕjt,z+Krcosϕjt,zcosϕjt,z,where g is a unit step function that determines whether the tooth is in or out of cut, Kt and Kr are the tangential and the radial cutting force coefficients, respectively, and ϕjt,z is the angular position of the jth tooth defined by(14)ϕjt,z=Ωt+∑i=1j−12πψi360−zγ,and the delay τjz can be derived from pitch angle ψj and helix angles βj as(15)τjz=30ψjπΩ−30zγjπRΩ,where γj=tanβj−1−tanβj1<j≤NtanβN−tanβ1j=1, Ω is the spindle speed in rpm, and R is the cutter radius.Figure 1
Schematic mechanical model of a two-DOF milling system with variable pitch/helix tool.DenotingM=mx00my, C=cx00cy, K=kx00ky, qt=xtyt, pt=Mq˙+cq/2, and xt=qtptT, the two-DOF milling model in equation (12) can be represented as(16)y˙t=A0yt+Atyt+∑j=1N∫0apCjt,zyt−τjt,zdz,where(17)A0=−M−1C2M−1CM−1C4−K−CM−12,At=00000000−∑j=1N∫0aphxx,jt,zdz−∑j=1N∫0aphxy,jt,zdz00−∑j=1N∫0aphyx,jt,zdz−∑j=1N∫0aphyy,jt,zdz00,Cjt,z=00000000hxx,jt,zhxy,jt,z00hyx,jt,zhyy,jt,z00.Here, it should be noted that the value ofh in equation (3) can be determined by helix angle β, cutter radius R, and discrete number of time period k and is equal to 2πR/tanβ/k.


## 3.1.2. Verification of the Proposed Method

To illustrate the performance of the proposed approach for the milling process with a variable helix cutter, the related studies [29, 32] are considered firstly. Based on the proposed approach and the methods in [29, 32], the stability charts are, respectively, calculated in cases of four types of milling cutters with different combinations of pitch angles ψ and helix angles β (see Table 1) and are plotted in Figure 2. The main parameters are as follows: down milling, the cutting-force coefficients Kt=679×106N/m2and Kr=256×106N/m2, the mode masses mx = 1.4986 kg and my = 1.199 kg, the natural frequencies ωx = 563.6 Hz and ωy = 516.21 Hz, the damping ratios ζx = 0.0558 and ζy = 0.025, the number of the cutter teeth N = 4, and the tool radius R=9.525×10−3m.Table 1
Tool parameters.
NumberTool typesPitch anglesψHelix anglesβ1Variable pitch(70, 110, 70, 110)°(30, 30, 30, 30)°2Variable helix 1(90, 90, 90, 90)°(30, 35, 30, 35)°3Variable helix 2(70, 110, 70, 110)°(30, 35, 30, 35)°4Variable helix 3(70, 110, 70, 110)°(35, 30, 35, 30)°Figure 2
The stability charts for 50% radial immersion down milling with different combinations of variable helix angle. (a)ψ = (70, 110, 70, 110)° and β = (30, 30, 30, 30)°; (b) ψ = (90, 90, 90, 90)° and β = (30, 35, 30, 35)°; (c) ψ = (70, 110, 70, 110)° and β = (30, 35, 30, 35)°; (d) ψ = (70, 110, 70, 110)° and β = (35, 30, 35, 30)°.
(a)(b)(c)(d)It can be seen from Figure2 that the results obtained via the proposed method in this paper are in basic agreement with those by the methods in [29, 32]. This indicates that the proposed method is effective. On the other hand, there are also some little differences among the three results. This is mainly attributed to the discrepancies in the method themselves or the selection of discrete parameters.


## 3.1.3. Comparisons of Methods

In Figure2, three typical methods are utilized to analyze the stability of milling with variable pith/helix cutters. In the following, the analysis of the aforementioned methods will be carried out through comparisons.(1)
The proposed method and the method in [32]:(i)
The proposed method is a numerical one as opposed to the frequency-domain one in [32] that is on the basis of the well-known zeroth-order approximation method [14]. Thus, compared with the method in [32], the proposed method has worse computational efficiency to predict stability charts naturally.(ii)
However, for a traditional milling, the method in [32] cannot accurately predict the stability lobes under small axial depth of cut because of the replacement of time-varying system matrices with time-averaged ones [14]. Moreover, if the nature corresponding to variable pith/helix cutter is also considered, the prediction error may be bigger regardless of low and high radial immersion due to the more markedly time-varying character [24, 25]. Fortunately, the aforementioned problems do not exist in the proposed method. In other words, the proposed method has a much better ability to accurately predict milling stability.(2)
The proposed method and the method in [29]:(i)
Both methods are motivated by the work in [10].(ii)
The method in [29] is on the basis of a new semidiscretization formulation that performs spatial and temporal discretization of the tool, but the proposed one is based on the work of Ding et al. [11], where one utilizes a different concept in the discretization scheme compared with SDM.(iii)
The rates of convergence of both methods are different. In [31], authors presented an updated SDM and pointed out that their method converges much faster than that in [29]. Actually, the method in [31] is a zeroth-SDM, which has similar rate of convergence to the basic algorithm of the proposed method, as known from [36]. By a comprehensive comparison, one can say that the proposed method has better rates of convergence than the method in [29].In addition, the biggest advantage of this method is that besides the aforementioned case related to distributed time delay, it is also effective to deal with the problem with regard to varying time delay, such as the variable spindle speed milling, that will be analyzed in the next section.


## 3.2. Milling with Varying Spindle Speed

### 3.2.1. Dynamic Model

Based on the nonlinear force model in [10], the one-DOF mathematical model for milling processes with spindle speed variation considering helix angles, as shown in Figure 3, can be written as(18)mxx¨t+cxx˙t+kxxt=−∑j=1N∫0aphjz,txt−xt−τjtdz,where hjt,z are the cutting force coefficients for the jth tooth defined as(19)hjt,z=apqfzq−1gjt,zsinqϕjt,zKtcosϕjt,z+Krsinϕjt,z,where q is the nonlinear parameter in cutting force, and the angular position of tooth j is(20)ϕjt,z=2π60∫0tΩsds+j2πN−zγ,where Ωsis the spindle speed and is assumed to change in the form of a sinusoidal wave [10], which is periodic at a time period T=60/Ω0/RVF, with a nominal value, Ω0, and an amplitude, Ω1=RVA×Ω0, as shown in Figure 3. For this sinusoidal modulation, the shape function is modeled as(21)Ωt=Ω0+Ω1sin2πTt=Ω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/Ω0T is the ratio of the speed variation frequency to the nominal spindle speed. Under these circumstances, the time delay τjt can be obtained in the following relationship:(22)τjt=ψjκ/6,where κ=Ω0+Ω1sinωmt.Figure 3
Schematic mechanical model of the one-DOF milling system with variable spindle speed.Letyt=mxx˙t+mxζxωxxt and xt=xtytT, then equation (18) can be represented as(23)y˙t=A0yt+Atyt+∑j=1N∫0apCjt,zyt−τjtdz,where(24)A0=−ζxωx1mxmxζωx2−mxωx2−ζxωx,At=00−1mx∑j=1N∫0aphjt,zdz0,Cjt,z=001mxhjt,z0.


### 3.2.2. Results and Discussion

To illustrate the effectiveness of the proposed method for milling process with variable spindle speed, the method and results in [10] are taken into consideration. The stability charts corresponding to the milling processes for RVA = 0.1 and for four different RVF values are calculated using the proposed method and the method in [10] and plotted in Figure 4. The main calculation parameters are as follows: up milling, the cutting-force exponent q = 0.75, the cutting-force coefficients kt=107×106N/m1+q and kr=40×106N/m1+q, the mode mass mx = 3.1663 kg, the natural frequency ωx = 400 Hz, damping ratios ζx = 0.02, the number of the cutter teeth N = 4, the tool radius R=9.525×10−3 m, and the feed per tooth fz = 0.1 mm/tooth. It can be seen from Figure 4 that the results obtained via the above two methods are in close agreement with each other. This indicates the effectiveness of the proposed method for the case of variable spindle speed.Figure 4
Comparison of stability charts for milling processes with sinusoidal spindle speed modulation with RVA = 0.1 in the high-speed domain: (a) RVF = 0.5; (b) RVF = 0.2; (c) RVF = 0.1; (d) RVF = 0.05.
(a)(b)(c)(d)Moreover, because the proposed method is based on the ones in [23, 31, 34], but suitable for variable spindle speed with helix angle, the stability charts associated with helix angle β = 30° are also plotted in Figure 4 for comparison. It can be seen that the effect of helix angle is mainly reflected in the high-speed region (e.g., among 12000–18000 rpm in every graphic of Figure 4), whereas this effect is relatively small for the low-speed region. The results are consistent with the conclusion in [33], where the effective suppression of period double chatter is investigated in the high-speed region for variable spindle speed milling.In addition, two points about the above two methods have to be emphasized here. First, they utilize different policies in the process of equation approximations. For the method in [10], the delayed term is approximated by a linear function of time and the periodic coefficient is approximated by a piecewise constant function. However, the delayed term yt−τ˜j,kt, the state term yt, and the periodic terms At and C˜j,kt in equation (2) are all discretized by linear interpolation for the proposed method. Second, the matrices Φ0, Φ1, Φ2, and Φ3 in equation (8) are dependent on spindle speed but on depth of cut; thus, the calculation time in the process of sweeping the range of the depth of cut will be decreased consequently [11].


## 3.2.1. Dynamic Model

Based on the nonlinear force model in [10], the one-DOF mathematical model for milling processes with spindle speed variation considering helix angles, as shown in Figure 3, can be written as(18)mxx¨t+cxx˙t+kxxt=−∑j=1N∫0aphjz,txt−xt−τjtdz,where hjt,z are the cutting force coefficients for the jth tooth defined as(19)hjt,z=apqfzq−1gjt,zsinqϕjt,zKtcosϕjt,z+Krsinϕjt,z,where q is the nonlinear parameter in cutting force, and the angular position of tooth j is(20)ϕjt,z=2π60∫0tΩsds+j2πN−zγ,where Ωsis the spindle speed and is assumed to change in the form of a sinusoidal wave [10], which is periodic at a time period T=60/Ω0/RVF, with a nominal value, Ω0, and an amplitude, Ω1=RVA×Ω0, as shown in Figure 3. For this sinusoidal modulation, the shape function is modeled as(21)Ωt=Ω0+Ω1sin2πTt=Ω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/Ω0T is the ratio of the speed variation frequency to the nominal spindle speed. Under these circumstances, the time delay τjt can be obtained in the following relationship:(22)τjt=ψjκ/6,where κ=Ω0+Ω1sinωmt.Figure 3
Schematic mechanical model of the one-DOF milling system with variable spindle speed.Letyt=mxx˙t+mxζxωxxt and xt=xtytT, then equation (18) can be represented as(23)y˙t=A0yt+Atyt+∑j=1N∫0apCjt,zyt−τjtdz,where(24)A0=−ζxωx1mxmxζωx2−mxωx2−ζxωx,At=00−1mx∑j=1N∫0aphjt,zdz0,Cjt,z=001mxhjt,z0.


## 3.2.2. Results and Discussion

To illustrate the effectiveness of the proposed method for milling process with variable spindle speed, the method and results in [10] are taken into consideration. The stability charts corresponding to the milling processes for RVA = 0.1 and for four different RVF values are calculated using the proposed method and the method in [10] and plotted in Figure 4. The main calculation parameters are as follows: up milling, the cutting-force exponent q = 0.75, the cutting-force coefficients kt=107×106N/m1+q and kr=40×106N/m1+q, the mode mass mx = 3.1663 kg, the natural frequency ωx = 400 Hz, damping ratios ζx = 0.02, the number of the cutter teeth N = 4, the tool radius R=9.525×10−3 m, and the feed per tooth fz = 0.1 mm/tooth. It can be seen from Figure 4 that the results obtained via the above two methods are in close agreement with each other. This indicates the effectiveness of the proposed method for the case of variable spindle speed.Figure 4
Comparison of stability charts for milling processes with sinusoidal spindle speed modulation with RVA = 0.1 in the high-speed domain: (a) RVF = 0.5; (b) RVF = 0.2; (c) RVF = 0.1; (d) RVF = 0.05.
(a)(b)(c)(d)Moreover, because the proposed method is based on the ones in [23, 31, 34], but suitable for variable spindle speed with helix angle, the stability charts associated with helix angle β = 30° are also plotted in Figure 4 for comparison. It can be seen that the effect of helix angle is mainly reflected in the high-speed region (e.g., among 12000–18000 rpm in every graphic of Figure 4), whereas this effect is relatively small for the low-speed region. The results are consistent with the conclusion in [33], where the effective suppression of period double chatter is investigated in the high-speed region for variable spindle speed milling.In addition, two points about the above two methods have to be emphasized here. First, they utilize different policies in the process of equation approximations. For the method in [10], the delayed term is approximated by a linear function of time and the periodic coefficient is approximated by a piecewise constant function. However, the delayed term yt−τ˜j,kt, the state term yt, and the periodic terms At and C˜j,kt in equation (2) are all discretized by linear interpolation for the proposed method. Second, the matrices Φ0, Φ1, Φ2, and Φ3 in equation (8) are dependent on spindle speed but on depth of cut; thus, the calculation time in the process of sweeping the range of the depth of cut will be decreased consequently [11].


## 4. Conclusion

In this work, an improved semidiscretization method for stability analysis of periodic DDEs with multiple distributed time-periodic delays is proposed. Then, two typical application examples, i.e., variable spindle speed milling system associated with time-varying delays and variable helix cutter milling system associated with multiple distributed delays, are utilized to demonstrate the effectiveness of the proposed algorithm. Through comparison with prior works, it is found that the results obtained by the presented method are in close agreement with those by the prior methods.Furthermore, for the milling process, the proposed method actually provides a generalized algorithm, which can be utilized to predict stability not only for the single milling process with variable pitch cutter, variable helix cutter, or variable spindle speed but also for combined processes, such as variable pitch angle and variable speed, variable helix angle and variable speed. Meanwhile, the application scope of the existing methods, e.g., the ones in [24, 31, 34], has been expanded to a certain extent.



---

*Source: 1982363-2020-07-25.xml*