Uncharted Stable Peninsula for Multivariable Milling Tools by High-Order Homotopy Perturbation Method

applied sciences Article Uncharted Stable Peninsula for Multivariable Milling Tools by High-Order Homotopy Perturbation Method Jose de la Luz Sosa 1, Daniel Olvera-Trejo 1,* , Gorka Urbikain 2,* , Oscar Martinez-Romero 1, Alex Elías-Zúñiga 1and Luis Norberto López de Lacalle 2 1Tecnológico de Monterrey, Escuela de Ingeniería y Ciencias, Av. Eugenio Garza Sada 2501, Monterrey, Nuevo León 64849, Mexico; [email protected] (J.d.l.L.S.); oscar[email protected] (O.M.-R.); [email protected] (A.E.-Z.) 2Department of Mechanical Engineering, University of the Basque Country, Alameda de Urquijo s/n, 48013 Bilbao, Spain; [email protected] *Correspondence: daniel.olvera.tr[email protected] (D.O.-T.); [email protected] (G.U.) Received: 9 October 2020; Accepted: 3 November 2020; Published: 6 November 2020   Abstract: In this work, a new method for solving a delay differential equation (DDE) with multiple delays is presented by using secondand third-order polynomials to approximate the delayed terms using the enhanced homotopy perturbation method (EMHPM). To study the proposed method performance in terms of convergency and computational cost in comparison with the first-order EMHPM, semi-discretization and full-discretization methods, a delay differential equation that model the cutting milling operation process was used. To further assess the accuracy of the proposed method, a milling process with a multivariable cutter is examined in order to find the stability boundaries. Then, theoretical predictions are computed from the corresponding DDE finding uncharted stable zones at high axial depths of cut. Time-domain simulations based on continuous wavelet transform (CWT) scalograms, power spectral density (PSD) charts and Poincar é maps (PM) were employed to validate the stability lobes found by using the third-order EMHPM for the multivariable tool. Keywords: chatter; multivariable tool; stable peninsula; homotopy perturbation method 1. Introduction There are many phenomena in different fields of science and engineering where the physical response of a variable involves not only the value at time t but also the effects that occur in an earlier state t−τ . Thus, delay systems appear in many engineering problems, such as in the shimmy effect (wheel vibration) [ 1 ], vehicle traffic models [ 2 ], feedback stabilization problems [ 3 ] and in the regenerative vibration of machine-tools better known as chatter [ 4 ]. In cases where the net force depends on the current values and some past values (history) such as position and speed, the system dynamic behavior can be modeled using a differential delay equation (DDE). It is well-known that during a milling process, unstable vibrations also known as self-excited vibration or chatter may occur. Chatter reduces the machining efficiency due to low material removal rate by reducing the workload and affects surface quality, shortens tool life and accelerates tool wear. Researchers are studying several ways to overcome this limitation. Kuljanic et al. [ 5 ] studied the incorporation of a chatter detection system based on multiple sensors to milling operations for industrial conditions, Zhuo et al. [ 6 ] used a method based on fractal dimension for the flank milling of a thin-walled blade, which can reflect the chatter severity level through the morphological change in signal. Paul and Morales [ 7 ], to mitigate chatter, presented an active controller based on the technique of discrete time sliding mode control (DSMC) blended with the type-2 fuzzy logic system. Moreover, Peng et al. [ 8 ] presented a method based on a dynamic cutting force simulation model and a Appl. Sci. 2020,10, 7869; doi:10.3390/app10217869 www.mdpi.com/journal/applsci Appl. Sci. 2020,10, 7869 2 of 22 machine learning approach based on statistical learning theory to predict and avoid the cutting chatter. In addition, to control and suppress chatter vibrations, the use of piezoelectric actuators embedded in the tool holder [ 9 ], electromagnetic actuators integrated into the spindle system [ 10 ] and tunable clamping table [ 11 ] has been analyzed. In the milling process, the use of variable pitch cutters has demonstrated to improve productivity [ 12 ]. Different from the uniform pitch cutter, when a variable pitch cutter is used the dynamics model of cutting vibration changes from DDEs with a single delay to DDEs with multiple delays [ 13 ]. A common technique offline to predict unstable vibrations is the so-called stability lobes of the DDE based on Floquet theory [ 14 ], in which a curve describes the limit of stable vibration under feasible range values of cutting parameters. The stability analysis of the milling process with multiple delays has been studied through different methods. Among all these methods, those with variable pitch tools play a critically important role [ 15 ]. Slavicek [ 16 ] was the one who first demonstrated the effectiveness of variable pitch cutters in suppressing vibrations in the milling process, he assumed a rectilinear tool motion for cutting teeth, and applied the theory of orthogonal stability to the irregular pitch of the tooth, by assuming an alternating step variation then, he obtained an expression of the stability limit as a function of the step angle variation. Budak [ 17 , 18 ] proposed an analytical method for nonconstant pitch milling cutters from a design point of view, showing for some applications how this variable effect helps to reduce self-excited vibrations, so he found that chatter stability can be improved significantly even at slow cutting speeds by properly designing the pitch angles. Altintas et al. [19] used the frequency domain method to analyze the milling stability of the variable pitch cutter and introduced a method to select the optimal pitch angles. Olgac and Sipahi proposed a mathematical approach, the cluster treatment of characteristic roots (CTCR), which optimizes the design of variable pitch cutters [ 20 ]. Jin et al. [ 21 ] presented an improved semi-discretization algorithm to predict the stability lobes for variable pitch cutters, which were verified and compared with previous works such as the Altintas analytical method (zero-order method) [ 19 ]. Comak and Budak [ 22 ] showed the optimal design of a tool for milling operations with variable geometry to widen the stability zones using the semi-discretization method, validating it experimentally. They also used a design methodology to determine the optimal pitch angle geometry for a given cutting condition, allowing increased stability. Zatarain et al. [ 23 ] extended the multifrequency solution proposed by Budak and Altintas [ 24 ] to include the helix effect, they pointed out that the variation of the helix angle plays an important role in stability graphs due to repetitive vibrations driven by impact (flip), they found that the flip lobes became closed curves that are separated by horizontal lines where the depth of cut is equal to a multiple of the helix pitch. A similar phenomenon was confirmed using the semi-discretization method (SDM) in [ 25 ], meanwhile, B.R. Patel et al. [ 26 ] considered the influence of the helix angle of the tool to obtain an analytical force model, they found that isolated islands of instability can occur in the milling processes, which are induced by the helix angle of the tool and lead to separate regions of period-doubling and quasi-period behavior. Sims et al. [ 27 ] by using an adapted and time-averaged version of the SDM analyzed both the influence on the variation of the helix angle and the pitch angle of the tool to improve the prediction of vibrations and estimate predictions of surface errors. They used the semi-discretization method, the time-averaged semi-discretization method and the temporal finite element method to predict vibration stability for variable helix and variable pitch milling tools. Turner et al. [ 28 ] modeled and compared stability for variable pitch and helix angle cutters, demonstrating that variable helix angle tools can have higher stability and productivity. Yusoffand Sims in [ 29 ] combined SDM with differential evolution to optimize variable helix milling tools to minimize vibration, their analysis predicted total vibration mitigation using the optimized variable helix milling tool at low radial immersion. Furthermore, Dombovari and Stepan [ 30 ] introduced a general mechanical model based on SDM to predict the linear stability of specialty cutters with optional continuous variation of the helix angle. Using an extended second-order SDM, Zhan et al. [ 15 ] predicted the stability lobe diagrams for tools with variable pitch angles. Meanwhile, Huang et al. [ 31 ] conducted a stability analysis for milling operations with variable pitch mills at variable speed, Appl. Sci. 2020,10, 7869 3 of 22 while Cai et al. [ 32 ] proposed an integrated process machine model based on the computer graphics method to simulate the milling process of a variable pitch cutter. On the other hand, Olvera and El í as-Zuñiga in [ 33 ] led to the development of the enhanced multistage homotopy perturbation method (EMHPM) to solve differential delay equations (DDEs) with constant and variable coefficients and then this EMHPM was applied to predict the stability of a multivariate milling tool in which they consider the helix angle and the pitch angle variation of the cutting edges [ 34 ]. Based on the Laplace formulation, Sims [ 35 ] studied the stability of milling operations with a variable helix angle. Using the multi-frequency solution, Otto et al. [ 36 ] derived a dynamic process model where the non-linear shear force and the runout effect are included for milling with non-uniform pitch and variable helix tools. Niu et al. [ 37 ] found that runout can significantly increase the stability limits regardless of spindle speed ranges, while Olvera et al. [ 38 ] in a study for a thin-walled workpiece demonstrated that by considering the effects of the runout, the helix angle and characterization dependent on the cutting speed, more precise stability boundaries are achieved. To demonstrate that one of the effective ways to suppress vibration in milling operations is to use tools with variable pitch and helix angle, Wang et al. [ 12 ] proposed an improved semi-discretization method based on Floquet 0 s theory. Since the delay between each cutting edge varies along with the axial depth of the tool in milling, they discretized the cutting tool in some axial layers to simplify the calculation. Iglesias et al. [ 39 ] presented a method to find the optimal angles between the inserts, and the stability diagrams were obtained through the iterative brute force (BF) method, which consists of an iterative maximization of stability through the semi-discretization method. They conclude that, if an optimal selection of the angle between the inserts is possible then, the material removal rate can be improved up to three times. Gou et al. [ 40 ] proposed an effective optimization method for the variable helical cutter introducing an index called “suppression factor” to measure stability quantitatively. Therefore, in the present work, the EMHPM developed in [ 33 ] and extended for analysis of multivariable tools in [ 34 ], is now expanded to solve the dynamics of the machining process in milling in which the approximation to the delay is performed with polynomials of degree two and three. In order to study the proposed method performance in terms of convergency and computational cost, a multivariable milling tool with a variable pitch cutter and helix angle is used to determine milling process in stability domains. This paper is summarized as follows. Section 2focuses on the development of secondand third-order EMHPM for stability analysis of DDE. Section 3studies the application of the secondand third-order EMHPM on the milling equation to demonstrate its improvement in the convergence rate. Section 4is focused on the use of the third-order EMHPM to compute the stability analysis in milling for multivariable tools, and theoretical predictions with time-domain simulations are performed. Finally, some conclusions are drawn. 2. Enhanced Multistage Homotopy Perturbation Method 2.1. Second-Order EMHPM Olvera et al. enhanced in [ 33 ] the multistage homotopy perturbation method (MHPM) proposed by Hashim [ 41 ]. The EMHPM considers the general case in which the nonlinear equation contains terms of the independent variable. This method is also useful to solve an n-dimensional DDE in the state-space form . x(t)=A(t)x+B(t)x(t−τ)(1) where A(t+τ)=A(t) , B(t+τ)=B(t) , x(t) , is the state vector, and τ is the time delay. Equation (1) can be written equivalently as: . xi(T)−Atxi(T)≈Btxiτ(T)(2) where xi(T) indicates the m -order solution for Equation (1) that satisfies the initial conditions xi(0)= xi−1 , At and Bt are the periodic matrix whose values vary with time t . In [ 42 ], Puma et al. applied Appl. Sci. 2020,10, 7869 4 of 22 the first-order EMHPM to estimate the delayed term xτ i(T) in Equation (2), in which the period [t0−τ,t0] was discretized in N equally spaced discrete state values, and the function that describes the delayed term xτ i(T) in the delayed interval [ti−N,ti−N+1] was approximated as a first-order polynomial representation. Defining xi≡xi(T)to simplify the notation, Equation (2) can be written as . xi(T)=Atxi(T) + Btxi−N+N−1 τ(xi−N+1−xi−N)T(3) Figure 1a shows the representation of the approximation of the delayed term with the first-order polynomial. In the second-order EMHPM, to approximate the function that describes the delayed term xτ i(T) in Equation (2), the Lagrange equation is used, making use of the discrete values xi−N,xi−N+1,xi−N+2as follows: fn(x)= n X i=0 Li(x)f(xi),Li(x)= n Y i=0,i,k x−xi xk−xi (4) to achieve a second-degree polynomial approximation, we have from Equation (4) that P2(x)=(x−∆t)(x−2∆t) (0−∆t)(0−2∆t)f(xi−N)+(x−0)(x−2∆t) (∆t−0)(∆t−2∆t)f(xi−N+1)+(x−0)(x−∆t) (2∆t−0)(2∆t−∆t)f(xi−N+2)(5) Appl. Sci. 2020, 10, x 4 of 23 () () () ititi TT T τ −≈xAx Bx  (2) where 𝐱(𝑇) indicates the 𝑚-order solution for Equation (1) that satisfies the initial conditions 𝐱(0) = 𝐱, 𝐀 and 𝐁 are the periodic matrix whose values vary with time 𝑡. In [42], Puma et al. applied the first-order EMHPM to estimate the delayed term 𝐱 (𝑇) in Equation (2), in which the period 󰇟𝑡−𝜏,𝑡 󰇠 was discretized in 𝑁 equally spaced discrete state values, and the function that describes the delayed term 𝐱 (𝑇) in the delayed interval 󰇟𝑡,𝑡󰇠 was approximated as a firstorder polynomial representation. Defining 𝐱≡𝐱 (𝑇) to simplify the notation, Equation (2) can be written as () () 1 1 ()xAxBx xx τ −−+− −  =++ −   i t i t iN iN iN N TT T (3) Figure 1a shows the representation of the approximation of the delayed term with the first-order polynomial. In the second-order EMHPM, to approximate the function that describes the delayed term 𝐱 (𝑇) in Equation (2), the Lagrange equation is used, making use of the discrete values 𝐱,𝐱,𝐱 as follows: () () ( ) () 00, , n n i nii iiik ki x x fx Lixfx Lx x x ==≠ − == − ∏ (4) to achieve a second-degree polynomial approximation, we have from Equation (4) that () ()( ) ()( ) () ()( ) ()( ) () ()( ) ()( ) () 212 202 0 002 02 202 xx x −−+ −+ −Δ − Δ − − Δ − −Δ =+ + −Δ − Δ Δ− Δ− Δ Δ− Δ−Δ iN iN iN xtx t x x t x xt Px f f f tt ttt ttt (5) Substituting 𝑥=𝑇 and ∆𝑡 = (𝑁 − 1)/𝜏, we obtain the function that describes the delayed interval as: () () 22 11212 13 1 1 22 22 2 iN iN iN iN iN iN iN iN NNT TT ττ −+ − − −+ −+ − −+ −+ −−   ≈+ − + − + − +     xx xxx xxx (6) When the delay is approximated by a second-degree polynomial it is called second-order EMHPM and should not be confused with the order of solution 𝑚 and which is determined by the last deformation taking into account the approximated solution. Notice that a polynomial of the second-degree requires three points. Likewise four points in the case of a third-degree polynomial as shown in Figure 1. (a) (b) Figure 1. (a) Scheme for the approximation of the delayed term by a first-order (solid black line), second-order (dashed blue line) and third-order (dotted red line) polynomial: (b) zoom in the time interval [𝑡,𝑡]. x t xi-N ti-N+1 ti-Nti-1 ti = (N-1) t xi-N+3 xi-N+2 xi-N+1 xi(t) t xi+1(t)xi+2(t) Figure 1. ( a ) Scheme for the approximation of the delayed term by a first-order (solid black line), second-order (dashed blue line) and third-order (dotted red line) polynomial: ( b ) zoom in the time interval [ti−N+1,ti−N+3]. Substituting x=T and ∆t=(N−1)/τ , we obtain the function that describes the delayed interval as: xi−N+1(T)≈xi−N+N−1 τT−3 2xi−N+2xi−N+1−1 2xi−N+2+N−1 τ2T2 2(xi−N−2xi−N+1+xi−N+2)(6) When the delay is approximated by a second-degree polynomial it is called second-order EMHPM and should not be confused with the order of solution m and which is determined by the last deformation taking into account the approximated solution. Notice that a polynomial of the second-degree requires three points. Likewise four points in the case of a third-degree polynomial as shown in Figure 1. The procedure to calculate the second-order EMHPM solution is based on the EMHPM procedure described in [33]. The solution for second-order EMHPM is recursively expressed of Xik(T)as Xik =Xa ik +Xb ik +Xc ik,k=1, 2, 3 . . . . (7) where Xa i0=xi−1, Xb i0=Xc i0=0 (8) Appl. Sci. 2020,10, 7869 5 of 22 and Xa ik =T kAtXa i(k−1)+g(k)Btxi−N Xb ik =T k+1AtXb i(k−1)+g(k)N−1 τBtT−3 2xi−N+2xi−N+1−1 2xi−N+2 Xc ik =T k+2AtXc i(k−1)+g(k)N−1 τ2BtT2 2(xi−N−2xi−N+1+xi−N+2) (9) So, the solution of Equation (1) is obtained by adding each of the approximations Xik of Equation (7). xi(T)≈ m X k=0 Xik(T)(10) 2.2. Third-Order EMHPM Solution For the polynomial representation of the third-degree, the function that describes the delayed term xτ i(T) is approximated by a polynomial of order three, then Equation (4) of the Lagrange interpolator is used accordingly. In this case, it is necessary to employ the xi−N , xi−N+1 , xi−N+2 , xi−N+3 discrete values. Following the same procedure described in Section 2.1, the function that describes the delayed interval is given as: xτ i(T)=xi−N+1(T)≈xi−N+N−1 τT−11 6xi−N+3xi−N+1−3 2xi−N+2+1 3xi−N+3+ N−1 τ2T2 2(2xi−N−5xi−N+1+4xi−N+2−xi−N+3)+N−1 τ3T3 6(−xi−N+3xi−N+1−3xi−N+2+xi−N+3) (11) Following the EMHPM procedure, the recursive solution of Equation (1) Xik(T)is expressed as Xik =Xa ik +Xb ik +Xc ik +Xd ik,k=1, 2, 3 . . . . (12) where Xa i0=xi−1,Xb i0=Xc i0=Xd i0=0 (13) and Xa ik =T kAtXa i(k−1)+g(k)Btxi−N Xb ik =T k+1AtXb i(k−1)+g(k)N−1 τBtT−11 6xi−N+3xi−N+1−3 2xi−N+2+1 3xi−N+3 Xc ik =T k+2AtXc i(k−1)+g(k)N−1 τ2BtT2 2(2xi−N−5xi−N+1+4xi−N+2−xi−N+3)) Xd ik =T k+3AtXd i(k−1)+g(k)N−1 τ3BtT3 6(−xi−N+3xi−N+1−3xi−N+2+xi−N+3)) (14) The approximate solution of Equation (1) can be obtained by substituting Equation (12) into Equation (10) adding each of the approximations Xik. 2.3. Stability Analysis To calculate the stability of the differential Equation (1) using the second-order EMHPM, the solution of Equation (10) for second-order EMHPM must be rewritten by grouping each of the discrete values xi,xi−N+2,xi−N+1,xi−N, resulting in xi(T)≈Pi(T)xi−1+Qi 0(T)xi−N+2+Qi(T)xi−N+1+Ri(T)xi−N(15) Appl. Sci. 2020,10, 7869 6 of 22 where Pi(T) = m P k=0 1 k!Ak tTk, Q0 i(T) = m P k=11 (k+2)!N−1 τ2Ak−1 tBtTk+2−1 2(k+1)!N−1 τAk−1 tBtTk+1 Qi(T) = m P k=1 1 (k+1)!N−1 τAk−1 tBtTk+1−2Q0 i Ri(T) = m P k=1 1 k!Ak−1 tBtTk−Q0 i−Qi (16) Similarly, to compute the stability lobes for the third-order EMHPM, the solution of the differential Equation (1) for third-order EMHPM is rewritten as xi(T)≈Pi(T)xi−1+Q”i(T)xi−N+3+Q0 i(T)xi−N+2+Qi(T)xi−N+1+Ri(T)xi−N(17) where Pi(T) = m P k=0 1 k!Ak tTk, Q”i(T) = m P k=11 (k+1)!N−1 τAk−1 tBtTk+11 3−1 (k+2)!N−1 τ2Ak−1 tBtTk+2+1 (k+3)!N−1 τ3Ak−1 tBtTk+3 Q0 i(T) = m P k=1  1 (k+1)!N−1 τAk−1 tBtTk+1+1 (k+2)!N−1 τ2Ak−1 tBtTk+2−7 2 +1 (k+3)!N−1 τ3Ak−1 tBtTk+39 2 −15 2Q0 0 i Qi(T) = m P k=11 (k+1)!N−1 τAk−1 tBtTk+1−3Q0 0 i−2Q0 i Ri(T) = m P k=1 1 k!Ak−1 tBtTk−Q0 0 i−Q0 i−Qi (18) The approximate solution obtained from Equation (17) was used to define a discrete map following the procedure described in [43]: wi=Diwi−1(19) where wi−1 is a vector with dimension equal to the total number of states (displacement and velocity) for all Ndiscrete intervals: wi−1= [x(i−1),. x(i−1),x(i−2),. . . ,x(i−N)]T(20) Diis a coefficient matrix and for the third-order EMHPM it has the form: Di=  P0 0 0 · · · 0Q00 iQ0iQiRi I0 0 0 · · · 0 0 0 0 0 0I0 0 · · · 0 0 0 0 0 0 0 I0· · · 0 0 0 0 0 . . .. . .. . ..... . .. . .. . .. . .. . .. . . 0 0 0 0 ...0 0 0 0 0 0 0 0 0 · · · I0 0 0 0 0 0 0 0 · · · 0I0 0 0 0 0 0 0 · · · 0 0 I0 0 0 0 0 0 · · · 0 0 0 I0  (21) Appl. Sci. 2020,10, 7869 7 of 22 It is important to point out that in the case of the second-order EMHPM, the matrix Di is like the matrix of the third-order EMHPM without the matrix Q00 i. Then, the Floquet transition matrix Φ is calculated over the main period τ=(N−1)/∆t , coupling each of the discrete maps Di,i=1, 2, . . . ,(N−1), to obtain: Φ=DN−1DN−2. . . D2D1(22) Thus, the stability of Equation (1) is determined by calculating the eigenvalues of the transition matrix given by Equation (22). The eigenvalues of the transition matrix are actually the Floquet multipliers which are the exponents of each complex exponential functions that describe the motion of Equation (1). If the modulus of greatest magnitude is greater than or equal to one, it implies that the system will behave in an unstable way and the amplitude of the vibration will increase exponentially, otherwise it will have a stable behavior. 3. Numeric Solution of the Milling Equation 3.1. Dynamic Model to the Milling Equation To validate the proposed EMHPM methods, the numerical solution of the delay differential equation analyzed by Olvera et al., in [ 33 ] was calculated, which describes the dynamic model of the milling process in one degree of freedom (DOF): .. x(t)+2ζωn . x(t)+ω2 nx(t) = − aphxx(t) mm (x(t)−x(t−τ)) (23) where ζ is the modal damping ratio, ωn is the natural frequency of the workpiece, ap is the axial depth of cut, mm is the modal mass, τ represents the time delay corresponding to the hitting period between each tooth of the tool and hxx(t) is the specific cutting force in the x-direction due to flexibility in x-direction, which was calculated depending on the position of the tool hxx(t) = zn X iz=1 g(φiz(t))sinφiz(t)(Ktc cos φiz(t) + Knc sin φiz(t))(24) zn is the number of edges of the tool, Ktc and Knc are the average specific cut coefficients in the tangential and normal direction, respectively, and φiz(t)is the angular position of each left edge described by φiz =(2πn/60)t+2πiz/zn(25) where n is the spindle speed in revolution per minute (rpm). The function g(φiz(t)) is a window function, which has the value of one when the current edge iz is cutting material, otherwise it takes the value zero. In up-milling φst = 0 and φex =cos−1(1−2ad) , conversely in down-milling φst =cos−1(2ad−1) and φex =π , ad is the radial immersion ratio of the cut and φst and φex are the angular positions where each edge enters and leaves the workpiece. The secondand third-order EMHPM is applied to obtain the solution of Equation (23) and it is compared with the solution given by the first-order EMHPM [ 33 ]; for a regular tool the matrix A and B are represented as: At= 0 1 −ω2 n−aphxx(t) mm −2ζωn,Bt= 0 0 aphxx(t) mm0(26) At and Bt correspond to the periodic matrix evaluated at time t . For demonstration purposes, time-domain simulations were computed for a full-immersion down-milling operation. We used the parameters employed by Insperger et al., in [ 43 ] where the stability lobes were also calculated. Appl. Sci. 2020,10, 7869 8 of 22 The modal parameters fn =922 Hz, ωn= 5793 rad/s, ζ= 0.011 and mm= 0.03993 kg corresponds to a single degree of freedom. The tangential and normal cutting coefficients are Ktc = 6 × 10 8 N/ m2 and Knc = 2 × 10 8 N/ m2 respectively for an end-mill with zn= 2. The time-domain solution was computed using the EMHPM considering N= 76 discrete intervals and m= 7. Two sets of cutting conditions were chosen for a fixed spindle speed value of n= 12,000 rpm where the axial depth of cut of ap=1.5 mm corresponds to a stable cutting operation while that for an unstable operation ap=3 mm was chosen. In Figure 2we plot the secondand third-order EMHPM solutions and compare it with the first-order EMHPM and the dde23 routine in Matlab, which is used to integrate DDE. Appl. Sci. 2020, 10, x 10 of 25 (a) (b) Figure 2. Numerical comparison of the enhanced homotopy perturbation method (EMHPM) solutions of the milling equation, Equation Error! Reference source not found., with the dde23 MATLAB routine. (a) Stable milling operation with 𝑎=1.5 mm, 𝑎=1 and 𝑛 = 12000 rpm and (b) unstable milling operation with 𝑎=3 mm, 𝑎=1 and 𝑛 = 12000 rpm. 3.2. Numerical Comparison between Methods In order to observe the rate of convergence of the first-, secondand third-order EMHPM, we chose the stable case with cutting conditions 𝑎=1.5 mm, 𝑎=1 and 𝑛 = 12000 rpm presented in Figure 3a, and the unstable case with cutting conditions 𝑎=3 mm, 𝑎=1 and 𝑛 = 12000 rpm showed in Figure 3b. The rate of convergence was analyzed by computing the absolute error between the solution with N discrete intervals and a converged solution. All methods were compared against itself using the solution provided with N = 200 discrete intervals, which are considered the converged solution. In Figure 3a it is observed that the convergence is better for the secondand third-order than the first-order, however, the difference of convergence between secondand third-order with the parameters used was negligible. On the other hand, Figure 3b shows that for few discrete intervals the third-order EMHPM had the fastest convergence in comparison with the secondand the firstorder EMHPM. However, the second-order and third-order curves behaved very similarly after N = 50 discrete intervals. It is important to mention that for a typical stability solution in the ranges of spindle speed 5000–10000 rpm, N = 40 discrete intervals will be enough to have accurate predictions. 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 Time ( s ) -1 -0.5 0 0.5 1 Displacement (m) 10-3 Figure 2. Numerical comparison of the enhanced homotopy perturbation method (EMHPM) solutions of the milling equation, Equation (23), with the dde23 MATLAB routine. ( a ) Stable milling operation with ap= 1.5 mm, ad= 1 and n= 12000 rpm and ( b ) unstable milling operation with ap= 3 mm, ad=1 and n=12000 rpm. 3.2. Numerical Comparison between Methods In order to observe the rate of convergence of the first-, secondand third-order EMHPM, we chose the stable case with cutting conditions ap= 1.5 mm, ad= 1 and n= 12,000 rpm presented in Figure 3a, and the unstable case with cutting conditions ap= 3 mm, ad= 1 and n= 12,000 rpm showed in Figure 3b. The rate of convergence was analyzed by computing the absolute error between the solution with Ndiscrete intervals and a converged solution. All methods were compared against itself using the solution provided with N=200 discrete intervals, which are considered the converged solution. In Figure 3a it is observed that the convergence is better for the secondand third-order than the first-order, however, the difference of convergence between secondand third-order with the parameters used was negligible. On the other hand, Figure 3b shows that for few discrete intervals the third-order EMHPM had the fastest convergence in comparison with the secondand the first-order EMHPM. However, the second-order and third-order curves behaved very similarly after N=50 discrete intervals. It is important to mention that for a typical stability solution in the ranges of spindle speed 5000–10,000 rpm, N=40 discrete intervals will be enough to have accurate predictions. Appl. Sci. 2020,10, 7869 9 of 22 Appl. Sci. 2020, 10, x 11 of 25 (a) (b) Figure 3. Convergence rate of absolute error between first-, secondand third-order EMHPM for down-milling operation. Cutting parameters for (a) 𝑎= 1.5 mm,𝑎=1 and 𝑛 = 12000 rpm and (b) 𝑎= 3 mm,𝑎=1 and 𝑛 = 12000 rpm. Since the rate of convergence was proved for time-domain simulations, we next explored the convergence of the methods applied to the stability analysis. The stability lobes computed with the secondand third-order EMHPM for regular milling tools were compared with its predecessor for radial immersion value of 𝑎 =1 and the other parameters indicated above as it was used in [44]. Figure 4 shows the stability diagrams for spindle speed in the range 2000–3000 rev/min where the precision of the method was compromised due to the higher value of the time delay. While the shaded gray area represents the stability lobes computed with N = 200 discrete intervals in all subfigures, in each subfigure solid black lines draw the stability frontier for a specific discrete interval and using the first-, secondor third-order EMHPM. In Figure 4 the first, second and third column represents the solution for the first-, the secondand the third-order EMHPM respectively, while the first and the second row was for N = 60 and N = 100 discrete intervals, respectively. It is observed that the error achieved in the third-order EMHPM was less than those attained for the first-order and second-order EMHPM solutions. This confirms that the third-order EMHPM had the highest rate of convergence. Absolute Error Figure 3. Convergence rate of absolute error between first-, secondand third-order EMHPM for down-milling operation. Cutting parameters for ( a ) ap= 1.5 mm , ad= 1 and n=12,000 rpm and (b)ap=3 mm, ad=1 and n=12,000 rpm. Since the rate of convergence was proved for time-domain simulations, we next explored the convergence of the methods applied to the stability analysis. The stability lobes computed with the secondand third-order EMHPM for regular milling tools were compared with its predecessor for radial immersion value of ad =1 and the other parameters indicated above as it was used in [ 44 ]. Figure 4shows the stability diagrams for spindle speed in the range 2000–3000 rev/min where the precision of the method was compromised due to the higher value of the time delay. While the shaded gray area represents the stability lobes computed with N=200 discrete intervals in all subfigures, in each subfigure solid black lines draw the stability frontier for a specific discrete interval and using the first-, secondor third-order EMHPM. In Figure 4the first, second and third column represents the solution for the first-, the secondand the third-order EMHPM respectively, while the first and the second row was for N=60 and N=100 discrete intervals, respectively. It is observed that the error achieved in the third-order EMHPM was less than those attained for the first-order and second-order EMHPM solutions. This confirms that the third-order EMHPM had the highest rate of convergence. Appl. Sci. 2020,10, 7869 16 of 22 tooth 0.05, 0.10, 0.015 and 0.20, so the resulting cutting coefficients Ktc for the tooth 1, 2, 3 and 4 were 1215 × 10 6 , 1369 × 10 6 , 897 × 10 6 and 1799 × 10 6 N/ m2 respectively, while that the coefficients Knc for the tooth 1, 2, 3 and 4 resulted 272 ×106, 520 ×106, 801 ×106and 859 ×106N/m2respectively. Table 2. Main geometric parameters of multivariable tool. Appl. Sci. 2020, 10, x 17 of 25 4.2. Experimental Characterization of One Degree of Freedom Milling Equation and Cutting Force Model 4.2.1. Experimental Modal Analysis An experimental workpiece was assembled with a 7075T6 aluminum block of 101 mm × 172 mm supported by two thin plates (walls) with a thickness of 4.5 mm. This assembly mimics a DOF as described in Equation Error! Reference source not found.. The workpiece assembly was rigidly fixed to the workbench of a Makino F3 machining center. For modal analysis, tap testing was performed using a 352C68 PCB Piezotronics accelerometer and an impact hammer model 9722A500. The signals were acquired with a Polytec VIB-E-220 data acquisition card and processed with VibSoft signal analyzer software as shown in Figure 5a. Using the CutPro 8 software, the modal parameters were fitted resulting the values 𝜁 = 0.068, 𝑚=3.8 kg, 𝑓  = 132 Hz and 𝜔= 829 rad/s. 4.2.2. Experimental Determination of Cutting Coefficients The force model in Equation Error! Reference source not found. was used to predict the cutting force magnitude for a given depth of cut. It is based on a mechanistic approach that assumes a relationship between forces and the uncut chip thickness by means of the cutting coefficients. The cutting force model was established by introducing cutting (shearing) and edge coefficients for the tangential and normal directions of the milling tool. The characterization procedure assumed the linear relationship between the averaged experimental cutting forces 𝐹  and the feed rate fz in xand ydirections. This relationship is established as follows: =+  z ce FfFF (42) Here, 𝐹   and 𝐹   are the cutting shear and edge components, respectively. The experimental forces at each feed rate are measured, and the cutting-edge components 𝐹   and 𝐹   were evaluated 4, 4 yc x c tc nc np np F F KK za za ==−  (43) A multivariable cutter provided by a local toolmaker was characterized by using Equation Error! Reference source not found. and the experimental setup shown in Figure 5b. Table 2 summarizes the main geometric characteristics of the multivariable tool. A total of five cuttings were performed for full radial immersion in aluminum 7075T6 during dry machining. The forces were recorded by using a dynamometer 9257B Kistler and the spindle speed was set at 3000 rpm based on the dynamometer’s natural frequency to avoid the amplification of milling forces. The force signals were acquired using a VibSoft-20 acquisition card at a sample rate of 48 kHz and processed in a custom-made MATLAB app to remove drift and noise. Cutting forces data were collected for the axial depth of cut of 2 mm and four values of feed per tooth 0.05, 0.10, 0.015 and 0.20, so the resulting cutting coefficients 𝐾 for the tooth 1, 2, 3 and 4 were 1215 × 10, 1369 × 10, 897 × 10 and 1799 × 10 N/m respectively, while that the coefficients 𝐾 for the tooth 1, 2, 3 and 4 resulted 272 × 10, 520 × 10, 801 × 10 and 859 × 10 N/m respectively. Table 2. Main geometric parameters of multivariable tool. Diameter 12.7 mm Cutting length 25 mm Coating type Uncoated Number of teeth 4 Helix angles 39°, 37°, 39°, 41° Pitch angles 80°, 100°, 70°, 110° Diameter 12.7 mm Cutting length 25 mm Coating type Uncoated Number of teeth 4 Helix angles 39◦, 37◦, 39◦, 41◦ Pitch angles 80◦, 100◦, 70◦, 110◦ 4.3. Stability Analysis of 1 DOF Milling with a Multivariable Tool The stability lobes computed for the multivariable tool using the third-order EMHPM with a mesh of 400 × 200 ( n×ap ) are shown in Figure 6together with stability lobes for a regular tool (angles of 90 ◦ and helix angles of 30 ◦ for all flutes). An approximation of order m= 7 was used with N=241 and ad= 1 mm. Notice from Figure 6that the stable zone obtained for the multivariable tool was significantly larger, meaning that the critical depth of cut was higher in most spindle speeds, which allowed having more global productivity. It is also observed in the range of spindle speed between 2000 and 3000 rpm, a stable peninsula formed with axial depth ranging from 11 to 20 mm or higher values of critical depth of cut. For instance, for the multivariable cutter at 2500 rpm, the critical depth of cut ap was 2.17 mm, however it became stable again as shown in Figure 6for the interval values between 11 and 20 mm. To validate this unexpected behavior, we performed several time-domain simulations using the third-order EMHPM solution described by Equation (37). Appl. Sci. 2020, 10, x 18 of 25 (a) (b) Figure 5. Scheme of the experimental setup for (a) the modal analysis and (b) cutting forces characterization. 4.3. Stability Analysis of 1 DOF Milling with a Multivariable Tool The stability lobes computed for the multivariable tool using the third-order EMHPM with a mesh of 400 ×200 (𝑛 × 𝑎𝑝) are shown in Figure 6 together with stability lobes for a regular tool (angles of 90° and helix angles of 30° for all flutes). An approximation of order 𝑚 = 7 was used with N=241 and 𝑎𝑑= 1 mm. Notice from Figure 6 that the stable zone obtained for the multivariable tool was significantly larger, meaning that the critical depth of cut was higher in most spindle speeds, which allowed having more global productivity. It is also observed in the range of spindle speed between 2000 and 3000 rpm, a stable peninsula formed with axial depth ranging from 11 to 20 mm or higher values of critical depth of cut. For instance, for the multivariable cutter at 2500 rpm, the critical depth of cut 𝑎𝑝 was 2.17 mm, however it became stable again as shown in Figure 6 for the interval values between 11 and 20 mm. To validate this unexpected behavior, we performed several time-domain simulations using the third-order EMHPM solution described by Equation Error! Reference source not found.. (a) (b) Figure 6. (a) Comparison of stability lobes for regular (black solid line) and multivariable (red solid line) cutters by using the third-order EMHPM and (b) zoom in on chosen cutting conditions for timedomain simulations. The selected points are marked as follows: unstable (cross mark), stable (circle mark) and transition (plus mark) cutting conditions. Furthermore, the simulated vibrations for the chosen cutting conditions were analyzed using the continuous wavelet transform (CWT), the power spectral density (PSD) and Poincaré maps (PM). The CWT is a time-frequency representation of a signal that offers the capability to observe how frequencies evolve in time. The scalograms display the absolute value of CWT of the simulated vibration and therefore, they were used to detect chatter phenomena that appeared when milling with a multivariable tool. The PSD is based on the Fourier transform that provides the transformation from the time-domain to the frequency-domain. Additionally, PSD is defined as the squared value of the signal and describes the power of a signal or time series distributed over different frequencies [46]. Moreover, a PM represents points in phase space, which are sampled every spindle rotation [47]. The frequencies 𝑓 of the CWT and PM were normalized 𝑓 𝑛= 𝑓 𝑓ℎ ⁄ according to the spindle Figure 6. ( a ) Comparison of stability lobes for regular (black solid line) and multivariable (red solid line) cutters by using the third-order EMHPM and ( b ) zoom in on chosen cutting conditions for time-domain simulations. The selected points are marked as follows: unstable (cross mark), stable (circle mark) and transition (plus mark) cutting conditions. Furthermore, the simulated vibrations for the chosen cutting conditions were analyzed using the continuous wavelet transform (CWT), the power spectral density (PSD) and Poincar é maps (PM). The CWT is a time-frequency representation of a signal that offers the capability to observe how frequencies evolve in time. The scalograms display the absolute value of CWT of the simulated vibration and therefore, they were used to detect chatter phenomena that appeared when milling with a multivariable tool. The PSD is based on the Fourier transform that provides the transformation from the time-domain to the frequency-domain. Additionally, PSD is defined as the squared value of the signal and describes the power of a signal or time series distributed over different frequencies [ 46 ]. Moreover, a PM represents points in phase space, which are sampled every spindle rotation [ 47 ]. Appl. Sci. 2020,10, 7869 17 of 22 The frequencies f of the CWT and PM were normalized fn=f/fh according to the spindle frequency fh . When milling with a regular milling tool the excitation frequency fe is equal to zn times frequencies of the spindle speed fh but in a multivariable tool, there are several excitation frequencies since the angular spacing between teeth change as a function of the axial depth of cut. Figure 7illustrates the CWT, PSD and PM for simulated vibrations using the multivariable tool with different axial depths denoted as cutting conditions A, B and C for the axial depths of cut of 1.0, 1.7 and 1.7 mm respectively. Figure 7a–c refers to the vibrations of the cutting conditions A marked in Figure 6, using a regular tool. The scalogram in Figure 7a identifies point A as a stable cutting since normalized cutting frequencies present a dominant value of fn =3.2, which corresponds to the natural frequency fm= 132 Hz. This is also confirmed by the PSD analysis shown in Figure 7b. The PM illustrated in Figure 7c shows a vibration that decreased with time and sampled data concentrated in the center confirmed a typical stable case. When the axial depth of cut was increased to 1.7 mm, the stability diagram predicted unstable cutting conditions according to the stability lobes for the regular tool. This case is denoted with cutting conditions B and the corresponding scalogram (shown in Figure 7d) illustrated how the intensity of the dominant frequency increased with time even when the excitation frequency was the same as the case in A. Appl. Sci. 2020, 10, x 19 of 25 frequency 𝑓ℎ. When milling with a regular milling tool the excitation frequency 𝑓 𝑒 is equal to 𝑧𝑛 times frequencies of the spindle speed 𝑓ℎ but in a multivariable tool, there are several excitation frequencies since the angular spacing between teeth change as a function of the axial depth of cut. Figure 7 illustrates the CWT, PSD and PM for simulated vibrations using the multivariable tool with different axial depths denoted as cutting conditions A, B and C for the axial depths of cut of 1.0, 1.7 and 1.7 mm respectively. Figure 7a–c refers to the vibrations of the cutting conditions A marked in Figure 6, using a regular tool. The scalogram in Figure 7a identifies point A as a stable cutting since normalized cutting frequencies present a dominant value of 𝑓 𝑛 = 3.2, which corresponds to the natural frequency 𝑓 𝑚=132 Hz. This is also confirmed by the PSD analysis shown in Figure 7b. The PM illustrated in Figure 7c shows a vibration that decreased with time and sampled data concentrated in the center confirmed a typical stable case. When the axial depth of cut was increased to 1.7 mm, the stability diagram predicted unstable cutting conditions according to the stability lobes for the regular tool. This case is denoted with cutting conditions B and the corresponding scalogram (shown in Figure 7d) illustrated how the intensity of the dominant frequency increased with time even when the excitation frequency was the same as the case in A. (a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 7. Analysis of cutting conditions A, B and C. Continuous wavelet transform (CWT) scalograms: (a,d,g); power spectral density (PSD): (b,e,h) and Poincaré maps (PM): (c,f,i) corresponds to the cutting conditions A, B and C respectively. Figure 7. Analysis of cutting conditions A, B and C. Continuous wavelet transform (CWT) scalograms: ( a , d , g ); power spectral density (PSD): ( b , e , h ) and Poincar é maps (PM): ( c , f , i ) corresponds to the cutting conditions A, B and C respectively. Appl. Sci. 2020,10, 7869 18 of 22 The PM diagram shown in Figure 7f exhibited a vibration far from zero. In fact, the PM diagram shows that the vibration amplitude grows exponentially because our equation of motion did not consider nonlinear effects such as those that appeared when the tool lost contact with the workpiece. Both cutting conditions A and B agreed with the stability boundaries in Figure 6. Now, the cutting conditions B were used but with a multivariable tool, which was referred to as cutting conditions C. The CWT plotted in Figure 7g described completely different results since there were no single dominant frequencies in comparison with cutting conditions A, but appeared several frequencies around fn =3.2 and close to fn =1 that reduced in intensity with time, suggesting a stable cutting. Figure 7i illustrates how the vibration amplitude approached to zero when using a multivariable tool in contrast to the PM obtained for the regular tool and exhibited in Figure 7f. This can be explained by observing that there were several excitation frequencies due to the irregular pitch and helix angles that break a single excitation frequency avoiding regenerative chatter phenomena. Figure 8illustrates the CWT, PSD and PM for simulated vibrations using the multivariable tool with different axial depths denoted as cutting conditions D, E, F and G for the axial depths of cut of 2.3, 3.0, 8.55 and 18 mm respectively. Notice that a stable case C was already validated when the axial depth was 1.7 mm in Figure 7g–i that corresponded to cutting conditions under the stability boundaries shown in Figure 6. For case D, a transient cutting condition was chosen very close to the critical axial depth of the cut. It is interesting to point out that transition cutting conditions in the CWT scalogram shown in Figure 8a not only shows frequencies with higher intensity in comparison with the stable case B, but also presents shifted frequencies that varied in intensity every single revolution. This shifting suggests a marginally stable cutting condition that was confirmed by the PM illustrated in Figure 8c, where circular trajectories were described close to the center point. Appl. Sci. 2020, 10, x 20 of 25 The PM diagram shown in Figure 7f exhibited a vibration far from zero. In fact, the PM diagram shows that the vibration amplitude grows exponentially because our equation of motion did not consider nonlinear effects such as those that appeared when the tool lost contact with the workpiece. Both cutting conditions A and B agreed with the stability boundaries in Figure 6. Now, the cutting conditions B were used but with a multivariable tool, which was referred to as cutting conditions C. The CWT plotted in Figure 7g described completely different results since there were no single dominant frequencies in comparison with cutting conditions A, but appeared several frequencies around 𝑓 𝑛 = 3.2 and close to 𝑓 𝑛 = 1 that reduced in intensity with time, suggesting a stable cutting. Figure 7i illustrates how the vibration amplitude approached to zero when using a multivariable tool in contrast to the PM obtained for the regular tool and exhibited in Figure 7f. This can be explained by observing that there were several excitation frequencies due to the irregular pitch and helix angles that break a single excitation frequency avoiding regenerative chatter phenomena. Figure 8 illustrates the CWT, PSD and PM for simulated vibrations using the multivariable tool with different axial depths denoted as cutting conditions D, E, F and G for the axial depths of cut of 2.3, 3.0, 8.55 and 18 mm respectively. Notice that a stable case C was already validated when the axial depth was 1.7 mm in Figure 7g–i that corresponded to cutting conditions under the stability boundaries shown in Figure 6. For case D, a transient cutting condition was chosen very close to the critical axial depth of the cut. It is interesting to point out that transition cutting conditions in the CWT scalogram shown in Figure 8a not only shows frequencies with higher intensity in comparison with the stable case B, but also presents shifted frequencies that varied in intensity every single revolution. This shifting suggests a marginally stable cutting condition that was confirmed by the PM illustrated in Figure 8c, where circular trajectories were described close to the center point. (a) (b) (c) (d) (e) (f) Figure 8. Cont. Appl. Sci. 2020,10, 7869 19 of 22 Appl. Sci. 2020, 10, x 21 of 25 (g) (h) (i) (j) (k) (l) Figure 8. Analysis of cutting conditions D, E, F and G. CWT scalogram:s (a,d,g,j); PSD: (b,e,h,k) and PM: (c,f,i,l) corresponds to the cutting conditions D, E, F and G respectively. Unstable vibrations that appeared for case E were because of the intensity of frequencies increased exponentially with time, see Figure 8d. Notice that other frequencies arose with time close to the values of 𝑓 𝑛 = 0.5 and 𝑓 𝑛 = 1.5. These frenquencies also occurred for cutting conditions D, which is an indication of the appearance of chatter phenomena. In contrast to Figure 7i for a stable case, Figure 8f exhibited few trajectories because the vibration amplitude was out of the range selected (±1 mm). The qualitative and quantitative dynamic behaviour due to cutting conditions F, and illustrated in Figure 8g-i, were classified as transition cutting behaviour. Here, a more severe shifting in frequencies was observed in the scalogram (Figure 8g). From Figure 8g, it is seen that drastic shifting occurred in the time domain in the range of normalized frequencies from 3.5 to 6. It was also evident in the PM showed in Figure 8i, that the amplitude of vibration remained below 1 mm during several revolutions of the tool but the amplitude of vibration never aproached to the center point, in contrast to the stable cutting condition C shown in Figure 7i in which the oscillation aplitudes aproached the center. An interesting dynamic behaviour was observed in the milling cutting process when the cutting conditions were selected in the middle of the stable peninsula, above unstable cutting conditions such as E cutting conditions. The axial depth of the cut was increased from the unstable axial depth of cut of 3–18 mm, 6 times higher of the stable cutting condition C and 2 times higher than the unstable condition E. Since the vibration quickly decreased in a few revolutions no dominant frequencies appeared in the CWT and PSD failed to clearly identify a dominant frequency since the vibration amplitude decreased to zero after few revolutions, as confirmed by the PM shown in Figure 8l. Figure 9 shows the normalized excitation frequencies that the multivariable tool produced for a fixed spindle speed of 2500 rpm. The total number of disks of 50 μm of thickness was grouped in sets of each millimeter in the axial direction. The waterfall plot in Figure 9 explains that a stable peninsula was formed above 11 mm because the workpiece was excited with several frequencies simultaneously. For instance, for a milling operation with the axial depth of cut of 1 mm (stable cutting), 80 discrete disks were cut with four normalized excitation frequencies values (3.3, 3.6, 4.5 and 5.1). On the other hand, when milling at 18 mm (stable cutting), there were 14 normalized Appl. Sci. 2020, 10, x 21 of 25 (g) (h) (i) (j) (k) (l) Figure 8. Analysis of cutting conditions D, E, F and G. CWT scalogram:s (a,d,g,j); PSD: (b,e,h,k) and PM: (c,f,i,l) corresponds to the cutting conditions D, E, F and G respectively. Unstable vibrations that appeared for case E were because of the intensity of frequencies increased exponentially with time, see Figure 8d. Notice that other frequencies arose with time close to the values of 𝑓 𝑛 = 0.5 and 𝑓 𝑛 = 1.5. These frenquencies also occurred for cutting conditions D, which is an indication of the appearance of chatter phenomena. In contrast to Figure 7i for a stable case, Figure 8f exhibited few trajectories because the vibration amplitude was out of the range selected (±1 mm). The qualitative and quantitative dynamic behaviour due to cutting conditions F, and illustrated in Figure 8g-i, were classified as transition cutting behaviour. Here, a more severe shifting in frequencies was observed in the scalogram (Figure 8g). From Figure 8g, it is seen that drastic shifting occurred in the time domain in the range of normalized frequencies from 3.5 to 6. It was also evident in the PM showed in Figure 8i, that the amplitude of vibration remained below 1 mm during several revolutions of the tool but the amplitude of vibration never aproached to the center point, in contrast to the stable cutting condition C shown in Figure 7i in which the oscillation aplitudes aproached the center. An interesting dynamic behaviour was observed in the milling cutting process when the cutting conditions were selected in the middle of the stable peninsula, above unstable cutting conditions such as E cutting conditions. The axial depth of the cut was increased from the unstable axial depth of cut of 3–18 mm, 6 times higher of the stable cutting condition C and 2 times higher than the unstable condition E. Since the vibration quickly decreased in a few revolutions no dominant frequencies appeared in the CWT and PSD failed to clearly identify a dominant frequency since the vibration amplitude decreased to zero after few revolutions, as confirmed by the PM shown in Figure 8l. Figure 9 shows the normalized excitation frequencies that the multivariable tool produced for a fixed spindle speed of 2500 rpm. The total number of disks of 50 μm of thickness was grouped in sets of each millimeter in the axial direction. The waterfall plot in Figure 9 explains that a stable peninsula was formed above 11 mm because the workpiece was excited with several frequencies simultaneously. For instance, for a milling operation with the axial depth of cut of 1 mm (stable cutting), 80 discrete disks were cut with four normalized excitation frequencies values (3.3, 3.6, 4.5 and 5.1). On the other hand, when milling at 18 mm (stable cutting), there were 14 normalized Figure 8. Analysis of cutting conditions D, E, F and G. CWT scalogram:s ( a , d , g , j ); PSD: ( b , e , h , k ) and PM: (c,f,i,l) corresponds to the cutting conditions D, E, F and G respectively. Unstable vibrations that appeared for case E were because of the intensity of frequencies increased exponentially with time, see Figure 8d. Notice that other frequencies arose with time close to the values of fn =0.5 and fn =1.5. These frenquencies also occurred for cutting conditions D, which is an indication of the appearance of chatter phenomena. In contrast to Figure 7i for a stable case, Figure 8f exhibited few trajectories because the vibration amplitude was out of the range selected ( ± 1 mm). The qualitative and quantitative dynamic behaviour due to cutting conditions F, and illustrated in Figure 8g-i, were classified as transition cutting behaviour. Here, a more severe shifting in frequencies was observed in the scalogram (Figure 8g). From Figure 8g, it is seen that drastic shifting occurred in the time domain in the range of normalized frequencies from 3.5 to 6. It was also evident in the PM showed in Figure 8i, that the amplitude of vibration remained below 1 mm during several revolutions of the tool but the amplitude of vibration never aproached to the center point, in contrast to the stable cutting condition C shown in Figure 7i in which the oscillation aplitudes aproached the center. An interesting dynamic behaviour was observed in the milling cutting process when the cutting conditions were selected in the middle of the stable peninsula, above unstable cutting conditions such as E cutting conditions. The axial depth of the cut was increased from the unstable axial depth of cut of 3–18 mm, 6 times higher of the stable cutting condition C and 2 times higher than the unstable condition E. Since the vibration quickly decreased in a few revolutions no dominant frequencies appeared in the CWT and PSD failed to clearly identify a dominant frequency since the vibration amplitude decreased to zero after few revolutions, as confirmed by the PM shown in Figure 8l. Figure 9shows the normalized excitation frequencies that the multivariable tool produced for a fixed spindle speed of 2500 rpm. The total number of disks of 50 µ m of thickness was grouped in sets of each millimeter in the axial direction. The waterfall plot in Figure 9explains that a stable peninsula was formed above 11 mm because the workpiece was excited with several frequencies simultaneously. For instance, for a milling operation with the axial depth of cut of 1 mm (stable cutting), 80 discrete disks were cut with four normalized excitation frequencies values (3.3, 3.6, 4.5 and 5.1). On the other hand, when milling at 18 mm (stable cutting), there were 14 normalized excitation frequencies (3.30, Appl. Sci. 2020,10, 7869 20 of 22 3.35, 3.39, 3.44, 3.49, 3.54, 3.60, 4.55, 4.64, 4.73, 4.82, 4.92, 5.02 and 5.13), most of them with at least 115 discrete disks. Appl. Sci. 2020, 10, x 23 of 25 Figure 9. The number of discrete disks and discrete excitation frequencies as a function of the axial depth of cut for the multivariable tool. 5. Conclusions In this work, quadratic and cubic polynomials were used to approximate the delayed terms of delay differential equations. Numerical simulations showed that using secondand third-order EMHPM improved the convergence rate and required less computational time when compared to the first-order EMHPM, and to semi-discretization and full-discretization methods, since fewer approximations or less discrete intervals were needed to reduce the computation time. To further assess the applicability of the proposed method, the third-order EMHPM was used for determining the stability bounds in one-degree-of-freedom milling operation with a multivariable tool, demonstrating that the stability zone improved in comparison with a regular tool. For instance, at 2500 rpm the critical axial depth of the cut was 1.3 mm using the regular milling tool. However, using the multivariable tool, the critical axial depth of the cut was increased until 2.17 mm but more interesting, a stable zone appeared above 8.55 mm. The CWT scalograms, PSD charts and PM were employed to validate the stability lobes found by using the third-order EMHPM for the multivariable tool. Numerical solutions confirmed the system dynamics behavior predicted by the third-order EMHPM. Based on the above results, this paper provided evidence the third-order EMHPM could be used to study dynamic phenomena that appeared at higher axial depths of cut due to the multivariable design of the tool, which broke the excitation frequencies at a lower depth of cut. Author Contributions: Conceptualization, J.d.l.L.S. and D.O.-T.; Methodology, J.d.l.L.S.; Resources, O.M.-R., A.E.-Z. and L.N.L.d.L.; Supervision, D.O.-T., A.E.-Z. and L.N.L.d.L.; Validation, J.d.l.L.S., D.O.-T. and G.U.; Writing—original draft, J.d.l.L.S. and D.O.-T.; Writing—review and editing, O.M.-R., G.U. and A.E.-Z. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Tecnológico de Monterrey through the Research Group of Nanotechnology for Devices Design, and by the Consejo Nacional de Ciencia y Tecnología de México (Conacyt), Project Numbers 242269, 255837, 296176, and National Lab in Additive Manufacturing, 3D Digitizing and Computed Tomography (MADiT) LN299129. Acknowledgments: The authors acknowledge the technical assistance from Oscar Escalera at Tecnologico de Monterrey. Conflicts of Interest: The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results. References 1. Mi, T.; Chen, N.; Stepan, G.; Takacs, D. Energy distribution of a vehicle shimmy system with the delayed tyre model. IFAC-PapersOnLine 2018, 51, 7–12, doi:10.1016/j.ifacol.2018.07.190. 2. Orosz, G.; Stépán, G. Subcritical hopf bifurcations in a car-following model with reaction-time delay. Proc. R. Soc. A Math. Phys. Eng. Sci. 2006, 462, 2643–2670, doi:10.1098/rspa.2006.1660. 3. Hu, H.; Zaihua, W. Dynamics of Controlled Mechanical Systems with Delayed Feedback; Springer Science & Business Media: Berlin, Germany, 2002; ISBN 9783642078392. Figure 9. The number of discrete disks and discrete excitation frequencies as a function of the axial depth of cut for the multivariable tool. 5. Conclusions In this work, quadratic and cubic polynomials were used to approximate the delayed terms of delay differential equations. Numerical simulations showed that using secondand third-order EMHPM improved the convergence rate and required less computational time when compared to the first-order EMHPM, and to semi-discretization and full-discretization methods, since fewer approximations or less discrete intervals were needed to reduce the computation time. To further assess the applicability of the proposed method, the third-order EMHPM was used for determining the stability bounds in one-degree-of-freedom milling operation with a multivariable tool, demonstrating that the stability zone improved in comparison with a regular tool. For instance, at 2500 rpm the critical axial depth of the cut was 1.3 mm using the regular milling tool. However, using the multivariable tool, the critical axial depth of the cut was increased until 2.17 mm but more interesting, a stable zone appeared above 8.55 mm. The CWT scalograms, PSD charts and PM were employed to validate the stability lobes found by using the third-order EMHPM for the multivariable tool. Numerical solutions confirmed the system dynamics behavior predicted by the third-order EMHPM. Based on the above results, this paper provided evidence the third-order EMHPM could be used to study dynamic phenomena that appeared at higher axial depths of cut due to the multivariable design of the tool, which broke the excitation frequencies at a lower depth of cut. Author Contributions: Conceptualization, J.d.l.L.S. and D.O.-T.; Methodology, J.d.l.L.S.; Resources, O.M.-R., A.E.-Z. and L.N.L.d.L.; Supervision, D.O.-T., A.E.-Z. and L.N.L.d.L.; Validation, J.d.l.L.S., D.O.-T. and G.U.; Writing—original draft, J.d.l.L.S. and D.O.-T.; Writing—review and editing, O.M.-R., G.U. and A.E.-Z. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Tecnol ó gico de Monterrey through the Research Group of Nanotechnology for Devices Design, and by the Consejo Nacional de Ciencia y Tecnolog í a de M é xico (Conacyt), Project Numbers 242269, 255837, 296176, and National Lab in Additive Manufacturing, 3D Digitizing and Computed Tomography (MADiT) LN299129. Acknowledgments: The authors acknowledge the technical assistance from Oscar Escalera at Tecnologico de Monterrey. Conflicts of Interest: The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results. References 1. Mi, T.; Chen, N.; Stepan, G.; Takacs, D. Energy distribution of a vehicle shimmy system with the delayed tyre model. IFAC-PapersOnLine 2018,51, 7–12. [CrossRef] 2. Orosz, G.; St é p á n, G. Subcritical hopf bifurcations in a car-following model with reaction-time delay. Proc. R. Soc. A Math. Phys. Eng. Sci. 2006,462, 2643–2670. [CrossRef] Appl. Sci. 2020,10, 7869 21 of 22 3. Hu, H.; Zaihua, W. Dynamics of Controlled Mechanical Systems with Delayed Feedback; Springer Science & Business Media: Berlin, Germany, 2002; ISBN 9783642078392. 4. Altintas, Y. Manufacturing Automation. Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design; Cambridge University Press: New York, NY, USA, 2012. 5. Kuljanic, E.; Totis, G.; Sortino, M. Development of an intelligent multisensor chatter detection system in milling. Mech. Syst. Signal Process. 2009,23, 1704–1718. [CrossRef] 6. Zhuo, Y.; Jin, H.; Han, Z. Chatter identification in flank milling of thin-walled blade based on fractal dimension. Procedia Manuf. 2020,49, 150–154. [CrossRef] 7. Paul, S.; Morales-Menendez, R. Chatter mitigation in milling process using discrete time sliding mode control with type 2-fuzzy logic system. Appl. Sci. 2019,9, 4380. [CrossRef] 8. Peng, C.; Wang, L.; Liao, T.W. A new method for the prediction of chatter stability lobes based on dynamic cutting force simulation model and support vector machine. J. Sound Vib. 2015,354, 118–131. [CrossRef] 9. Li, D.; Cao, H.; Chen, X. Fuzzy control of milling chatter with piezoelectric actuators embedded to the tool holder. Mech. Syst. Signal Process. 2020,148, 107190. [CrossRef] 10. Wan, S.; Li, X.; Su, W.; Yuan, J.; Hong, J.; Jin, X. Active damping of milling chatter vibration via a novel spindle system with an integrated electromagnetic actuator. Precis. Eng. 2019,57, 203–210. [CrossRef] 11. Munoa, J.; Sanz-Calle, M.; Dombovari, Z.; Iglesias, A.; Pena-Barrio, J.; Stepan, G. Tuneable clamping table for chatter avoidance in thin-walled part milling. CIRP Ann. 2020,69, 313–316. [CrossRef] 12. Wang, Y.; Wang, T.; Yu, Z.; Zhang, Y.; Wang, Y.; Liu, H. Chatter prediction for variable pitch and variable helix milling. Shock Vib. 2015,2015. [CrossRef] 13. Mei, Y.; Mo, R.; Sun, H.; He, B.; Bu, K. Stability analysis of milling process with multiple delays. Appl. Sci. 2020,10, 3646. [CrossRef] 14. Wan, M.; Zhang, W.H.; Dang, J.W.; Yang, Y. A unified stability prediction method for milling process with multiple delays. Int. J. Mach. Tools Manuf. 2010,50, 29–41. [CrossRef] 15. Zhan, D.; Jiang, S.; Niu, J.; Sun, Y. Dynamics modeling and stability analysis of five-axis ball-end milling system with variable pitch tools. Int. J. Mech. Sci. 2020,182, 105774. [CrossRef] 16. Slavicek, J. The effect of irregular tooth pitch on stability of milling. Proc. 6th MTDR Conf. 1965,1, 15–22. 17. Budak, E. An analytical design method for milling cutters with nonconstant pitch to increase stability, Part I: Theory. J. Manuf. Sci. Eng. Trans. ASME 2003,125, 29–34. [CrossRef] 18. Budak, E. An analytical design method for milling cutters with nonconstant pitch to increase stability, Part 2: Application. J. Manuf. Sci. Eng. Trans. ASME 2003,125, 35–38. [CrossRef] 19. Altintas, Y.; Engin, S.; Budak, E. Analytical stability prediction and design of variable pitch cutters. J. Manuf. Sci. Eng. Trans. ASME 1999,121, 173–178. [CrossRef] 20. Olgac, N.; Sipahi, R. Dynamics and stability of variable-pitch milling. J. Vib. Control 2007 ,13, 1031–1043. [CrossRef] 21. Jin, G.; Zhang, Q.; Hao, S.; Xie, Q. Stability prediction of milling process with variable pitch cutter. Math. Probl. Eng. 2013,2013. [CrossRef] 22. Comak, A.; Budak, E. Modeling dynamics and stability of variable pitch and helix milling tools for development of a design method to maximize chatter stability. Precis. Eng. 2017,47, 459–468. [CrossRef] 23. Zatarain, M.; Muñoa, J.; Peigne, G.; Insperger, T. Analysis of the influence of mill helix angle on chatter stability. CIRP Ann. Manuf. Technol. 2006,55, 365–368. [CrossRef] 24. Budak, E.; Altintas, Y. Analytical prediction of chatter stability in milling—Part I: General formulation. Am. Soc. Mech. Eng. Dyn. Syst. Meas. Control 1998,120, 22–30. [CrossRef] 25. Insperger, T.; Muñoa, J.; Zatarain, M.; Peign é , G. Unstable islands in the stability chart of milling processes due to the helix angle. In Proceedings of the CIRP 2nd International Conference High Performance Cutting, Vancouver, BC, Canada, 12–13 June 2006. 26. Patel, B.R.; Mann, B.P.; Young, K.A. Uncharted islands of chatter instability in milling. Int. J. Mach. Tools Manuf. 2008,48, 124–134. [CrossRef] 27. Sims, N.D.; Mann, B.; Huyanan, S. Analytical prediction of chatter stability for variable pitch and variable helix milling tools. J. Sound Vib. 2008,317, 664–686. [CrossRef] 28. Turner, S.; Merdol, D.; Altintas, Y.; Ridgway, K. Modelling of the stability of variable helix end mills. Int. J. Mach. Tools Manuf. 2007,47, 1410–1416. [CrossRef] Appl. Sci. 2020,10, 7869 22 of 22 29. Yusoff, A.R.; Sims, N.D. Optimisation of variable helix tool geometry for regenerative chatter mitigation. Int. J. Mach. Tools Manuf. 2011,51, 133–141. [CrossRef] 30. Dombovari, Z.; Stepan, G. The effect of helix angle variation on milling stability. J. Manuf. Sci. Eng. Trans. ASME 2012,134. [CrossRef] 31. Huang, J.; Deng, P.; Li, H.; Wen, B. Stability analysis for milling system with variable pitch cutters under variable speed. J. Vibroeng. 2019,21, 331–347. [CrossRef] 32. Cai, S.; Yao, B.; Feng, W.; Cai, Z.; Chen, B.; He, Z. Milling process simulation for the variable pitch cutter based on an integrated process-machine model. Int. J. Adv. Manuf. Technol. 2020 ,106, 2779–2791. [CrossRef] 33. Olvera, D.; El í as-Z ú ñiga, A.; L ó pez De Lacalle, L.N.; Rodr í guez, C.A. Approximate solutions of delay differential equations with constant and variable coefficients by the enhanced multistage homotopy perturbation method. Abstr. Appl. Anal. 2015, 1–12. [CrossRef] 34. Compe á n, F.I.; Olvera, D.; Campa, F.J.; L ó pez De Lacalle, L.N.; El í as-Z ú ñiga, A.; Rodr í guez, C.A. Characterization and stability analysis of a multivariable milling tool by the enhanced multistage homotopy perturbation method. Int. J. Mach. Tools Manuf. 2012,57, 27–33. [CrossRef] 35. Sims, N.D. Fast chatter stability prediction for variable helix milling tools. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2016,230, 133–144. [CrossRef] 36. Otto, A.; Rauh, S.; Ihlenfeldt, S.; Radons, G. Stability of milling with non-uniform pitch and variable helix Tools. Int. J. Adv. Manuf. Technol. 2017,89, 2613–2625. [CrossRef] 37. Niu, J.; Ding, Y.; Zhu, L.M.; Ding, H. Mechanics and multi-regenerative stability of variable pitch and variable helix milling tools considering runout. Int. J. Mach. Tools Manuf. 2017,123, 129–145. [CrossRef] 38. Olvera, D.; Urbikain, G.; El í as-Zuñiga, A.; L ó pez de Lacalle, L.N. Improving stability prediction in peripheral milling of Al7075T6. Appl. Sci. 2018,8, 1316. [CrossRef] 39. Iglesias, A.; Dombovari, Z.; Gonzalez, G.; Munoa, J.; Stepan, G. Optimum selection of variable pitch for chatter suppression in face milling operations. Materials 2018,12, 112. [CrossRef] [PubMed] 40. Guo, Y.; Lin, B.; Wang, W. Optimization of variable helix cutter for improving chatter stability. Int. J. Adv. Manuf. Technol. 2019,104, 2553–2565. [CrossRef] 41. Hashim, I.; Chowdhury, M.S.H. Adaptation of homotopy-perturbation method for numeric-analytic solution of system of ODEs. Phys. Lett. Sect. A 2008,372, 470–481. [CrossRef] 42. Puma-Araujo, S.D.; Olvera-Trejo, D.; Mart í nez-Romero, O.; Urbikain, G.; El í as-Z ú ñiga, A.; Lacalle, L.N.L. de Semi-active magnetorheological damper device for chatter mitigation during milling of thin-floor components. Appl. Sci. 2020,10, 5313. [CrossRef] 43. Insperger, T.; St é p á n, G. Updated semi-discretization method for periodic delay-differential equations with discrete delay. Int. J. Numer. Methods Eng. 2004,61, 117–141. [CrossRef] 44. Insperger, T. Full-discretization and semi-discretization for milling stability prediction: Some comments. Int. J. Mach. Tools Manuf. 2010,50, 658–662. [CrossRef] 45. Ding, Y.; Zhu, L.M.; Zhang, X.J.; Ding, H. A full-discretization method for prediction of milling stability. Int. J. Mach. Tools Manuf. 2010,50, 502–509. [CrossRef] 46. Lee, E.T.; Eun, H.C. Structural damage detection by power spectral density estimation using output-only measurement. Shock Vib. 2016,2016. [CrossRef] 47. Fries, R.H. Fundamentals of Vibrations; Waveland Press: Long Grove, IL, USA, 2000; ISBN 9781482270372. Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).