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

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

#### 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:where is the state, is the input, is a constant matrix, and , respectively, are and periodic coefficient matrices that satisfy and , *j* = 1, 2, …, *N*, and , is the time period, and is the number of time delay. Considering the idea in [10] and defining , equation (1) can be written in the following form:where

, with ceil being the ceiling function (i.e., ceil (*x*) is the smallest integer not less than *x*) and defines the discretization step for integral interval . Then, equation (2) is approximated [10, 11] in interval bywhere

Substituting equation (5) into equation (4) leads towhere

Clearly, , , , and can be expressed as follows:where **I** denotes the identity matrix. Let and

Then, combining equation (6) with (9), one can recast it into a discrete map aswhere each matrix is given by

It should be noted that the horizontal position of the discrete input matrices and in equation (11) depends on the value of corresponding to in equation (5) and, respectively, begins from the column of 2 − 1 and 2 + 1 for a single degree of freedom (DOF) system ( and for a two-DOF system).

Then, the approximate Floquet transition matrix is obtained as . 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 aswith the modal mass , the spring stiffness , and the damping coefficient , where and , are the damping ratios, and are the natural frequencies. , , , and are the cutting force coefficients for the *j*th tooth, which can be defined aswhere is a unit step function that determines whether the tooth is in or out of cut, and are the tangential and the radial cutting force coefficients, respectively, and is the angular position of the *j*th tooth defined byand the delay can be derived from pitch angle and helix angles aswhere , is the spindle speed in rpm, and is the cutter radius.

Denoting , , , , , and , the two-DOF milling model in equation (12) can be represented aswhere

Here, it should be noted that the value of in equation (3) can be determined by helix angle , cutter radius , and discrete number of time period and is equal to .

###### 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 and , the mode masses = 1.4986 kg and = 1.199 kg, the natural frequencies = 563.6 Hz and = 516.21 Hz, the damping ratios = 0.0558 and = 0.025, the number of the cutter teeth = 4, and the tool radius *m*.

**(a)**

**(b)**

**(c)**

**(d)**

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]. 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 Figure 2, 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 aswhere are the cutting force coefficients for the *j*th tooth defined aswhere *q* is the nonlinear parameter in cutting force, and the angular position of tooth *j* iswhere is the spindle speed and is assumed to change in the form of a sinusoidal wave [10], which is periodic at a time period , with a nominal value, , and an amplitude, , as shown in Figure 3. For this sinusoidal modulation, the shape function is modeled aswhere is the ratio of the speed variation amplitude to the nominal spindle speed and is the ratio of the speed variation frequency to the nominal spindle speed. Under these circumstances, the time delay can be obtained in the following relationship:where .

Let and , then equation (18) can be represented aswhere

###### 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 = 0.75, the cutting-force coefficients and , the mode mass = 3.1663 kg, the natural frequency = 400 Hz, damping ratios = 0.02, the number of the cutter teeth = 4, the tool radius m, and the feed per tooth = 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.

**(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 , the state term , and the periodic terms and in equation (2) are all discretized by linear interpolation for the proposed method. Second, the matrices , , , and 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.

#### Data Availability

The 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

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (11702192 and 51405343), Tianjin Science and Technology Planning Project (18JCZDJC10050), School Project of Tianjin Vocational and Technical Normal University (KJ1909), and Innovation Team Training Plan of Tianjin Universities and Colleges (TD13-5096).