Milling Stability Prediction with Multiple Delays via the Extended Adams-Moulton-Based Method

The occurrence of machining chatter may undermine the workpiece surface quality, accelerate the tool wear, and even result in serious damage to the machine tools. Consequently, it is of great importance to predict and eliminate the presence of such unstable and detrimental vibration. In this paper, we present an extended Adams-Moulton-based method for the stability prediction of milling processes withmultiple delays. Taking the nonuniform pitch cutters or the tool runout into account, the regenerative chatter for milling operations can be formulated as delay differential equations with multiple delays. The dynamics model for milling regenerative chatter is rewritten in the state-space form. Dividing the spindle rotation period equally into small time intervals, the delay terms are approximated by Lagrange interpolation polynomials, and the Adams-Moulton method is adopted to construct the Floquet transition matrix. On this basis, the milling stability can be derived from the spectral radius of the transition matrix based on Floquet theory. The calculation efficiency and accuracy of the proposed algorithm are verified through making comparisons with the semidiscretization method (SDM) and the enhanced multistage homotopy perturbation method (EMHPM). The results show that the proposed method has both high computational efficiency and accuracy.


Introduction
In machining operations, chatter vibration is still one of the main constraints to high productivity and part quality.It is a typical kind of self-excited vibration between the cutter and the workpiece and can occur in almost every machining process.The onset of such detrimental instability may result in poor surface roughness, rapid tool wear, and large reduction of tool life.Therefore, many researches on the modelling, prediction, and avoidance of milling chatter have been conducted [1][2][3].Theoretically speaking, modecoupling, frictional, thermomechanical, and regenerative mechanisms can induce chatter to occur [4].With regard to milling operations, the regenerative chatter is considered as the most common unstable situation.The corresponding mathematical model can be described by delay differential equations (DDEs) with time-periodic coefficients [5,6].
To achieve stable milling operations, one effective and significant technique is selecting proper cutting parameters based on the stability lobe diagrams, which can be acquired via the milling stability prediction.Therefore, many approaches have been proposed to approximate the DDEs to derive the milling stability lobe diagrams, such as numerical methods [7][8][9][10], analytical methods [11][12][13], and semianalytical methods .Analytical and semianalytical methods have significant advantages over numerical methods on the computation efficiency.Consequently, they are widely used in industrial applications and gained extensive attention from the academic fields.Altintas ¸and Budak [11,12] proposed the zero-order approximation method for the stability prediction of milling operations in the frequency domain.The prominent advantage of this approach is its low computational cost.Nevertheless, the zero-order approximation method is not quite competent for low radial immersion conditions.To improve the zero-order approximation method, Merdol and Altintas [13] presented the so-called multifrequency solution through including more numbers of harmonics.With the help of the weighted residual method, Bayly et al. [14] introduced the temporal finite element analysis (TFEA) method to predict milling stability.Butcher et al. [15] 2 Mathematical Problems in Engineering developed the Chebyshev collocation method for the stability analysis of milling processes.Urbikain and coworkers extended the Chebyshev collocation method for various turning processes, including turning of nonrigid parts [16], turning of rigid parts [17], turning with low rotational speeds [18], and heavy-duty turning processes [19].Insperger et al. [20][21][22] presented the well-accepted semidiscretization methods (0st SDM and 1st SDM) based on Floquet theory for periodic delayed systems.The semidiscretization methods can effectively be utilized for predicting stability of different milling processes [6].Consequently, the SDMs are commonly utilized as benchmarks methods for other timedomain semianalytical methods.To obtain higher accuracy and efficiency, Jiang et al. [23] presented a second-order semidiscretization method by utilizing Newton interpolation polynomials and improved precise time-integration algorithm.Different from the SDMs, Ding et al. [24,25] developed the full-discretization methods (1st FDM and 2nd FDM) for milling stability prediction, in which the state term, the delay term, and the parameter matrix are approximated by linear interpolations, respectively.Three-order and hyperthird-order FDMs were subsequently investigated by Quo et al. [26] and Ozoegwu et al. [27], respectively.However, these methods lead to an increase in computational time since the structures are more complex.Ding et al. [28,29] presented the numerical integration method (NIM) for milling operations with both constant and variable spindle speeds.Niu and coworkers [30] further developed the variable-step numerical integration method (VNIM) for milling stability prediction with periodic spindle speed variation.Different from SDMs and FDMs, Li et al. [31] developed the complete discretization scheme (CDS) by discretizing all parts of DDE and utilizing Euler's method.Xie [32] presented an improved complete discretization method for more efficient stability prediction of milling processes, in which the coefficient matrixes were approximated by linear interpolation.To improve the accuracy of the CDS, Li et al. [33] developed the Runge-Kuttabased complete discretization method.By approximating the term with time-derivative with a weighted linear sum of the corresponding function values, Ding et al. [34] suggested the so-called differential quadrature method (DQM) for the stability analysis of milling operations.Olvera and coworkers [35] combined homotopy method with simulated annealing algorithm for fast prediction of milling stability lobes.Niu et al. [36] proposed the Runge-Kutta methods (CRKM and GRKM).Ding et al. [37] developed the wavelet-based approach for stability analysis of periodic delay differential systems.Zhang et al. [38] suggested Simpson based method (SBM) for the milling stability prediction.Applying the finite difference method and extrapolation method, Zhang et al. [39] presented the numerical differentiation method for stability analysis of milling processes.More recently, Qin et al. [40] presented the Adams-Moulton-based method (AMM) for the stability prediction of milling processes, which has high efficiency and accuracy compared with the 1st SDM and the SBM.
Nevertheless, the above works were mostly conducted based on the ideal milling operations with regular uniform pitch cutters, in which there exists only single time delay.
Taking the nonuniform pitch cutters or the case tool runout into account, the regenerative chatter models for milling operations are described by delay differential equations with multiple delays.As a consequence, many efforts have been made to extend the above algorithms to the multiple delays case.Altintas ¸et al. [41] developed an analytical method for stability prediction of milling process with variable pitch cutters, which was well validated by extensive milling experiments.The results of their research demonstrated the significant influence of the pitch angles on the stability domain.With the aid of the analytical method, Budak [42,43] developed an optimal pitch angles design method for increasing the milling stability.On the basis of cluster treatment of characteristic roots method, Olgac and Sipahi [44] studied the stability boundary of milling process with unequally pitched cutters and presented an optimization procedure for the geometry design of variable pitch cutters.Sims et al. [45] employed the modified SDM, the time-averaged SDM, and the TFEA method to investigate the stability of variable pitch and variable helix milling cutters.Their works showed that under small radial immersions condition the cyclic-fold bifurcations can arise for both nonuniform pitch and variable helix milling tools.Based on the updated semidiscretization method from [21], Wan et al. [46] developed a unified method to predict milling stability with multiple delays arising in variable pitch cutters or cutter runout.Zhang et al. [47,48] developed an improved FDM and a variable-step NIM for the stability prediction of milling operations with multiple delays.Compeán and coworkers [49] developed the enhanced multistage homotopy perturbation method (EMHPM) for milling stability analysis with multiple delays.The so-called spectral element approach was introduced by Khasawneh and Mann [50] for stability analysis of timeperiodic delay equations with multiple delays.Jin et al. [51][52][53] presented an improved SDM to investigate the effect of the tool geometries on the stability trends for variable pitch or variable helix milling.Sims [54] introduced an efficient approach to variable helix tool stability based upon the Laplace transform.Ding and coworkers [55] extended the differential quadrature method for stability analysis of dynamic systems with multiple delays.By combining threeorder FDM and variable interpolation technique, Guo et al. [56] proposed a time-domain semianalytical method for prediction of milling stability lobes with nonuniform helix tools.
In recent years, efficient and accurate milling stability prediction has been a key issue both in academic and industrial fields.However, it is difficult to achieve both high computational accuracy and efficiency simultaneously.Based on our previous work [40], this paper develops an extended Adams-Moulton-based method for the milling stability prediction with multiple delays.The remainder of this paper is organized as follows.After the introduction, Section 2 gives a concise description of the dynamics model for milling operations with multiple delays.The extended Adams-Moulton-based method (EAMM) is proposed in Section 3. Section 4 validates the computation accuracy and efficiency of the proposed method by a two-DOF milling operation.Finally, the conclusion is drawn in Section 5.

Mathematical Model of Milling Operations
Theoretically speaking, when the constant pitch cutter is employed and the cutter runout is neglected, there exists only one delay term in the system dynamics equation.Based on [1,6,21], the dynamics model of milling operations that include the regenerative effect can be modelled by the following delay differential equation: where q() denotes the modal vector of the cutter.M, C, and K represent the modal parameter matrixes.  denotes the depth of cut.K  () is the periodic coefficient matrix, that is, K  () = K  ( + ). is the time delay.For milling operations,  equals the tooth passing period.However, the actual milling cutters cannot be always completely symmetrical.Therefore, there may exist a certain deviation between the spindle rotation axis and the tool geometry axis, which constitutes the so-called cutter runout.In addition, when the variable pitch cutter without cutter runout is considered, any cutting point will always remove the surface left by the first previous tooth.In this case, due to the unevenly pitched space angle, the delays that equal relevant tooth passing period will be different.Consequently, in regard to the practical milling process, the governing dynamics equation should be modelled as delay differential equations with multiple delays instead.Based on [41,46,47], the governing equation of the milling process with multiple delays can be modelled by where  denotes the number of flutes, and the periodic coefficient matrix K  () is defined as Without loss of generality, we utilize milling with nonconstant pitch cutter to illustrate the mathematical model of milling with multiple delays.Based on [41,46,47], the tool is firstly divided into a finite number of disk elements along the axial direction.Then the resultant cutting forces in the  and  directions are acquired by numerically summing the force components acting on each individual element.In such conditions, ℎ  (), ℎ  (), ℎ  (), and ℎ  () can be finally obtained by where   and   are the tangential and normal cutting force coefficients, respectively.  (, ) denotes the angular position of the th tooth, given by where  denotes the helix angle,   represents the pitch angle between the th tooth and the ( − 1)th tooth,  is the radius of the cutter, and Ω represents the spindle speed.The screen function (  ()) is utilized to determine whether the tool is cutting the part, given by where  st and  ex represent the start and exit angles of tool.For downmilling,  st = arcos(2/-1), and  ex = ; for upmilling,  st = 0, and  ex = arcos(1-2/), where / is the radial immersion ratio.

Extended Adams-Moulton-Based Method
To numerically determine the milling stability, the governing equation ( 2) should be reexpressed in the state-space form.Specifically, define p() = M q() + Cq()/2 and y() = [q(), p()]  ; (2) can be rewritten as where By applying the state-space theory, the analytical response of ( 10) can be deduced as where y( 0 ) denotes initial state point and  0 represents the initial time instant.
Firstly, the spindle rotation period  is divided equally into  small time intervals, which makes  = ℎ.Obviously, these ( + 1) discretized points can be expressed as During the th time interval, that is,   <  <  +1 , ( 12) can be equivalently reexpressed as By substituting  =  +1 , y( +1 ) is defined as In order to simplify the derivation process, we will use some abbreviated expressions; that is, y  denotes y(  ), y − denotes y(  − ), and A  denotes A(  ).It is noted that (15) is derived from the analytical solution.Consequently, the key to solving (15) is to approximate the Duhamel term with high accuracy and numerical stability.Based on the two-step Adams-Moulton method, the state term y +1 can be approximated by When the time delay   is not equal to integer multiples of the step length ℎ, the delay term y −  should be interpolated by using the relevant boundary values.Define   = fix(  /ℎ); then   =   ℎ +   ,   =   −   ℎ, where the function fix() denotes the integer part of .For example, we can approximate y +1−  linearly by using the two boundary values, that is, y +1−  and y −  , resulting in Similarly, the delay terms y −  and y −1−  can be approximated linearly by By substituting ( 17)-( 19) into ( 16), one can read where In addition, y 2 can be approximated with the one-step Adams-Moulton formula Let  = 1 and substituting ( 17) and ( 18) into (22) yields where Obviously, y 1 and y +1− satisfy By combining ( 20), (23), and ( 25 ) , ) , Finally, the transition matrix Ψ with the EAMM can be written as According to Floquet theory, the stability of periodic systems depends on the spectral radius of the transition matrix.Consequently, the stability of milling operations can be obtained from the following criterion: where (Ψ) presents the spectral radius; that is, It should be noted that the construction of the Floquet transition matrix should be based on the period of the coefficient matrix rather than on that of the delay terms.For milling processes with single delay, the time period  is equal to the tooth passing period.However, it is equal to the spindle speed period for milling processes with multiple delays.Consequently, the spindle speed period is discretized.On the other hand, the ratio of time period to time delay (ROTPTD) is equal to one for milling processes with single delay, while it can be arbitrary for milling with multiple delays.In the proposed method, the delay terms are approximated by Lagrange interpolation polynomials with corresponding boundary values.Consequently, the proposed algorithm can be utilized for stability analysis of periodic delay systems with an arbitrary ROTPTD.

Validation and Comparison
In this section, the computation efficiency and accuracy of the proposed method are verified by a two-DOF milling operation with variable pitch cutter in [41].Note that Wan et al. [46] had extended the semidiscretization method in [21] to the milling stability analysis with multiple delays.To evaluate the effectiveness of the proposed method, we will make comparisons with the semidiscretization method (SDM) and the enhanced multistage homotopy perturbation method (EMHPM), in which the same program structure and the same model parameters are adopted.
Based on [41], the dynamics model of two-DOF milling operations can be modelled as In (30),   ,   ,   ,   ,   , and   represent the modal parameters of the system.The nonconstant pitch cutter has four flutes, a diameter of 19.05 mm, and a helix angle of 30 ∘ .There are three modes in the  direction and one mode in the  direction.In-depth analysis showed that the second mode in the  direction has dominant influence on the system stability.Consequently, the other two modes in the  direction are not taken into consideration.The modal parameters of the downmilling system are the same as [41]: the modal masses are   = 1.4986 kg and   = 1.199 kg, the natural angular frequencies are   = 563.6 × 2 rad/s and   = 516.21× 2 rad/s, and the damping ratios are   = 0.0558 and   = 0.025.The coefficient matrix K  () is the same as that given by ( 4).The time delay   represents the pitch period related to the pitch angle.The pitch angles of the cutter are 70 ∘ -110 ∘ -70 ∘ -110 ∘ , which results in four time delays; that is,  1 =  3 = 7/36, and  2 =  4 = 11/36.The part material is Al356 alloy, and the cutting force coefficients are   = 6.79 × 10 8 N/m 2 and   = 2.56 × 10 8 N/m 2 .
With the same matrix transformation presented in Section 2, (30) can be rewritten in the state-space form: where ) , To verify the feasibility of the proposed method, both large and low radial immersion milling conditions need to be investigated.First, the radial immersion ratio / is set as 1 and 0.6 to examine large radial immersions conditions.The stability lobe diagrams are constructed over a 200 × 120 sized grid.The time interval  is selected as 55 and 75, respectively.The cutting parameter combinations are set as follows: the spindle speed Ω ranges from 2 krpm to 12 krpm, and the depths of cut   ranges from 0 to 15 mm.The reference stability limits demoded by the red line are computed by the SDM with  = 600.Figures 1  and 2 show the stability lobe diagrams under large radial immersion conditions and the efficiency of these methods.The results demonstrate that the proposed method achieves a much higher computational efficiency than the SDM and the EMHPM.Indeed, the computation time of the proposed method can be reduced by 67-70% and 45-52%, compared with those of the SDM and the EMHPM, respectively.The stability lobe diagrams computed by the proposed method agree well with the exact stability lobe diagrams.In general, the accuracy of the stability lobe diagrams with the EAMM is higher than those with the other two methods based on the same computational parameters.The relative errors of these methods over the spindle speed range from 2000 rpm to 5000 rpm are presented in Figure 3.It demonstrates that the EAMM is of the highest accuracy.For instance, when compared with the SDM, the accuracy can be improved by up to 23% at the spindle speed Ω = 2000 rpm with the radial immersion ratio / = 1 and by up to 38% with the radial immersion ratio / = 0.6.
Meanwhile, we set the radial immersion ratio / set as 0.3 and 0.1 to examine low radial immersions milling.The time interval  is also selected as 55 and 75.The stability lobe diagrams are constructed over a 200 × 120 sized grid.The range of the spindle speed Ω remains the same as the previous case, while the depths of cut   range from 0 to 20 mm.The exact stability limits demoded by the red line are also calculated by the SDM with  = 600.Figures 4 and 5 show the stability lobe diagrams under low radial immersion conditions and the computational time of the SDM, the EMHPM, and the EAMM.It shows that the stability charts computed by the proposed method also agree well with the reference stability lobe diagrams for low radial immersions conditions.Moreover, the proposed method achieves a much higher computational efficiency than the other two methods.Compared with the SDM and the EMHPM, the computation time of the proposed method can be reduced by 67-69% and 45-53%, respectively.In addition, the results show that under the same computational parameters the accuracy of the EAMM is better than the SDM and the EMHPM.The relative errors of these three methods over the spindle speed range from 2000 rpm to 5500 rpm are presented in Figure 6.It shows that the accuracy of the EAMM is much better than the other two methods over this range of spindle speeds.For instance, when compared with the SDM, the accuracy can be improved by up to 22% at the spindle speed Ω = 2000 rpm and the radial immersion ratio / = 0.3 and by up to 19% at the spindle speed Ω = 3200 rpm and the radial immersion ratio / = 0.1.

Conclusion
In this work, an efficient and accurate semianalytical algorithm is proposed for the stability prediction of milling processes with multiple delays.Firstly, the milling dynamics model for regenerative chatter is rewritten in the state-space form.After the spindle rotation period is equally discretized, the delay terms are approximated by Lagrange interpolation polynomials, and the Adams-Moulton method is employed to construct the Floquet transition matrix.Finally, the stability of milling operations can be predicted by examining the spectral radius of the Floquet transition matrix.A two-DOF milling model with variable pitch tool has been adopted to demonstrate the proposed method.The numerical results demonstrate that under the same computational condition the proposed method achieves a higher computational efficiency than the SDM and the EMHPM.Compared with the SDM and the EMHPM, the computation time of the proposed method can be reduced by 67-70% and 45-53%, respectively.In general, the accuracy of the EAMM is higher than the SDM and the EMHPM.In addition, the accuracy of the stability lobe diagrams computed by the proposed method can be improved significantly over the spindle speed range from 2000 rpm to 5000 rpm.