Stability of Milling Process with Variable Spindle Speed Using Runge–Kutta-Based Complete Method

The variable-spindle-speed (VSS) technique is effective in preventing regenerative chatter in milling processes. However, spindle-speed-modulation parameters should be deliberately selected to augment the material removal rate. Stability-prediction algorithms of stability predicting play an important role in this respect, as they allow the prediction of stability for all ranges of a given spindle speed. The increase in calculation time in variable-spindle-speed milling, which is caused by the modulation frequency, hinders its practical use in the workshop. In this paper, a Runge–Kutta-based complete discretization method (RKCDM) is presented to predict the stability of milling with variable spindle speeds, which is described by a set of delay differential equations (DDEs) with time-periodic coefficients and time-varying delay. The convergence and calculation efficiency are compared with those of the semidiscretization method (SDM) under different testing configurations and milling conditions. Results show that RKCDM is more accurate and saves at least 50% of the calculation time of SDM. The effects of modulation parameters on the stability of VSS milling are explored through stability lobe diagrams produced from RKCDM.


Introduction
In the milling process, chatter caused by vibrations between the cutter and workpiece hinders manufacturing production. It cannot only reduce machining performance, but can cause a poor workpiece-surface finish, increase tool wear, and decrease the lifetime of the machine tool itself.
Chatter is generally classified in two categories. Primary chatter is usually triggered by friction between the cutter and the workpiece. Secondary chatter, which is caused by a regenerative mechanism, is most common factor. at is why many publications discuss only regenerative chatter.
Several strategies have been proposed to suppress or avoid chatter. ey can be roughly categorized as follows: (1) Strategies based on a stability lobe diagram (SLD) and the selection of cutting parameters (2) Active chatter-suppression techniques to decrease the vibration of a machine tool system (3) Semiactive chatter-suppression strategies (4) Passive chatter-suppression strategies accounting for a machine tool's compliance or configuration (5) Other strategies based on signal surveillance in the milling process for the adjustment of milling parameters SLD is the most logical and widely accepted method to select cutting parameters to avoid chatter. Researchers have proposed numerous methods to predict machining stability.
Because regenerative chatter is a result of improper phase differences between undulate surfaces caused by successive cuts by a tool, alternative methods usually vary the time lag between cutter teeth. Slavicek [13] demonstrated that this can be achieved through variable-pitch cutters. Altintas et al. [14], Olgac and Sipahi [15], Jin et al. [16], and Sims et al. [17] all demonstrated that variable-pitch cutters can avoid or significantly suppress chatter. e idea of changing the spindle speed to suppress or avoid chatter in machining dates back to the 1970s and the pioneering research of Stoferle and Grab [18]. After that, Takemura et al. [19] studied the stability of variable-spindlespeed turning, and experimental tests showed only small improvements in the depth of cuts, although the stability of VSS was higher than that of constant-spindle-speed (CSS) milling. Altintas and Chan [20] investigated the stability of VSS with sinusoidal modulation using a time-simulation method. Tsao et al. [21] described the VSS milling system in the angle domain and showed that chatter could be suppressed or eliminated to a great extent by properly selecting the modulation parameters. Long and Balachandran [22] analyzed the stability of up-and down-milling with sinusoidal modulation using the semidiscretization method. In the frequency domain, Zatarain et al. [23] proposed a general theory to analyze the stability of VSS milling. Seguy et al. [24,25] used a semidiscretization method to study the stability of VSS milling in the high-speed domain and reported that VSS was an effective method to reduce perioddoubling chatter. Based on the full-discretization method, Xie and Zhang [26,27] presented an improved SDM to predict the stability of VSS milling. Ding et al. [28] studied the stability of the VSS milling system using the constantstep numerical method. Niu et al. [29] unified different periodic time-variant spindle-speed-modulation schemes in a unique framework, derived and calculated the timevarying delay based on the fourth-order Runge-Kutta method, and calculated and compared the stability of VSS milling with different modulation schemes using variablestep numerical integration. Totis et al. [30] studied the stability of milling with VSS in the angular domain. Jin et al. [31] studied the stability of milling with constant pitch, variable pitch, and VSS, using a full-discretization method.
As mentioned above, much research has been done to predict the stability of VSS milling. However, the required computational time increases greatly due to frequency modulation, and this hinders the application of VSS [29]. Niu et al. [32] and Li et al. [11] recently proposed a Runge-Kutta-based complete discretization method (RKCDM) using complete discretization scheme [33], but their research was based exclusively on CSS machining. In this paper, we present a two-degree-of-freedom dynamic milling model expressed by DDEs and a Runge-Kutta-based complete discretization method for efficient analysis of VSS milling. e convergence and computational efficiency are compared with a state-of-the-art semidiscretization method. Finally, we discuss the VSS modulation effects on the stability of milling with the aid of maps and contour plots of axial depth of cut (ADOC) calculated by RKCDM.

Modeling of VSS Process
2.1. VSS. In this paper, periodic spindle-speed variation with sinusoidal modulation will be considered and analyzed. As shown in Figure 1, a sinusoidal spindle-speed variation Ω(t) is periodic at a time period T � 60/Ω 0 /RVF, with a nominal value Ω 0 , and an amplitude of the speed variation Ω A � RVA × Ω 0 . e shape function of the sinusoidal modulation is modeled as where RVA � Ω A /Ω 0 is the ratio of the speed-variation amplitude to the nominal spindle speed, RVF � 60/(Ω 0 T) is the ratio of the speed-variation frequency to the nominal spindle speed, and η is the phase angle of the sinusoidal modulation at time t � 0.

Structural Model.
In the machining process, relatively more flexible parts will undergo dynamic deflection and generate an undulating surface. So, the system can be simplified to an oscillator motivated by the force of milling.
We show an example model [34] in Figure 2. e governing equation of the model is where ζ is the relative damping, ω n is the angular natural frequency, w is the axial depth of cut, and m t is the modal mass of the teeth. e dynamic characteristic is assumed to be identical in the x and y directions, so the model mass m t , the relative damping ζ, and the natural frequency ω n are the same in both directions. e time-varying cutting force coefficients, h xx (t), h xy (t), h yx (t), and h yy (t), are expressed as 2 Mathematical Problems in Engineering where N is the number of cutter teeth, and K t and K n are, respectively, the tangential and normal linearized cuttingforce coefficients. e j-th tooth angular position ϕ j (t) is defined as e window function g(ϕ j (t)) can be defined as where ϕ st and ϕ ex are, respectively, the start and exit angles of the j-th cutter tooth. For up-milling, ϕ st � 0 and ϕ ex � arccos(1 − 2 · aD). For down-milling, ϕ st � arccos(1 − 2 · aD) and ϕ ex � 0, where aD � r/D is the radial immersion ratio, r is the radial depth of cut, and D is the cutter diameter.

Calculation of Variable Time Delays.
Substituting (1) into (4), the angular position of the j-th tooth can be expressed as where ω m � 2π/T is the angular velocity. In this case, the time delay τ(t) cannot be given in an explicit form, as it is determined by an equation: To calculate τ(t) directly is time-consuming [28], so we simplify the integral equation. By substituting (1) into (7), we achieve an implicit relationship: e function τ(t) cannot be solved in an explicit form; hence, it must be calculated numerically. However, if RVA and RVF are small enough; then, τ(t) is approximately [22] where τ 0 � 60/(NΩ 0 ) and τ 1 /τ 0 � Ω 1 /Ω 0 . As shown in Figure 3, the exact time delays which is calculated using numerical solution from (8) and the approximate values from (9) is consistency. So, the approximation in (9) is effective. Figures 3(a) and 3(b) show that the approximated error for RVA � 0.1 is less than RVF � 0.2 because a small value of RVA is a condition of (9). Below, we will calculate time-varying delay through (9).

Stability of Milling with VSS
It must be pointed out that although VSS with sinusoidal modulation is discussed here, the proposed RKCDM can deal with other modulation types of VSS with no further difficulty. rough algebraic transformation, (2) can be written as

Mathematical Problems in Engineering 3
e first step of discretization is construction of the timeinterval discretization of t i , t i+1 with length Δt, i � 0, 1, . . ., such that T � kΔt, where k is an integer that can be considered an approximation parameter referring to the time period. e time-step Δt can be expressed as e average delay for the discretization interval t i , t i+1 is defined as e series of integers m i is introduced as follows: since and where int ( * ) is the function that rounds positive numbers toward zero (e.g., int (7.875) � 7). Since the delay varies periodically over time, integer m i might vary for different semidiscretization steps.
Note that M can be considered an approximation parameter regarding the length of the time delay. e most popular version of the Runge-Kutta algorithm is [35] x i+1 � x i + 1 6 where x i denotes x(iΔt), x t+1 denotes x((i + 1)Δt), and t i represents iΔt, i ∈ Z and 0 ≤ i ≤ m i . e following equations can be deduced in sequence to solve (10): where By substituting (19)- (21) into (18), we can obtain the following iterative formula: where the coefficients can be expressed as where I is an n × n identity matrix and n is the dimension of vector x(t).
To obtain the transition matrix, we first define a new n × (M + 1) dimensional vector z i as Finally, the resulting discrete map is expressed as where the coefficient matrix D i can be constructed as a (2M + 4)-dimensional matrix: erefore, the approximate Floquet transition matrix can be described as en, the stability of the system can be calculated based on the Floquet theory. If the moduli of all eigenvalues of the transition matrix Φ are less than unity, then the system is stable; otherwise, it is unstable.

Simulation and Discussion
In this section, an experimental example [3] is used to verify the accuracy and efficiency of RKCDM. e modal parameters of the milling systems are tool diameter D � 12.7, number of teeth N � 2, tangential cutting force coefficient K t � 6 × 10 8 N/m 2 , radial cutting force coefficient K n � 2 × 10 8 N/m 2 , natural frequency ω nx � ω ny � 922.0 Hz, modal damping ratio ζ x � ζ y � 0.011, and modal mass m tx � m ty � 0.03993 kg. All calculations are performed in MATLAB R2014a on a desktop computer (Intel Core (TM) i7-6700K central processing unit, 4.0 GHz, 16 GB). We first discuss the convergence of RKCDM and SDM.

Convergence and Computation Efficiency.
e convergence rates of RKCDM and SDM with different conditions are compared in Figure 4. e discrete number of time delays τ (τ(t) for VSS milling) is set as m for constant-speed milling and M for VSS milling and is determined according to (17), and the maximal module of the eigenvalues of the transition matrix Φ is set as λ. To compare the convergence of RKCDM and SDM, the value of λ under condition of m � 500 is written as λ 0 as the exact value to assess its convergence rate. When m increases, the approximation error of RKCDM is lower than SDM. e same phenomenon can be obtained for VSS milling with the increase of M. is means that RKCDM converges faster than SDM with the increase of the discrete number of time delays, because for SDM, the discretization error is O(h) 2 [36]. However, discretization error for RKCDM is O(h) 5 . Niu et al. [32] used two-point barycentric Lagrange points to approximate the middle points, so the remainder of the

two-point interpolation is
To increase the accuracy of RKCDM, 2k discrete points are taken in one toothpassing interval to calculate cutting-force coefficients h xx (t), h xy (t), h yx (t), and h yy (t) for the following examples. So, at time interval [t i , t i+1 ], A i , A i+1 , B i , and B i+1 , can be calculated using even points, and A i+0.5 and B i+0.5 are calculated using the left points. e benefit of this is that the accuracy of RKCDM is the same with the classical fourthorder Runge-Kutta method. To improve the prediction accuracy, the discrete number M of time-varying delays τ(t) is always greater than 80, and according to Figure 4, RKCDM can have the minimum approximation error.
To compare the computational efficiency, we refer to Figure 5, which shows the time consumption using RKCDM and SDM. For VSS milling, the consumption time to achieve stability increases as RVA decreases. For small RVA, more time can be saved with RKCDM than with SDM. Hence, we conclude that the time consumption increases with the number of discrete time delays. e results of Figure 5 indicate that 54% of time consumption can be saved using RKCDM instead of SDM.  Figure 6 shows the SLD of VSS and CSS milling. Using VSS, the stability boundary obviously increases, especially at low spindle speeds. Although the depth of cut increases using VSS, the modulation also has an important effect on the stability. Taking Figure 6(a) as an example, the depth of cut around nominal spindle speed 22,000 rpm decreases abruptly. To comprehensively explore the effects of modulation parameters on the stability of milling, the relationships among depth of cut, RVA, and RVF are shown in the 3D plots of Figures 7 and 8. According to Figure 6(a), the axial depth of cut has a minimum value of 1 mm for CSS milling.

Analysis of Effect of Modulation Parameters on Stability of VSS Milling.
After VSS is used, the ADOC at these two nominal spindle speeds increases significantly, as illustrated in Figures 7 and 8. Figure 7 shows that the ADOC changes with the variation of sinusoidal modulation parameters RVA and RVF around a nominal spindle speed of 4,900 rpm. Obviously, when VSS is used, the lowest point is higher than that of CSS milling. From Figure 7, we know that the ADOC is more sensitive to RVA than to RVF; so, for the machinist, the selection of modulation parameters is more interesting than that of spindle speed. Taking a nominal spindle speed of 4,900 rpm as an example, with modulation parameters RVA � 0.10 and RVF � 0.17, ADOC will increase to 4.16 mm. e contour plot in Figure 7(b) shows that not only at the vicinity of point (0.10, 0.17) but also all stable ADOC limit in these two figures (Figures 7 and 8  show a three-dimensional stability diagram at a nominal spindle speed of 22,000 rpm in Figure 8. With VSS, ADOC is always higher than CSS in Figure 7. However, this is different at a nominal spindle speed of 22,000 rpm, where ADOC does not always increase with VSS milling over CSS milling. Figure 6 shows that ADOC is 7.83 mm without VSS. In Figure 8, ADOC is always lower than 7.83 mm with RVA in the range [0.2, 0.3]. So, it is important to choose modulation parameters appropriately to improve stability with the VSS technique. In Figure 8, ADOC increases to a local summit with the variation of RVA and RVF, such as around RVA � 0.06 and RVF � 0.16, RVA � 0.08 and RVF � 0.13, and RVA � 0.14 and RVF � 0.15. Taking RVA � 0.06 and RVF � 0.16 as an example, the ADOC can reach 9.99 mm in this area. So, ADOC can still be increased with appropriate selection of modulation parameters.

Conclusion
A new algorithm for stability predictions of VSS milling was developed, and the following conclusions were drawn: (1) e governing DDE of the time-periodic milling process including the regenerative effect is described by a set of algebraic equations, and the stability of the system is predicted by the Floquet theory. (2) Comparisons of convergence rate and computational efficiency were conducted. RKCDM shows a higher convergence rate than SDM, so RKCDM is more accurate. (3) Simulation results indicate that RKCDM can be used for parameter selecting in VSS milling. e maps and contour plots of ADOC produced from RKCDM show the stability of VSS milling with sinusoidal modulation changes with RVA and RVF, and the stability is more sensitive to variation in RVA. e maps and contour plots of ADOC can be used to select the optimal modulation parameters for VSS milling.
Finally, although the optimal parameters of the VSS technique have been selected through the method proposed here, it is obvious that the parameters should be further tested, and this will be our future work.

Data Availability
e data used to support the findings of this study are included within the article. Mathematical Problems in Engineering 9