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.

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 the -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 [1618] 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 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

The general form of linear, time-periodic DDEs with multiple time-periodic delays can be expressed aswhere is the state, is the input, is a constant matrix, and , respectively, are and periodic coefficient matrices that satisfy and ,  , and , 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 .

Consider that the period is divided into number of discrete time intervals, such that each interval length . Introduce symbol to represent the th time interval, . Here, means the th time node and is equal to . Thus, considering the idea in [2], the averaged delay for the discretization interval is defined as follows:where

The number of intervals related to the delay item can be approximately obtained bywhere 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 with initial condition , the following equation is derived:

Substituting into (6), then it can be equivalently expressed asIn ,  ,  ,  , and are defined as follows: where ,  ,  ,  , and . Here, it should be noted that the approximation of in (8) is the same as that for the so-called zero-order SDM in [2].

Substituting (8) into (7) leads towhereClearly, , , , and can be expressed as follows:where denotes the identity matrix. Let andthen combining (9) and (12), one can be recast into a discrete map aswhere each matrix is given bywhere . The horizontal position of the discrete input matrices and in (14) depends on the value of corresponding to and they, respectively, begin from the column of and for a single DOF system as opposed to the column of and 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 and . 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 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 can be decreased to for a single DOF system and 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 [2630], with variable-spindle speed [3133], 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 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 formwhere 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 iswhere is the spindle speed and is assumed to change in the form of a sinusoidal wave, which is periodic at a time period , with a nominal value, , and an amplitude , as shown in Figure 1. 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.

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 , the state term , and the periodic terms and 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 with zero helix angles are considered under half-immersion up-milling. The cutting force coefficients are and . The mode mass is kg, the natural frequency is = 400 Hz, and damping ratios is . 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   ( is the tooth passing period) with and being relatively prime [2] and the equation obtained consequently, if the resolution of is adopted as 40 and , the resolutions of period are 320, 800, 1600, and 3200 for RVF = 0.5, 0.2, 0.1, and 0.05, respectively. For a 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 , , , and 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.

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 asthat 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 , , and in (21) have the formwhere 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

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 , the natural frequencies which are  Hz and  Hz, the damping ratios which are and , the modal masses which are  kg and  kg, and the cutting force coefficients which are and . 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 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.

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.

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 (51405343, 51505336, and 11702192), Tianjin Research Program of Application Foundation and Advanced Technology (15JCQNJC05000 and 15JCQNJC05200), Innovation Team Training Plan of Tianjin Universities and Colleges (TD12-5043), and Tianjin Science and Technology Planning Project (15ZXZNGX00220).