Abstract
The fulldiscretization method introduced by Ding has high computational efficiency without loss of any numerical precision regardless of low and high radial depth of cut for the stability analysis of ordinary milling operation in time domain. Contrast to ordinary milling, it exists the separation phenomenon between each cutter tooth and the workpiece in feed ultrasonic milling operation. So before analyzing the system stability by the fulldiscretization method, the contact state between the cutter tooth and the workpiece should be made clear. In this paper, the trajectory method is proposed by the authors to judge the contact state between the cutter tooth and the workpiece. The stability analysis is developed for separated feed ultrasonic and ordinary end milling operations by the above methods and is validated via time domain simulations and experiments for both operations. The analyses show excellent agreement with both the time domain simulations and the experiments. Further, several end milling experiments were conducted that demonstrate the result ultrasonic vibration in feed direction can suppress chatter in machining operations.
1. Introduction
High speed milling technology with high production efficiency, high machining accuracy, high machining surface quality and low cost advantages, has gradually become the mainstream of machining technology, widely used in aviation, aerospace, mold, automotive and other fields [1].
However, in traditional machining technology the cutting vibration displacements between the last time and this time have produced phase difference, thus resulting in the occurrence of machining chatter. The chatter reduces the surface machining accuracy, increases tool wear, and shortens the life of the spindle and so on. So, it is one of the main factors of influencing the system stability. Therefore, it is particularly important to research the problem of machining chatter.
Now there are lots of chatter control methods of which basically have ultrasonic machining method, variable speed cutting method, changing tool geometric angles method and active control laws, intelligent control and so on [2]. Ultrasonic machining method is a technology of combination ultrasonic machining with traditional cutting way. It can improve the traditional cutting machining performance. In recent years, scholars have researched into ultrasonic turning technology and ultrasonic grinding technology deeply and widely. Because the trajectory of milling cutter tooth relative to the workpiece is very complex, the study on ultrasonic vibration assisted milling is little and less on stability problems of ultrasonic vibration assisted milling.
Separated feed ultrasonic milling is continuous separation system between the cutter tooth and the part because of applying ultrasonic vibration in feed direction at low radial depth of cut. For stability analysis of milling system at low radial depth, the scholars made greater efforts to seek more accurate and timesaving methods. Altintas and Budak described the dynamic milling process into delayed differential equations including time periodic matrix considering the regenerative effect [3, 4]. They solved the milling stability by the method that is based on Fourier series called Zeroth Order Approximation (ZOA). But, the ZOA method lacks enough accurate to predict stability at low radial immersions. In response to the shortcomings of ZOA, the method of MultiFrequency Solution including of higher order harmonics of periodic coefficients was developed by them [5, 6]. After the ZOA method, time domain methods that compute milling stability from eigenvalue analysis of systems’ Floquet transition matrix gained more attention by the scholars. Insperger and Stepan proposed semidiscretization method to solve stability analysis of milling [7, 8]. Henninger and Eberhard improved computational efficiency and accuracy of the semidiscretization method for periodic delaydifferential equations [9]. Ding Ye and coworkers introduced a new time domain method called fulldiscretization that is based on the direct integration scheme for prediction of milling stability [1013]. This method has high computational efficiency without loss of any numerical precision regardless of low and high radial depth of cut. Indepth application and study of the fulldiscretization method in milling chatter stability analysis can be seen in the works [1421]. Ozoegwu and coworkers found that accuracy of stability result rose with order of fulldiscretization method peaking at fourth order then declined at the fifth order [22]. They extended the accuracy of this method by proposal of higher other vector numerical integration schemes in the formulation of discrete maps [23]. Iglesias and Munoa proposed new analytical formulae related to the parameter domains of both Hopf and period doubling type stability boundaries emerging in the regenerative mechanical model of time periodical milling processes. These formulae are useful to enrich and speed up the currently used numerical methods [24].
In the above of stability analysis, the majority ignored the effects of feed on milling stability. The facts prove that because the contact angle between the cutter tooth and the workpiece is relatively small at low radial depth of cut, the static chip thickness due to feed per tooth generating relative to the dynamic chip thickness can't be ignored. As for milling system, applying ultrasonic vibration in feed direction is the equivalent of feed per tooth of milling system in periodic change, which will lead to cutting separation phenomenon and change of static and dynamic chip thickness. So, on the discussion of stability of ultrasonic vibration assisted milling, the influence of feed per tooth can't be ignored.
This paper has proposed trajectory method to judge cutting separation or not, found out the real instantaneous chip thickness, established time delay differential equations of feed ultrasonic milling system, used the fulldiscretization method to study this system stability, realized numerical simulation by MATLAB7.1 software, obtained this system stability lobe diagrams, and compared with ordinary milling.
2. The motion analysis of feed ultrasonic milling
According to velocity coefficient $K$, the feed ultrasonic milling operation is divided into separated ($K<1$) and unseparated ($K\ge 1$) feed ultrasonic milling. The velocity coefficient is a key parameter. It is shown as follows:
where $v$ is the instantaneous cutting velocity of cutter tooth relative to workpiece, ${v}_{c}$ is the critical velocity, $f$ is ultrasonic frequency, $A$ is ultrasonic amplitude in cutting direction.
In fact, it is not sufficient for only depending on the velocity coefficient to judge cutting separation or not. On the basis of [25] and the results of numerical simulation, the necessary conditions of cutting separation are: $\lambda /N$ is noninteger and $A>0.5{f}_{z}$ (where $\lambda =\omega /{\omega}_{r}$, $\omega $ is ultrasonic transducer vibration angular frequency, ${\omega}_{r}$ is spindle rotational speed of the cutter, $N$ is the number of cutter teeth, ${f}_{z}$ is the feed per tooth.
According to the motion of cutter tip in feed ultrasonic milling, trajectory equations can be expressed by:
where ${v}_{f}$ is feed velocity, $r$ is the cutter radius. The motion simulations of two subsequent cutter teeth in feed ultrasonic milling operation have been acquired by Eq. (2). The results are shown in Fig. 1. It can be clearly seen that the trajectories of two subsequent teeth ($j$ tooth and $j1$ tooth) are continuous intersection, which demonstrate that the current tooth and the workpiece are in state of continuous separation under high frequency.
Fig. 1Milling trajectory simulations of two subsequent cutter teeth
3. Dynamic model of the feed ultrasonic milling
For example end milling, the milling operation is an interaction between the cutter tooth and the workpiece. This is shown in the block diagram of Fig. 2. Considering the vibrations in $x$ and $y$ directions, the model of tool vibrations can be simplified to a massdampingspring mechanical system of two degrees of freedom. The dynamical equations are given by:
where ${m}_{ti}$, ${c}_{ti}$ and ${k}_{ti}$ ($i=x,y$) are the system mass, damping and stiffness in $x$ and $y$ directions respectively, ${F}_{tx}\left(t\right)$ and ${F}_{ty}\left(t\right)$ are the cutting forces on tool in $x$ and $y$ directions at $t$ time, $x\left(t\right)$ and $y\left(t\right)$ are the displacement of the cutter in $x$ and $y$ directions at $t$ time.
The cutting forces on the cutter tooth in $x$ and $y$ directions can be expressed:
where ${\varphi}_{j}\left(t\right)$ is the rotation angel of the cutter tooth $j$at $t$ time: ${\varphi}_{j}\left(t\right)=\frac{2\pi \mathrm{\Omega}}{60}(j1)\frac{2\pi}{N}$, ($j=\mathrm{1,2},...N$), $\mathrm{\Omega}$ is the spindle speed, ${F}_{tj}$ and ${F}_{rj}$ are the tangential and radial forces on the $j$ cutter tooth ,respectively.
The tangential and radial forces on the cutter tooth $j$ are given by:
where ${K}_{t}$ and ${K}_{r}$ are tangential and radial cutting parameters respectively, ${a}_{p}$ is axial cutting depth, ${K}_{te}$ and ${K}_{re}$ are tangential and radial of the cutting edge parameters, $q$ is the cutting constant, ${h}_{j}\left(t\right)$ is the dynamic chip thickness considering regeneration effect, $g\left({\varphi}_{j}\right(t\left)\right)$ is the function that judge whether a tooth is in or out of cut. The tooth is in cut if ${\varphi}_{st}\le {\varphi}_{j}\left(t\right)\le {\varphi}_{ex}$.
This function is given by:
where ${\varphi}_{st}$ and ${\varphi}_{ex}$ are the start and exit angles, respectively. According to the coordinate frame shown in Fig. 2, for up milling, ${{\varphi}_{s}}_{t}=0$ and ${\varphi}_{ex}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}(12{a}_{D})$; for down milling, ${\varphi}_{st}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}(12{a}_{D})$ and ${\varphi}_{ex}=\pi $, where ${a}_{D}$ is the radial immersion ratio.
The instantaneous chip thickness is composed of static and dynamic chip thickness. ${h}_{jst}$ is static chip thickness that the system feed and ultrasonic amplitude in radial direction result in, ${h}_{jdn}$is the dynamic chip thickness that dynamic displacements of two subsequent teeth result in.
The instantaneous chip thickness, static chip thickness and dynamic chip thickness can be expressed as follows.
where $\tau $ is the cutter tooth passing interval, $w\left(t\right)$ is window function that judges whether the interact between the cutter tip and the workpiece or not. Its value is one when the cutter tip contacts with the workpiece, or else its value is zero.
Fig. 2The dynamic model of feed ultrasonic vibration milling
The model of feed ultrasonic milling forces adopts exponential model. To predict the system stability, using linear theory of nonlinear periodic function [26], the Eq. (5) can be linearized around $x=0$ and $y=0$ which results in $F=F(x=0,y=0)+\mathrm{\Delta}F$. Then the linearized force $\mathrm{\Delta}F$ can be written as:
with the matrix ${K}_{h}\left(t\right)$ defined by:
where $u\left(t\right)={\left[\begin{array}{ll}x\left(t\right)& y\left(t\right)\end{array}\right]}^{T}$, $s=\mathrm{s}\mathrm{i}\mathrm{n}\left({\varphi}_{j}\right(t\left)\right)$ and $c=\mathrm{c}\mathrm{o}\mathrm{s}\left({\varphi}_{j}\right(t\left)\right)$.
So, the dynamic model of feed ultrasonic milling can be expressed by:
where $M$, $C$ and $K$ are the mass, damping and stiffness matrices, respectively.
4. Stability prediction of the milling system
4.1. The trajectory method
Fig. 1 shows local amplification trajectories of two subsequent cutter teeth in feed ultrasonic milling system. It can be seen that the instantaneous chip thickness is the length along radial direction of the two subsequent teeth trajectories. In Fig. 1, the instantaneous chip thickness is shown in red solid lines, for example ${P}_{i,j1}{P}_{i,j}$. Because of continuous separation between the current cutter tooth and the workpiece in machining, their contact state should be judged. If separation, instantaneous chip thickness doesn’t exist and loses the sense. Here, on the basis of fulldiscretization method, the authors put forward the trajectory method to judge actual instantaneous chip thickness. The window function is given by:
where ${x}_{i,j}$, ${y}_{i,j}$_{}and ${x}_{i,j1}$, ${y}_{i,j1}$ are the coordinate values in $x$ and $y$ directions of the $j$ and $j1$ teeth at $i$ and $i\tau $ time, respectively. So, it is convenient to use the trajectory method for predicting stability of feed ultrasonic milling system.
4.2. Stability analysis of feed ultrasonic milling system
Substitute $v\left(t\right)=M\dot{u}+Cu/2$ and $q\left(t\right)={\left[\begin{array}{ll}u\left(t\right)& v\left(t\right)\end{array}\right]}^{T}$ into Eq. (11), it is given by:
where:
$M=\left[\begin{array}{cc}{m}_{tx}& 0\\ 0& {m}_{ty}\end{array}\right],C=\left[\begin{array}{cc}2{m}_{tx}{\zeta}_{x}{\omega}_{nx}& 0\\ 0& 2{m}_{ty}{\zeta}_{y}{\omega}_{ny}\end{array}\right],K=\left[\begin{array}{cc}{m}_{tx}{\omega}_{nx}^{2}& 0\\ 0& {m}_{ty}{\omega}_{ny}^{2}\end{array}\right],$
with ${m}_{ti}$, ${\zeta}_{i}$ and ${\omega}_{ni}$ ($i=x,y$) the model mass, the relative damping and the angular natural frequency. ${h}_{xx}$, ${h}_{xy}$, ${h}_{yx}$ and ${h}_{yy}$ are four projections of the specific cutting force coefficient defined as:
${h}_{xy}=w\left(t\right)\sum _{j=1}^{N}g\left({\varphi}_{j}\left(t\right)\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\varphi}_{j}\left(t\right)\right)\left({K}_{t}\mathrm{c}\mathrm{o}\mathrm{s}\left({\varphi}_{j}\left(t\right)\right)+{K}_{r}\mathrm{s}\mathrm{i}\mathrm{n}\left({\varphi}_{j}\left(t\right)\right)\right),$
${h}_{yx}=w\left(t\right)\sum _{j=1}^{N}g\left({\varphi}_{j}\left(t\right)\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\varphi}_{j}\left(t\right)\right)\left({K}_{t}\mathrm{s}\mathrm{i}\mathrm{n}\left({\varphi}_{j}\left(t\right)\right)+{K}_{r}\mathrm{c}\mathrm{o}\mathrm{s}\left({\varphi}_{j}\left(t\right)\right)\right),$
${h}_{yy}=w\left(t\right)\sum _{j=1}^{N}g\left({\varphi}_{j}\left(t\right)\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\varphi}_{j}\left(t\right)\right)\left({K}_{t}\mathrm{s}\mathrm{i}\mathrm{n}\left({\varphi}_{j}\left(t\right)\right)+{K}_{r}\mathrm{c}\mathrm{o}\mathrm{s}\left({\varphi}_{j}\left(t\right)\right)\right).$
Combination the trajectory method above with the fulldiscretization method, the transition matrix of feed ultrasonic milling system has been achieved. Now, the stability investigation is reduced to the problem, whether the eigenvalues of the transition matrix are in modulus less than one by Floquet theory. The compute parameter is chosen as 40. The radial immersion ratio is 0.05. The stability charts are calculated over a 400×200 sized grid of parameters. Computer programs of the proposed approach are all written in MATLAB7.1 and implemented on a personal computer [Intel Core (TM) i53470CPU, 3.2 GHz, 4 GB]. The compute time is about 276.294 seconds.
5. Test verification and results analysis
The chatter tests had been made on VMC850 CNC machining center. The end milling with 4 cutter teeth is solid carbide, the cutter diameter is 8mm, helical angle is 30°. For each experiment the part material is Ti6Al4V, $E=$110 GPa, $\rho =$4505 kg/m^{3}, its size is 30 mm×10 mm×15 mm. The test equipment has Kistler 9257B three component dynamometer, NI PCI9205 data acquisition analyzer, Kistler 5070 charge amplifier, 8207 impact hammer and the piezoelectric acceleration sensor made in Orient Institute. The model parameters of this system have been achieved by impact hammer test using the above device as shown in Fig. 3. These parameters are shown in Table 1. The coefficients of cutting force have been solved by the method in [27]. The coefficients of cutting force are given as follows: ${K}_{t}=$ 9.23×10^{8} N/m^{2}, ${K}_{r}=$4.95×10^{8} N/m^{2}, ${K}_{te}=$12.64 N/m, ${K}_{re}=$ 7.95 N/m.
Table 1Model parameters
$M$ (kg)  $C$ (Ns/m)  $K$ (N/m)  
$\left[\begin{array}{cc}0.06641& 0\\ 0& 0.06586\end{array}\right]$  $\left[\begin{array}{cc}1.758& 0\\ 0& 1.769\end{array}\right]$  $\left[\begin{array}{cc}4.562\times 1{0}^{5}& 0\\ 0& 4.435\times 1{0}^{5}\end{array}\right]$ 
Fig. 3Photograph of the setup of modal impact tests
According to the above methods, substitute model parameters and ultrasonic frequency 19.58 KHz into Eq. (13), the stability lobe diagrams are acquired at the different ultrasonic amplitude and shown in Fig. 4. It can be seen that the stability areas of feed ultrasonic milling system are more than that of ordinary milling. As the ultrasonic amplitude increases, the limiting axial depth of cut increases. At about 3200 rpm, the limiting axial depth of cut in stability diagram of ultrasonic amplitude 12 μm increases 15.56 percent compared with that of ordinary milling. However, the increase is very slight since the influence of ultrasonic amplitude as a part of feed on cutting forces is slight. But, the different ultrasonic amplitude can change the contact state of the toolworkpiece.
Fig. 4Stability lobe at different ultrasonic amplitude
Chatter tests had been conducted as Fig. 5(a) and (b) show. The proving results are shown in Fig. 6. The signs □, ○ and ⊗ are expressed by chatter, stability and slight chatter in the feed ultrasonic milling, respectively. In these experiments, ultrasonic frequency and ultrasonic amplitude are 19.58 KHz and 8 μm. Chatter is determined by examining the force signals and their frequency spectral plots. Tested fivepoint A, B, C, D, E, they are instability point under ordinary milling but stability point under feed ultrasonic milling. Two examples of the points B and D under ordinary and feed ultrasonic milling respectively are shown in Figs. 710. At the point of B, the spindle speed, the radial depth and feed are 1600 rpm, 0.2 mm, and 0.04 mm, respectively. Highfrequency variations are evident in the force plot. The frequency spectral plot shows significant energy in the $Y$ force plot at 128.7 Hz and 235.4 Hz. Since these frequencies are not the tool rotation frequencies or a multiple of the tool rotation frequencies, chatter was present. So, the process of ordinary milling at point B is not stable. But, under feed ultrasonic milling as Fig. 8 shows, there are no highfrequency variations in the force plots. These frequencies that the frequency spectral plot shows significant energy are the tool rotation frequencies and multiples of the tool rotation frequency. So, the operation at point B is stable. At the point of D, the spindle speed, the radial depth and feed are 2500 rpm, 0.4 mm, and 0.06 mm, respectively. In Fig. 9, the frequency spectral plot shows significant energy in the $Y$ force plot at 1156.2 Hz and 1654.67 Hz. Chatter was present in this operation. So, the operation at point D is not stable under ordinary milling. The same experiment at point D under feed ultrasonic milling had been done and the results are shown in Fig. 10. It can be seen that the operation is stable. The frequency spectral plot shows significant energy at the tool rotation frequencies or a multiple of the tool rotation frequencies. In addition, the pictures of machining surface under ordinary milling and feed ultrasonic milling are shown in Fig. 11. It can be seen that the machining surface at feed ultrasonic milling is more smooth than that of ordinary milling. Under the condition of ordinary milling, the machining surface of the workpiece presents chatter waves. So, the analysis and experimental results demonstrate that under the same machining parameters the feed ultrasonic milling operation is more stable than that of ordinary milling operation.
Test results agree with the simulation results very well by chatter tests. The analysis and experimental results indicate that for all spindle speeds the limiting depth of cut increases as the ultrasonic amplitude increases. The results demonstrate that ultrasonic vibration in feed direction can dramatically suppress the presence of chatter for this end milling operation. The results provide empirical evidence that ultrasonic vibration assisted end milling can improve the stability of the system.
Fig. 5The setup for the chatter tests
a) Photograph of processing site of the chatter tests
b) Schematic of the setup for the chatter tests
Fig. 6Stability comparison between prediction and experiment results
Fig. 7Timedomain signals and frequency spectrum of cutting forces of point B under normal milling
a) Timedomain signals of cutting forces
b) Frequency spectrum of cutting forces
Fig. 8Timedomain signals and frequency spectrum of cutting forces of point B under feed ultrasonic milling
a) Timedomain signals of cutting forces
b) Frequency spectrum of cutting forces
Fig. 9Timedomain signals and frequency spectrum of cutting forces of point D under normal milling
a) Timedomain signals of cutting forces
b) Frequency spectrum of cutting forces
Fig. 10Timedomain signals and frequency spectrum of cutting forces of point D under feed ultrasonic milling
a) Timedomain signals of cutting forces
b) Frequency spectrum of cutting forces
The reasons of individual errors are: one is the modal parameters without considering the coupling of $x$ and $y$ directions; two is the ultrasonic amplitude that can’t guarantee constant in machining, three is to assume that the cutting force coefficients unchanged in the process. However, the above errors do not affect the overall stability trend which is achieved by the 4.1 methods. Because the stability limit has uncertainty, actually the stability limit should be stability belt [28]. So, it is feasible and reliable for applying trajectory method to stability prediction of the separated feed ultrasonic milling system.
It can be also seen from Fig. 6 that with the increase of speed, the limiting axial depth of cut of the slight chatter point is gradually increasing. The higher the speed, the higher is the stability probability of the test points on the border of the stability prediction diagram. With the increase of speed, the cutter tooth passing frequency rises and the time of the cutter tooth left on the workpiece lengthens, which is convenient for the full contact of the cutter tooth and the workpiece and add the instantaneous chip thickness. In addition, ultrasonic vibration that is the varying accelerated motion, makes surrounding contact materials of the workpiece and cutting tool quickly shattered and is conducive to chip shed from the workpiece and cutting tool.
Fig. 11Workpiece surface quality after milling
a) Unstability
b) Stability
6. Conclusions
The fulldiscretization method introduced by Ding [10] was extended in this paper to consider the effect of ultrasonic vibration in feed direction on end milling operation. Trajectory method has been put forward based on motion characteristics of feed ultrasonic milling. The method was compared to time domain simulation and experimental results for end milling operations. The analyses and experiments results clearly demonstrate that the trajectory method is a simple, practical and high computational efficiency method. It opens up a new way of judging the stability of intermittent cutting.
The analyses provide insight into stability of the separated feed ultrasonic milling system. The results demonstrate that stability region of feed ultrasonic milling system compared to ordinary milling system improves. With the increase of ultrasonic amplitude, stability region is continuous improving, but improving slight. The stability belts have been achieved at different ultrasonic amplitudes. The analyses developed for feed ultrasonic end milling operations compared very well to both time domain simulation and experimental results. The results show that chatter was suppressed under ultrasonic vibration and with the increase of ultrasonic amplitude, the effect of suppressing chatter is better.
References

Ai X. High Speed Machining Technology. National Defense Industry Press, Beijing, 2003.

Tang J., Zhao B. Stability analysis of the separated longitudinaltorsional composite ultrasonic milling. Acta Armamentarii, Vol. 36, Issue 7, 2015, p. 13191325.

Altintas Y., Budak E. Analytical prediction of stability lobes in milling. Annals of the CIRP Journal of Manufacturing Science and Technology, Vol. 44, Issue 1, 1995, p. 357362.

Altintas Y. Analytical prediction of threedimensional chatter stability in milling. Japan Society of Mechanical Engineers, International, Vol. 44, Issue 3, 2001, p. 717723.

Merdol S. D., Altintas Y. Multi frequency solution of chatter stability for low immersion milling. Journal of Manufacturing Science and EngineeringTransactions of the ASME, Vol. 126, Issue 3, 2004, p. 459466.

Altintas Y., Stepan G., Merdol D., Dombovari Z. Chatter stability of milling in frequency and discrete time domain. CIRP Journal of Manufacturing Science and Technology, Vol. 1, Issue 1, 2008, p. 3544.

Insperger T., Mann B. P. Stability of upmilling and downmilling, Part 1: alternative analytical methods. International Journal of Machine Tools and Manufacture, Vol. 43, Issue 1, 2003, p. 2534.

Insperger T., Stepan G. Updated semidiscretization method for periodic delay differential with discrete delay. International Journal for Numerical Methods in Engineering, Vol. 61, Issue 1, 2004, p. 117141.

Henninger C., Eberhard P. Improving the computational efficiency and accuracy of the semidiscretization method for periodic delaydifferential equations. European Journal of Mechanics A – Solids, Vol. 27, Issue 6, 2008, p. 975985.

Ding Y., Zhu L. M., Zhang X. J., et al. A fulldiscretization method for prediction of milling stability. International Journal of Machine Tools and Manufacture, Vol. 50, Issue 5, 2010, p. 502509.

Ding Y., Zhu L. M., Zhang X. J., Ding H. Secondorder fulldiscretization method for milling stability prediction. International Journal of Machine Tools and Manufacture, Vol. 50, Issue 10, 2010, p. 926932.

Ding Y., Zhu L. M., Zhang X. J., Ding H. Numerical integration method for prediction of milling stability. Journal of Manufacturing Science and Engineering, Transactions of the ASME, Vol. 133, Issue 3, 2011, p. 3100531009.

Zhang X. J., Xiong C. H., Ding Y., et al. Milling stability analysis with simultaneously considering the structural mode coupling effect and regenerative effect. International Journal of Machine Tools and Manufacture, Vol. 53, Issue 1, 2012, p. 127140

Ding Y., Zhu L. M., Zhang X. J., Ding H. Stability analysis of milling via the differential quadrature method. Journal of Manufacturing Science and Engineering, Transactions of the ASME, Vol. 135, Issue 4, 2013, p. 044502.

Liu Y. L., Zhang D. H., Wu B. H. An efficient full discretization method for prediction of milling stability. International Journal of Machine Tools and Manufacture, Vol. 63, 2012, p. 4448.

Insperger T. Fulldiscretization and semidiscretization for milling stability prediction: some comments. International Journal of Machine Tools and Manufacture, Vol. 50, 2010, p. 658662.

Guo Q., Sun Y. W., Jiang Y. On the accurate calculation of milling stability limits using thirdorder fulldiscretization method. International Journal of Machine Tools and Manufacture, Vol. 62, 2012, p. 6166.

Ozoegwu C. G., Omenyi S. N. Wave attenuation effects on the chatter instability of endmilling. Noise Control Engineering Journal, Vol. 61, Issue 4, 2013, p. 436444.

Huang T., Zhang X. M., Zhang X. J., Ding H. An efficient linear approximation of acceleration method for milling stability prediction. International Journal of Machine Tools and Manufacture, Vol. 74, 2013, p. 5664.

Ozoegwu C. G. Least squares approximated stability boundaries of milling process. International Journal of Machine Tools and Manufacture, Vol. 79, 2014, p. 2430.

Ding Y., Niu J. B., Zhu L. M., Ding H. Numerical integration method for stability analysis of milling with variable spindle speeds. Journal of Vibration and AcousticsTransactions of the ASME, Vol. 138, 1, p. 201611010.

Ozoegwu C. G., Omenyi S. N., Ofochebe S. M. Hyperthird order fulldiscretization methods in milling stability prediction. International Journal of Machine Tools and Manufacture, Vol.92, Issue.3, 2015, p.1 9.

Ozoegwu C. G. High order vector numerical integration schemes applied in state space milling stability analysis. Applied Mathematics and Computation, Vol. 273, 2016, p. 10251040.

Iglesias A., Munoa J., Ciurana J., et al. Analytical expressions for chatter analysis in milling operations with one dominant mode. Journal of Sound and Vibration, Vol. 375, Issue 4, 2016, p. 403421.

Shen X. H. Study on the Technology and Mechanism of Ultrasonic Vibration Assisted Milling. Shandong University, Shandong, 2011.

Faassen R. P. H., van de Wouw N., Oosterling J. A. J., et al. Prediction of regenerative chatter by modeling and analysis of highspeed milling. International Journal of Machine Tools and Manufacture, Vol. 43, 2003, p. 14371446.

Altintas Y. Manufacturing Automation, Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design. Cambridge University Press, New York, 2000.

Mohammad H. K. Robust Multicriteria Optimization of Surface Location Error and Material Removal Rate in Highspeed Milling under Uncertainty. University of Florida, Gainesville, 2005.
About this article
The authors gratefully acknowledge the financial support of the National Science Foundation of China (No. 51475148 and 51175153). The authors also would like to thank the anonymous reviewers for the valuable comments and suggestions to improve the manuscript.