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

Dynamic stability problems leading to delay differential equations (DDEs) are found inmany 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 semidiscretizationmethod. In order to verify the correctness of the proposedmethod, 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.*en, comparisons with prior methods for stability prediction aremade to verify the accuracy and efficiency of the proposed approach. As far as themilling process is concerned, the proposedmethod 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.


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. erefore, 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][9][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][15][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. e 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.
It should be noted that although the methods in [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][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. e 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. e 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
Without loss of generality, n-dimensional linear, timeperiodic DDEs with multiple distributed and timevarying delays can be expressed by the following statespace form: _ y(t) � A 0 y(t) + A(t)y(t) where y(t) ∈ R n is the state, u(t) ∈ R m is the input, A 0 is a constant matrix, A(t) and B j (t, ϕ), respectively, are n × n and n × m periodic coefficient matrices that satisfy A(t) � A(t + T)and B j (t, ϕ) � B j (t + T, ϕ), j � 1, 2, . . ., N, and τ j (t, ϕ) � τ j (t + T, ϕ) > 0, T is the time period, and N is the number of time delay. Considering the idea in [10] and defining C j (t, ϕ) � B j (t, ϕ)D, equation (1) can be written in the following form: where 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 ]. en, equation (2) is approximated [10,11] Substituting equation (5) into equation (4) leads to 2 Mathematical Problems in Engineering Clearly, Φ 0 , Φ 1 , Φ 2 , and Φ 3 can be expressed as follows: where I denotes the identity matrix. Let M � max m j.k,i and en, combining equation (6) with (9), one can recast it into a discrete map as where each D i matrix is given by It should be noted that the horizontal position of the discrete input matrices J i R j,k,i and J i S j,k,i in equation (11) depends on the value of m j,k,i corresponding to τ j,k,i in equation (5) and, respectively, begins from the column of 2m j,k,i − 1 and 2m j,k,i + 1 for a single degree of freedom (DOF) system (4m j,k,i − 3 and 4m j,k,i + 1 for a two-DOF system).
en, the approximate Floquet transition matrix is 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.

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][30][31][32], whereas the other is corresponding to variable spindle speed where multiple time-varying delays exist [33,34].

Milling with Variable Helix Cutter
e two-DOF mathematical model for the milling process with a variable helix cutter [32], as shown in Figure 1, can be written as with the modal mass m (•) , the spring stiffness , where (•) � x and y, ζ (•) are the damping ratios, and ω (•) are the natural frequencies. h j,xx (t, z), h j,xy (t, z), h j,yx (t, z), and h j,yy (t, z)are the cutting force coefficients for the jth tooth, which can be defined as where g is a unit step function that determines whether the tooth is in or out of cut, K t and K r are the tangential and the radial cutting force coefficients, respectively, and ϕ j (t, z) is the angular position of the jth tooth defined by and the delay τ j (z) can be derived from pitch angle ψ j and helix angles β j as Mathematical Problems in Engineering where , Ω is the spindle speed in rpm, and R is the cutter radius.
, p(t) � M _ q + cq/2, and x(t) � q(t) p(t) T , the two-DOF milling model in equation (12) can be represented as where Here, it should be noted that the value of h 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.

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. e main parameters are as follows: down milling, the cutting-force coefficients K t � 679 × 10 6 N/m 2 and K r � 256 × 10 6 N/m 2 , the mode masses m x � 1.4986 kg and m y � 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 − 3 m.
It can be seen from Figure 2 that the results obtained via the proposed method in this paper are in basic agreement with those by the methods in [29,32]. is indicates that the proposed method is effective. On the other hand, there are also some little differences among the three results. is is mainly attributed to the discrepancies in the method themselves or the selection of discrete parameters. Figure 2, three typical methods are utilized to analyze the stability of milling with variable pith/helix cutters. In the following, the analysis of    Figure 1: Schematic mechanical model of a two-DOF milling system with variable pitch/helix tool.  The method in this paper The method in Ref. [29] The method in Ref. [32]  The method in this paper The method in Ref. [29] The method in Ref. [32]  The method in this paper The method in Ref. [29] The method in Ref. [32]  The method in this paper The method in Ref. [29] The method in Ref. [32] (d) the aforementioned methods will be carried out through comparisons.

Comparisons of Methods. In
(1) e proposed method and the method in [32]: (i) e proposed method is a numerical one as opposed to the frequency-domain one in [32] that is on the basis of the well-known zerothorder approximation method [14]. us, 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 timevarying 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) e proposed method and the method in [29]: (i) Both methods are motivated by the work in [10].
(ii) e 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) e 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.

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 where h j (t, z) are the cutting force coefficients for the jth tooth defined as where q is the nonlinear parameter in cutting force, and the angular position of tooth j is where Ω(s)is 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 where RVA � Ω 1 /Ω 0 is the ratio of the speed variation amplitude to the nominal spindle speed and RVF � 60/(Ω 0 T) is the ratio of the speed variation frequency to the nominal spindle speed. Under these circumstances, the time delay τ j (t) can be obtained in the following relationship: where κ � Ω 0 + Ω 1 sin(ω m t).
Let y(t) � m x _ x(t) + m x ζ x ω x x(t) and x(t) � x(t) y(t)] T , then equation (18) can be represented as where Workpiece Milling tool Feed

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. e 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. e main calculation parameters are as follows: up milling, the cutting-force exponent q � 0.75, the cutting-force coefficients k t � 107× 10 6 N/m 1+q and k r � 40 × 10 6 N/m 1+q , the mode mass m x � 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 f z � 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. is indicates the effectiveness of the proposed method for the case of variable spindle speed.
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. e 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 y(t − τ j,k (t)), the state term y(t), and the periodic terms A(t) and C j,k (t) 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].

Conclusion
In this work, an improved semidiscretization method for stability analysis of periodic DDEs with multiple distributed time-periodic delays is proposed. en, 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. rough 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.

Data Availability
e MATLAB data used to support the findings of this study have been deposited in the "Baidu online disk" repository (https://pan.baidu.com/s/10mA5uypwbzUd2QGIO8drlA, extraction code: aakg).

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.