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

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 proposedmethod 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.


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 infinitedimensional 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 the -subdivision method, Bhatt and Hsu [10] determined stability criteria for second-order scalar DDEs.Budak and Altıntas ¸ [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][17][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 timeperiodic delays is proposed to obtain the stability chart of DDEs.The structure of the paper is as follows.In Section 2, 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.

Mathematical Model
The general form of linear, time-periodic DDEs with multiple time-periodic delays can be expressed as where y() ∈ R  is the state, u() ∈ R  is the input, A 0 is a constant matrix, A() and B  (), respectively, are  ×  and × periodic coefficient matrices that satisfy A() = A(+) and B  () = B  ( + ),  = 1, 2, . . ., , and   () =   ( + ) > 0, D is an  ×  constant matrix,  is the time period, and  is the number of time-delays.Note that (1) can also be written in the form with C  () = B  ().
Consider that the period  is divided into  number of discrete time intervals, such that each interval length Δ = /.Introduce symbol [  ,  +1 ] to represent the th time interval,  ∈ Z + .Here,   means the th time node and is equal to Δ.Thus, considering the idea in [2], the averaged delay for the discretization interval [  ,  +1 ] is defined as follows: where The number of intervals  , related to the delay item  , can be approximately obtained by 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 [  ,  +1 ] with initial condition (  ) =   , the following equation is derived: Substituting  =  +1 into (6), then it can be equivalently expressed as In [  ,  +1 ], A(), C  (), y(), and y( −  , ) are defined as follows: where A  = A(  ), C , = C  (  ), y  = y(  ),  , = ( , Δ + Δ/2 −  , )/Δ, and  , = 1 −  , .Here, it should be noted that the approximation of y( −  , ) in ( 8) is the same as that for the so-called zero-order SDM in [2].Substituting ( 8) into ( 7) leads to where Clearly, Φ 0 , Φ 1 , Φ 2 , and Φ 3 can be expressed as follows: where I denotes the identity matrix.Let  = max( . ) and then combining ( 9) and ( 12), one can be recast into a discrete map as where each D  matrix is given by where H +1 = (I−P  ) −1 .The horizontal position of the discrete input matrices H +1  , R , and H +1  , R , in ( 14) depends on the value of  , corresponding to  , and they, respectively, begin from the column of 2 , − 1 and 2 , + 1 for a single DOF system as opposed to the column of 4 , −3 and 4 , +1 for a two-DOF one.
Based on ( 13) and ( 14), the following mathematical expressions can be established by coupling the solutions of the  successive time intervals in period : where Φ is the Floquet transition matrix that gives the connection between V  and V 0 .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 matrix V  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 V  can be decreased to  + 2 for a single DOF system and 2 + 4 for a two-DOF system.This can give some additional improvement in the computational time for the proposed method.

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][27][28][29][30], with variable-spindle speed [31][32][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.

Milling with
Varying Spindle Speed.Generally, the mathematical models for milling processes with spindle speed variation can be written as 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 where  is the mode mass,   is the natural frequency,  is the damping ratios, and () is the specific directional factor and has the form where   is the axial depth of cut,  is the number of teeth,   and   are the linearized cutting coefficients in tangential and radial directions,   () denotes whether the th tooth is cutting, and the angular position of tooth  is T/4 T/2 3T/4 T 0 t (s) where Ω() is the spindle speed and is assumed to change in the form of a sinusoidal wave, which is periodic at a time period  = 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 where RVA = Ω 1 /Ω 0 is the ratio of the speed variation amplitude to the nominal spindle speed and RVF = 60/(Ω 0 ) is the ratio of the speed variation frequency to the nominal 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( −  , ), the state term y(), and the periodic terms A() and C  () 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 ( = 4) with zero helix angles are considered under half-immersion up-milling.The cutting force coefficients are   = 800 × 10 6 N/m 2 and   = 300 × 10 6 N/m 2 .The mode mass is  = 3.1663 kg, the natural frequency is   = 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].
Meanwhile, the computational times corresponding to every graphic in Figure 2 are also recorded to evaluate the efficiency of proposed method.Here, considering the assumption  1  =  2  0 ( 0 is the tooth passing period) with  1 and  2 being relatively prime [2] and the equation  1 / 2 = RVF/ obtained consequently, if the resolution of  0 is adopted as 40 and  1 = 1, the resolutions of period  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   ×   size grid, the method in [2] must be calculated   ×   ×  times to obtain a stability chart, but only   times for the method in this paper.

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 that is, the original system of ( 2) is degenerated into one with multiple delays in this case.  is the pitch period corresponding to the pitch angle  (see the right graphic of Figure 3), and the matrices A 0 , A(), and B  () in ( 21) have the form where   and   are the damping ratios,   and   are the natural frequencies, and   and   are the modal masses of the cutter.ℎ , (), ℎ , (), ℎ , (),, and ℎ , () are the cutting force coefficients for the th tooth defined as The method in [2] The proposed method ×10 The method in [2] The proposed method The method in [2] The proposed method The method in [2] The proposed method To illustrate the performance of the proposed approach on uniform and variable pitch milling tools, the frequencydomain 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  The proposed method The method in [29] 3000 4000 5000 6000 7000 8000 9000 10000 2000 Spindle speed (rev/min)   = 256 × 10 6 N/m 2 .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 (  = 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 timeperiodic 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.

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, variablespindle 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.

Figure 1 :
Figure 1: Schematic drawing of the sinusoidal modulation of the spindle speed.

Figure 4 .
The main system parameters are down-milling, half-immersion, the number of the cutter teeth which is  = 4, the natural frequencies which are   = 563.6Hz and   = 516.21Hz, the damping ratios which are   = 0.0558 and   = 0.025, the modal masses which are   = 1.4986 kg and   = 1.199 kg, and the cutting force coefficients which are   = 679 × 10 6 N/m 2 and Mathematical Problems in Engineering Tool circumference (variable pitch tool) Tool circumference (regular tool)

Figure 3 :
Figure 3: Schematic mechanical model of a system for milling process with variable pitch cutter.