# A Method for Stability Analysis of Periodic Delay Differential Equations with Multiple Time-Periodic Delays

Authors: Gang Jin; Houjun Qi; Zhanjie Li; Jianxin Han; Hua Li
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Delay differential equations (DDEs) are widely utilized as the mathematical models in engineering fields. In this paper, a method is proposed to analyze the stability characteristics of periodic DDEs with multiple time-periodic delays. Stability charts are produced for two typical examples of time-periodic DDEs about milling chatter, including the variable-spindle speed milling system with one-time-periodic delay and variable pitch cutter milling system with multiple delays. The simulations show that the results gained by the proposed method are in close agreement with those existing in the past literature. This indicates the effectiveness of our method in terms of time-periodic DDEs with multiple time-periodic delays. Moreover, for milling processes, the proposed method further provides a generalized algorithm, which possesses a good capability to predict the stability lobes for milling operations with variable pitch cutter or variable-spindle speed.


---

## Body

## 1. Introduction

Time-delay systems widely exist in engineering and science, where the rate of change of state is determined by both present and past state variables, such as machining processes [1, 2], wheel dynamics [3, 4], feedback controller [5, 6], gene expression dynamics [7], and population dynamics [8, 9]. However, for some of above applications, the time delay in the dynamic system may lead to instability, poor performance, or other types of potential damage. Therefore, it is necessary for engineers and scientists to research the dynamics of these systems to reduce or avoid such problems.Compared to the finite dimensional dynamics for systems without time delay, time-delay systems have infinite-dimensional dynamics and are usually described by delay differential equations (DDEs). Their stability properties can be analyzed through obtaining the stability charts that show the stable and unstable domains. For example, a stable milling process can be realized by choosing the corresponding parameter from a stability lobe diagram (SLD), which is a function of spindle speed and depth of cut parameters. Thus, more and more attention has been paid on this issue and many analytical and numerical methods have been developed to derive the stability conditions for the system parameters.By using theD-subdivision method, Bhatt and Hsu [10] determined stability criteria for second-order scalar DDEs. Budak and Altıntaş [11, 12] and Merdol and Altintas [13] proposed a method in frequency domain called multifrequency solution. By employing a shifted Chebyshev polynomial approximation, Butcher et al. [14, 15] presented a new technique to study the stability properties of dynamic systems by obtaining an approximate monodromy matrix. Insperger and Stepan [16–18] proposed a known method called semidiscretization method (SDM), which is based on the discretization of the DDEs and approximates their infinite-dimensional phase space by a finite discrete map in time domain. Bayly et al. [19] carried out a temporal finite element analysis for solving the DDEs, which are written in the form of a state space model and discretizing the time interval of interest into a finite number of temporal elements. Based on the direct integration scheme, Ding et al. [20], Liu et al. [21], and Jin et al. [22] used a full-discretization method to gain stability chart efficiently. Recently, Khasawneh and Mann [23] and Lehotzky et al. [24] presented a numerical algorithm called spectral element method. This method has good efficiency because of its highly accurate numerical quadratures for the integral terms.SDM is a known and widely used method to determine stability charts for general time-periodic DDEs arising in different engineering problems. In this paper, based on SDM, a generalized method for periodic DDEs with multiple time-periodic delays is proposed to obtain the stability chart of DDEs. 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

The general form of linear, time-periodic DDEs with multiple time-periodic delays can be expressed as(1)y˙t=A0yt+Atyt-∑j=1NBjtut-τjtut=Dyt,where y(t)∈Rn is the state, u(t)∈Rm is the input, A0 is a constant matrix, A(t) and Bj(t), respectively, are n×n and n×m periodic coefficient matrices that satisfy A(t)=A(t+T) and Bj(t)=Bj(t+T),  j=1,2,…,N, and τj(t)=τj(t+T)>0, D is an m×n constant matrix, T is the time period, and N is the number of time-delays. Note that (1) can also be written in the form (2)y˙t=A0yt+Atyt-∑j=1NCjtyt-τjtwith Cj(t)=Bj(t)D.Consider that the periodT is divided into k number of discrete time intervals, such that each interval length Δt=T/k. Introduce symbol [ti,ti+1] to represent the ith time interval, i∈Z+. Here, ti means the ith time node and is equal to iΔt. Thus, considering the idea in [2], the averaged delay for the discretization interval [ti,ti+1] is defined as follows:(3)τi,j=1Δt∫titi+1τjtdt=τ0,j-τ1,jci,where(4)ci=k2π∫i2π/ki+12π/ksintdt,i=0,1,…,k-1.The number of intervalsmi,j related to the delay item τi,j can be approximately obtained by(5)mi,j=int⁡τi,j+Δt/2Δt,where int(∗) indicates the operation that rounds positive number towards zero.Substituting (3) into (2) and solving it as an ordinary differential equation over the discretization period [ti,ti+1] with initial condition x(ti)=xi, the following equation is derived:(6)yt=eA0t-tiyi+∫titeA0t-ξAξyξ-∑j=1NCjξyξ-τi,jdξ.Substitutingt=ti+1 into (6), then it can be equivalently expressed as(7)yti+1=eA0Δtyi+∫titi+1eA0ti+1-ξAξyξ-∑j=1NCjξyξ-τi,jdξ.In [ti,ti+1],  A(t),  Cj(t),  y(t), and y(t-τi,j) are defined as follows: (8)At=Ai+Ai+1-AiΔtt-ti,Cjt=Ci,j+Ci+1,j-Ci,jΔtt-ti,yt=yi+yi+1-yiΔtt-ti,yt-τi,j=βi,jyi-mi,j+αi,jyi+1-mi,j,where Ai=A(ti),  Ci,j=Cj(ti),  yi=y(ti),  αi,j=(mi,jΔt+Δt/2-τi,j)/Δt, and βi,j=1-αi,j. Here, it should be noted that the approximation of y(t-τi,j) in (8) is the same as that for the so-called zero-order SDM in [2].Substituting (8) into (7) leads to(9)yi+1=Φ0+Fiyi+Piyi+1-∑j=1Nβi,jRi,jyi-mi,j+αi,jRi,jyi+1-mi,j,where(10)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,Ri,j=Φ1-1ΔtΦ2Ci,j+1ΔtΦ2Ci+1,j.Clearly, Φ0, Φ1, Φ2, and Φ3 can be expressed as follows:(11)Φ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=max⁡(mi.j) and(12)Zi=col⁡yi,yi-1⋯yi-M;then combining (9) and (12), one can be recast into a discrete map as(13)Zi+1=DiZi,where each Di matrix is given by(14)Di=Hi+1Φ0+Fi0⋯000I0⋯0000I⋯000⋮⋮⋱⋮⋮⋮00⋯I0000⋯0I0+∑j=1N0⋯-Hi+1αi,jRi,j-Hi+1βi,jRi,j⋯00⋯00⋯00⋯00⋯0⋮⋮⋱⋮⋮⋮0⋯00⋯00⋯00⋯0,where Hi+1=(I-Pi)-1. The horizontal position of the discrete input matrices Hi+1αi,jRi,j and Hi+1βi,jRi,j in (14) depends on the value of mi,j corresponding to τi,j and they, respectively, begin from the column of 2mi,j-1 and 2mi,j+1 for a single DOF system as opposed to the column of 4mi,j-3 and 4mi,j+1 for a two-DOF one.Based on (13) and (14), the following mathematical expressions can be established by coupling the solutions of the k successive time intervals in period T: (15)Vk=ΦV0=Dk-1Dk-2⋯D1D0V0,where Φ is the Floquet transition matrix that gives the connection between Vk and V0. According to the Floquet theory, the stability of the system is determined using the following criterion. If the moduli of all the eigenvalues of the transition matrix Φ are less than unity, the system is stable. Otherwise, it is unstable.Here, it should be noted that the matrixVi can be reduced because only the delayed positions show up in the governing equation of the milling process. Thus, the size of the approximation vector in (14) could be reduced by removing the delayed values of the velocities, such that the size of vector Vi can be decreased to M+2 for a single DOF system and 2M+4 for a two-DOF system. This can give some additional improvement in the computational time for the proposed method.


## 3. Verification of Method

There are several numerical and semianalytical techniques to determine the stability conditions for periodic DDEs. However, most of them were developed with the aim of constructing stability charts for milling processes, such as the analysis of the milling system with runout [25], with variable pitch/helix cutter [26–30], with variable-spindle speed [31–33], or with serrated cutter [34, 35]. In order to verify the proposed method, two typical milling operations are chosen and considered. The first is the varying spindle speed process, which can be described by a DDE with time-periodic delay in general. The other is the milling process with variable pitch cutters, which is often characterized by a DDE with multiple delays. Both methods are known means to influence and to prevent chatter vibration in milling.


### 3.1. Milling with Varying Spindle Speed

Generally, the mathematical models for milling processes with spindle speed variation can be written as(16)y˙t=A0yt+Atyt-Btyt-τt;that is, (2) is degenerated into one with one-time-periodic delay. For a single DOF system in [2], the matrices in (16) have the form (17)A0=01-ωn2-2ζωn,At=Bt=00-Gtm0,where m is the mode mass, ωn is the natural frequency, ζ isthe damping ratios, and G(t) is the specific directional factor and has the form(18)Gt=ap∑j=1Ngjtsin⁡ϕjtKtcos⁡ϕjt+Krsin⁡ϕjt,where ap is the axial depth of cut, N is the number of teeth, Kt and Kr are the linearized cutting coefficients in tangential and radial directions, gj(t) denotes whether the jth tooth is cutting, and the angular position of tooth j is(19)ϕjt=2π60∫0tΩsds+j2πN,where Ω(s) is the spindle speed and is assumed to change in the form of a sinusoidal wave, 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 1. For this sinusoidal modulation, the shape function is modeled as(20)Ωt=Ω0+Ω1sin⁡2πTt=Ω01+RVA·sin⁡RVF·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.Figure 1
Schematic drawing of the sinusoidal modulation of the spindle speed.To illustrate the effectiveness of the approach method for milling with spindle speed, the method and results in [2] are taken into consideration. Here, it should be noted that the delayed term is approximated by a linear function of time and the periodic coefficient is approximated by a piecewise constant function for the method in [2]. However, for the proposed method, the delayed term y(ξ-τi,j), the state term y(ξ), and the periodic terms A(ξ) and Cj(ξ) in (7) are all discretized by linear interpolation (see (8)). Thus, different policies are utilized in the process of equation approximations for two methods. Figure 2 illustrates the stability charts that correspond to the milling processes for RVA = 0.1 and for four different RVF values using the proposed method and the method in [2]. The parameters are as follows. The cutting processes using a 4-flute tool (N=4) with zero helix angles are considered under half-immersion up-milling. The cutting force coefficients are Kt=800×106N/m2 and Kr=300×106N/m2. The mode mass is m=3.1663kg, the natural frequency is ωn = 400 Hz, and damping ratios is ξ=0.02. It can be seen from Figure 2 that the results obtained via the proposed method in this paper are in close agreement with those in [2].Figure 2
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)Meanwhile, the computational times corresponding to every graphic in Figure2 are also recorded to evaluate the efficiency of proposed method. Here, considering the assumption q1T=q2τ0  (τ0 is the tooth passing period) with q1 and q2 being relatively prime [2] and the equation q1/q2=RVF/N obtained consequently, if the resolution of τ0 is adopted as 40 and q1=1, the resolutions of period T are 320, 800, 1600, and 3200 for RVF = 0.5, 0.2, 0.1, and 0.05, respectively. For a 200×100 grid of the spindle speed and the depth of cut and a personal computer (Intel(R) Core(TM) i5-2300, 2.8 GHz, 3 GB), the computational times are, respectively, 436 s, 1026 s, 2012 s, and 4132 s corresponding to RVF = 0.5, 0.2, 0.1, and 0.05 for the proposed method as opposed to 1953 s, 4757 s, 9185 s, and 18341 s for the method in [2] using our own codes. Time costs reduce nearly by 70% for every case. Obviously, the low computational cost of our method is illustrated. The reason about the cost reduction can be explained as follows. The matrices Φ0, Φ1, Φ2, and Φ3 in (11) are dependent on spindle speed but on depth of cut. Consequently, they are not needed to calculate in the process of sweeping the range of the depth of cut for the proposed method. However, this is also necessary for the method in [2]. Thus, for a parameter plane formed by the spindle speed and the cutting depth and divided into a Ns×Nd size grid, the method in [2] must be calculated Ns×Nd×k times to obtain a stability chart, but only Ns times for the method in this paper.


### 3.2. Milling with Variable Pitch Cutter

Considering a system for milling process with variable pitch cutter [29, 30] as shown in Figure 3, the mathematical models can be written as(21)y˙t=A0yt+Atyt-∑j=1NBjtyt-τj;that is, the original system of (2) is degenerated into one with multiple delays in this case. τj is the pitch period corresponding to the pitch angle ψ (see the right graphic of Figure 3), and the matrices A0, A(t), and Bj(t) in (21) have the form(22)A0=00100001-ωnx20-2ζxωnx00-ωny20-2ζyωny,At=∑j=1NBjt,Bjt=00000000-apmxhxx,jt-apmyhxy,jt00-apmxhyx,jt-apmyhyy,jt00,where ζx and ζy are the damping ratios, ωx and ωy are the natural frequencies, and mx and my are the modal masses of the cutter. hxx,j(t), hxy,j(t), hyx,jt,, and hyy,j(t) are the cutting force coefficients for the jth tooth defined as(23)hxx,jt=gϕjtKtcos⁡ϕjt+Krsin⁡ϕjtsin⁡ϕjt,hxy,jt=gϕjtKtcos⁡ϕjt+Krsin⁡ϕjtcos⁡ϕjt,hyx,jt=gϕjt-Ktsin⁡ϕjt+Krcos⁡ϕjtsin⁡ϕjt,hyy,jt=gϕjt-Ktsin⁡ϕjt+Krcos⁡ϕjtcos⁡ϕjt.Figure 3
Schematic mechanical model of a system for milling process with variable pitch cutter.To illustrate the performance of the proposed approach on uniform and variable pitch milling tools, the frequency-domain method published in [29] is considered. The comparison of the results using the proposed approach and the method [29] is carried out for both uniform and variable pitch cutter milling, as shown in Figure 4. The main system parameters are down-milling, half-immersion, the number of the cutter teeth which is N=4, the natural frequencies which are ωx=563.6 Hz and ωy=516.21 Hz, the damping ratios which are ζx=0.0558 and ζy=0.025, the modal masses which are mx=1.4986 kg and my=1.199 kg, and the cutting force coefficients which are Kt=679×106N/m2and Kt=256×106N/m2. It can be seen from Figure 4(a) that the proposed method agrees closely with the results of analytical method for the uniform pitch cutter milling. For the variable pitch cutter as shown in Figure 4(b), two methods gain consistent predicting results, except for the high-speed domain at approximately 8500 rpm, where a clear deviation occurred. Reference [27] chose one point (ap = 5 mm and Ω = 8500 rpm) in this deviation and showed its stabilization by time-domain simulations. The reason of this phenomenon is as follows. The proposed method is based on the time-periodic cutting force coefficients (see (23)), rather than the simplified time-averaged ones in [29]. Thus, the stability prediction by our method is more reasonable and has better accuracy.Figure 4
Comparison of the predicted stability lobes: (a) uniform pitch cutter withψ= [90°, 90°, 90°, 90°] and (b) variable pitch cutter with ψ= [70°, 110°, 70°, 110°].
(a)
(b)


## 3.1. Milling with Varying Spindle Speed

Generally, the mathematical models for milling processes with spindle speed variation can be written as(16)y˙t=A0yt+Atyt-Btyt-τt;that is, (2) is degenerated into one with one-time-periodic delay. For a single DOF system in [2], the matrices in (16) have the form (17)A0=01-ωn2-2ζωn,At=Bt=00-Gtm0,where m is the mode mass, ωn is the natural frequency, ζ isthe damping ratios, and G(t) is the specific directional factor and has the form(18)Gt=ap∑j=1Ngjtsin⁡ϕjtKtcos⁡ϕjt+Krsin⁡ϕjt,where ap is the axial depth of cut, N is the number of teeth, Kt and Kr are the linearized cutting coefficients in tangential and radial directions, gj(t) denotes whether the jth tooth is cutting, and the angular position of tooth j is(19)ϕjt=2π60∫0tΩsds+j2πN,where Ω(s) is the spindle speed and is assumed to change in the form of a sinusoidal wave, 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 1. For this sinusoidal modulation, the shape function is modeled as(20)Ωt=Ω0+Ω1sin⁡2πTt=Ω01+RVA·sin⁡RVF·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.Figure 1
Schematic drawing of the sinusoidal modulation of the spindle speed.To illustrate the effectiveness of the approach method for milling with spindle speed, the method and results in [2] are taken into consideration. Here, it should be noted that the delayed term is approximated by a linear function of time and the periodic coefficient is approximated by a piecewise constant function for the method in [2]. However, for the proposed method, the delayed term y(ξ-τi,j), the state term y(ξ), and the periodic terms A(ξ) and Cj(ξ) in (7) are all discretized by linear interpolation (see (8)). Thus, different policies are utilized in the process of equation approximations for two methods. Figure 2 illustrates the stability charts that correspond to the milling processes for RVA = 0.1 and for four different RVF values using the proposed method and the method in [2]. The parameters are as follows. The cutting processes using a 4-flute tool (N=4) with zero helix angles are considered under half-immersion up-milling. The cutting force coefficients are Kt=800×106N/m2 and Kr=300×106N/m2. The mode mass is m=3.1663kg, the natural frequency is ωn = 400 Hz, and damping ratios is ξ=0.02. It can be seen from Figure 2 that the results obtained via the proposed method in this paper are in close agreement with those in [2].Figure 2
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)Meanwhile, the computational times corresponding to every graphic in Figure2 are also recorded to evaluate the efficiency of proposed method. Here, considering the assumption q1T=q2τ0  (τ0 is the tooth passing period) with q1 and q2 being relatively prime [2] and the equation q1/q2=RVF/N obtained consequently, if the resolution of τ0 is adopted as 40 and q1=1, the resolutions of period T are 320, 800, 1600, and 3200 for RVF = 0.5, 0.2, 0.1, and 0.05, respectively. For a 200×100 grid of the spindle speed and the depth of cut and a personal computer (Intel(R) Core(TM) i5-2300, 2.8 GHz, 3 GB), the computational times are, respectively, 436 s, 1026 s, 2012 s, and 4132 s corresponding to RVF = 0.5, 0.2, 0.1, and 0.05 for the proposed method as opposed to 1953 s, 4757 s, 9185 s, and 18341 s for the method in [2] using our own codes. Time costs reduce nearly by 70% for every case. Obviously, the low computational cost of our method is illustrated. The reason about the cost reduction can be explained as follows. The matrices Φ0, Φ1, Φ2, and Φ3 in (11) are dependent on spindle speed but on depth of cut. Consequently, they are not needed to calculate in the process of sweeping the range of the depth of cut for the proposed method. However, this is also necessary for the method in [2]. Thus, for a parameter plane formed by the spindle speed and the cutting depth and divided into a Ns×Nd size grid, the method in [2] must be calculated Ns×Nd×k times to obtain a stability chart, but only Ns times for the method in this paper.


## 3.2. Milling with Variable Pitch Cutter

Considering a system for milling process with variable pitch cutter [29, 30] as shown in Figure 3, the mathematical models can be written as(21)y˙t=A0yt+Atyt-∑j=1NBjtyt-τj;that is, the original system of (2) is degenerated into one with multiple delays in this case. τj is the pitch period corresponding to the pitch angle ψ (see the right graphic of Figure 3), and the matrices A0, A(t), and Bj(t) in (21) have the form(22)A0=00100001-ωnx20-2ζxωnx00-ωny20-2ζyωny,At=∑j=1NBjt,Bjt=00000000-apmxhxx,jt-apmyhxy,jt00-apmxhyx,jt-apmyhyy,jt00,where ζx and ζy are the damping ratios, ωx and ωy are the natural frequencies, and mx and my are the modal masses of the cutter. hxx,j(t), hxy,j(t), hyx,jt,, and hyy,j(t) are the cutting force coefficients for the jth tooth defined as(23)hxx,jt=gϕjtKtcos⁡ϕjt+Krsin⁡ϕjtsin⁡ϕjt,hxy,jt=gϕjtKtcos⁡ϕjt+Krsin⁡ϕjtcos⁡ϕjt,hyx,jt=gϕjt-Ktsin⁡ϕjt+Krcos⁡ϕjtsin⁡ϕjt,hyy,jt=gϕjt-Ktsin⁡ϕjt+Krcos⁡ϕjtcos⁡ϕjt.Figure 3
Schematic mechanical model of a system for milling process with variable pitch cutter.To illustrate the performance of the proposed approach on uniform and variable pitch milling tools, the frequency-domain method published in [29] is considered. The comparison of the results using the proposed approach and the method [29] is carried out for both uniform and variable pitch cutter milling, as shown in Figure 4. The main system parameters are down-milling, half-immersion, the number of the cutter teeth which is N=4, the natural frequencies which are ωx=563.6 Hz and ωy=516.21 Hz, the damping ratios which are ζx=0.0558 and ζy=0.025, the modal masses which are mx=1.4986 kg and my=1.199 kg, and the cutting force coefficients which are Kt=679×106N/m2and Kt=256×106N/m2. It can be seen from Figure 4(a) that the proposed method agrees closely with the results of analytical method for the uniform pitch cutter milling. For the variable pitch cutter as shown in Figure 4(b), two methods gain consistent predicting results, except for the high-speed domain at approximately 8500 rpm, where a clear deviation occurred. Reference [27] chose one point (ap = 5 mm and Ω = 8500 rpm) in this deviation and showed its stabilization by time-domain simulations. The reason of this phenomenon is as follows. The proposed method is based on the time-periodic cutting force coefficients (see (23)), rather than the simplified time-averaged ones in [29]. Thus, the stability prediction by our method is more reasonable and has better accuracy.Figure 4
Comparison of the predicted stability lobes: (a) uniform pitch cutter withψ= [90°, 90°, 90°, 90°] and (b) variable pitch cutter with ψ= [70°, 110°, 70°, 110°].
(a)
(b)


## 4. Conclusion

In this work, an improved semidiscretization algorithm is proposed to obtain the stability char for DDEs with multiple time-periodic delays. Two milling examples, variable-spindle speed milling system with one-time-periodic delay and variable pitch cutter milling system with multiple delays, are utilized to demonstrate effectiveness of the proposed algorithm. Through the comparison with prior works, it is found that the results gained by the presented method in this paper are in close agreement with those existing in the past literature. Moreover, the proposed method also has good computational efficiency. Here, it should be noted that if discussing the milling process only, the proposed method is a generalized algorithm, which can consider the milling processes with variable pitch cutter and variable-spindle speed simultaneously.



---

*Source: 9490142-2017-09-17.xml*