Third-order least squares modelling of milling state term for improved computation of stability boundaries

Abstract The general least squares model for milling process state term is presented. A discrete map for milling stability analysis that is based on the third-order case of the presented general least squares milling state term model is first studied and compared with its third-order counterpart that is based on the interpolation theory. Both numerical rate of convergence and chatter stability results of the two maps are compared using the single degree of freedom (1DOF) milling model. The numerical rate of convergence of the presented third-order model is also studied using the two degree of freedom (2DOF) milling process model. Comparison gave that stability results from the two maps agree closely but the presented map demonstrated reduction in number of needed calculations leading to about 30% savings in computational time (CT). It is seen in earlier works that accuracy of milling stability analysis using the full-discretization method rises from first-order theory to second-order theory and continues to rise to the third-order theory. The present work confirms this trend. In conclusion, the method presented in this work will enable fast and accurate computation of stability diagrams for use by machinists.


Introduction
Regenerative chatter has been a cog in the wheel of productive manufacturing of dimensionally accurate components in the machining industry. Machine tool chatter has been encountered and described by Taylor as the most obscure and delicate of all problems facing the machinist (Taylor, 1907). The current understanding that machine tool vibrations are mainly composed of regenerative chatter means that Taylor was more often than not battling with regenerative chatter. The other less-frequent causes of chatter which are respectively, treated in the works (Tlusty & Polacek, 1963;Wiercigroch & Budak, 2001;Wiercigroch & Krivtsov, 2001) are mode coupling effects, frictional effects and thermo-mechanical effects. Regenerative chatter is triggered by random cutting force variation which in turn is caused by regenerative effects. Regenerative effects are the effects of waviness on a machined surface that are caused by disturbances. These disturbances are not deterministic thus may cause KEYWORDS chatter; milling process; least squares approximation; full-discretization method; machine tool chatter stability

OPEN ACCESS
two consecutive tool passes to be out of phase resulting in cutting force variation that excites the tool into growing motion (instability) if the structural damping capacity of the system is not strong enough to dissipate the excitation energy. Works abound in the literature on modelling and stability analysis of regenerative chatter (Altintas, 2001;Altintas & Budak, 1995;Bayly, Halley, Mann, & Davies, 2003;Bayly et al., 2002;Bobrenkov, Khasawneh, Butcher, & Mann, 2010;Ding, Zhu, Zhang, & Ding, 2010a, 2010bFaassen, van de Wouw, Oosterling, & Nijmeijer, 2003;Hanna & Tobias, 1974;Insperger, 2002Insperger, , 2010Insperger & Stépán, 2004;Mann, Bayly, Davies, & Halley, 2004;Merdol & Altintas, 2004;Merrit, 1965;Moon, 1998;Ozoegwu, 2011Ozoegwu, , 2014, 2014, 2016, 2015Patel, Mann, & Young, 2008;Quo, Sun, & Jiang, 2012;Sellmeier & Denkena, 2011;Stepan, 1989;Yi, Nelson, & Ulsoy, 2007). The out-of-process approach that has been widely pursued aims at differentiating the milling process plane of spindle speed and depth of cut into stable and unstable subdomains. The frequency domain methods of Zeroth-Order Approximation (ZOA; Altintas, 2001;Altintas & Budak, 1995) and Multi-Frequency Solution (MFS; Merdol & Altintas, 2004) were proposed by Altintas and co-workers. The ZOA method involves retaining only the first term in the Fourier series expansion of periodic coefficients of the milling model while MFS involves inclusion of higher order harmonics. The ZOA method performed poorly in stability analysis of low-radial immersion milling but the MFS method performed excellently well in both the high and low radial immersion cases. The semi-discretization method is a timedomain method introduced by Insperger (2002). The semi-discretization method is known to predict stability of both the high and low radial immersion milling processes (Insperger & Stépán, 2004). The spatial finite-element concept was adapted to the so-called temporal finite-element analysis (TFEA) for milling stability prediction (Bayly et al., 2002;Patel et al., 2008). The TFEA method originally failed to predict accurate stability at continuous tool immersion (Bayly et al., 2003) and cuts involving simultaneous tooth engagement (Bobrenkov et al., 2010). Increasing the size of Floquet transition matrix of TFEA by finer discretization of discrete delay solved the problem (Bobrenkov et al., 2010;. Some other efficient but not yet popular methods can be seen (Sellmeier & Denkena, 2011;Yi et al., 2007). A recently developed method that is receiving acceptance is the full-discretization method (Ding et al., 2010a). The full-discretization method was seen to be much more computationally faster than the semi-discretization method because the needed matrix exponentials are calculated once at every spindle speed step and not needed to be updated while sweeping the range of the depth of cut at the speed step (Ding et al., 2010a). Linear interpolation theory is used in handling the arising integration problem in (Ding et al., 2010a;Insperger, 2010) to produce a discrete map used for stability analysis. The interpolation theory was later extended to second order (Ding et al., 2010b) and later to third order (Quo et al., 2012). Accuracy improved with increase of order from one to three in the interpolation theory. Use of first-and second-order least squares methods for modelling state term was first explored in the work (Ozoegwu, 2014) in which it was found that it offered reduction in number of needed numerical calculations thus reducing CT relative to the methods of interpolation theory of same order. This beneficial effect was justified in (Ozoegwu, 2014) to grow with order of approximation due to the structure of the arising monodromy matrix. The least squares approximated polynomial for state term of the tool is presented generally and can systematically be extended to higher order cases.
In the present study, the third-order case of the generalized least squares approximated state term is applied in chatter stability analysis and directly compared with that of third-order interpolation theory in (Quo et al., 2012) in order to further strengthen the conclusions of the work (Ozoegwu, 2014).

Mathematical modelling of milling process
The dynamical model shown in Figure 1 is a two degree of freedom (2DOF) depiction of an end-milling tool that vibrates in the x − y plane (horizontal plane). This is the more relevant model when the tool is compliant in both the feed direction (x direction) and feed normal direction (y direction). The modal parametersk x , m x and c x are for x-vibration, while k y , m y and c y are for y-vibration. The tool is symmetric when the x and y modes of vibration have identical parameters. The governing pair of equations of motion becomes It is seen from Figure 1 that the x and y-components of cutting force on the tool are resolved from the normal F n,j (t) and tangential forces F t,j (t) to become The non-linear cutting force law for the 2DOF tool becomes where f a,x and f a,x and f a,y are respective actual feeds in the x and y directions given by and w is depth of cut, C t and C n are the tangential and normal cutting coefficients associated with the workpiece material properties and tool shape and Ӽ is the ratio C n ∕C t . Use of non-linear cutting force law for 1DOF modelling of milling process is directly seen in (Insperger, 2002;Ozoegwu & Omenyi, 2014). The work (Bayly et al., 2002) utilized the linear cutting force law with γ = 1 in its 2DOF modelling of milling process. The vibration of the tool has two components, namely the τ-periodic response due to periodic force of tool-workpiece interaction and the regenerative perturbations. The total motions in the x and y directions are respectively, given as sum of feed motions and vibrations thus where the notations z x (t) and z y (t) are the respective regenerative perturbations in x and y directions, the notations x t (t) and y t (t) are the respective τ-periodic responses in x and y directions and vt is the feed motion in x direction. Inserting equations (3, 4 and 5) in Equation (2) gives where The Q(t) term is expanded in linear Taylor series about v sin j (t) to become Setting z x (t) and z y (t) to zero in the light of Equation (8) The periodic forces −F px (t) and −F py (t), respectively, drive the two orthogonal τ-periodic tool responses x t (t) and y t (t). Putting Equation (9) into (1) and simplifying leads to the coupled 2DOF DDE model for milling process as with the specific periodic cutting force variations Equation (11) put in matrix form reads The coupling in the model which stems from the two-dimensional feed of force law is contained in the τ-periodic matrix because of the presence of non-zero off diagonal elements. Making use of state variables 1 (t) = z x (t) , 2 (t) =̇z x (t) , 3 (t) = z y (t) and 4 (t) =̇z y (t) results in the state space form The natural frequencies and damping ratios of the tool are given in terms of modal parameters k x , m x , c x , k y , m y and c y, respectively, as It should also be noted that when the feed normal motion of the tool is ignored the terms h xy (t), h yx (t) and h yy (t) varnish leaving behind h xx (t) which is now designated with h(t) to give Equation (15) for 1DOF system with.
T . Only the vibration in the feed direction is considered in 1DOF system. The instantaneous angular position of jth tooth j (t) is given as The screen function of the jth tooth is derived in detail in (Ozoegwu, Omenyi, Ofochebe, & Achebe, 2013) to be, respectively, given for up-milling and down-milling in terms of j (t) and ρ (radial immersion) as It should be noted that linearity of cutting force means that γ = 1 which when substituted in equation (12) gives specific force variations identical with those in (Insperger & Stépán, 2004) given directly on the bases of linear force law. The advantage of the form of the presented specific force variations is that non-linear cutting force cases like γ = 0.75 and 0.8 can be investigated straightaway.

The generalized p-order least squares modelling of the milling state term
The discrete delay τ is divided into k equal time intervals leading to equation (15) becoming equation (23) The discretization integer k is normally called the approximation parameter in the literature.
The local delay differential equation (23) is solved by being integrated between the limits t i and t i+1 to become Least squares method is then used to approximate the state term y(s), the delayed state y(s − ) and the time periodic coefficient (s) to allow for evaluation of the Duhamel term of equation (24). Taking a cue from the generalized least squares method given in (Ozoegwu, 2014), it is considered that the state term is approximated as a polynomial y(s) = [a(s)] T b where a(s) = s 0 s 1 s 2 … … s p T is the polynomial of basis vector of p-order is the vector of coefficients also of p-order. Since y(s) is known at every discrete time t i as i the p-order least squares approximation of y(s) results from mini- (23) y(t) = y(t) + (t)y(t) − (t)y(t − ) a s l y l Without these simplifying substitutions introduced at the stage, recommended analytical requirements for handling equation (26) could become prohibitive with rise in order of approximation. This recommendation is the gateway to the use of equation (26) beyond the second order. It should be noted that if the substitutions were made at this stage in the earlier work (Ozoegwu, 2014), second-order least squares modelling of milling state term would have been much shorter.

Re-arranged, Equation (31) becomes
where The discrete map of the system is created from equation (34) to have the monodromy matrix (32g) Since milling discrete map based on state term modelled by third-order interpolation theory exists in the literature (Quo et al., 2012) it becomes pertinent at this point to compare its rate of convergence, accuracy and stability results with the third-order least squares approximated discrete map proposed in this section. This is done in the next two subsections of this section after which higher order least squares approximated discrete maps are considered in the subsequent sections.

Numerical estimate of rate of convergence
The parameters used in all numerical computations in this work which are identical with the experimentally determined parameters of (Bayly et al., 2002) are compiled in Table 1. Table 1 reflects symmetry of the 2DOF tool in which parameters are the same in the feed and feed normal directions. The properties of the laptop computer used for computations are as follows: manufacturer; Hewlett-Packard, Model; Presario CQ61 notebook PC, Rating; 3.9 (Windows Experience Index), Prosessor; AMD Sempron(tm) M120 2.10 GHz, Installed Memory (RAM); 2.00 GB (1.75 GB usable), System type; 64-bit Operating system. natural frequency of 1dof tool ω n 5793 rads −1 Mass of 1dof tool m 0.03993 kg natural frequencies of 2dof symmetric tool ω nx = ω ny 5793 rads −1 Mass of 2dof tool m x = m y 0.03993 kg damping ratio of 1dof tool ζ 0.011 damping ratio of 2dof tool ζ x = ζ y 0.011 tangential cutting coefficient C t 6 × 10 8 nm −1−γ tangential to normal cutting coefficient ratio Ӽ 2/6 feed exponent in the cutting force law γ 1 number of teeth . it is seen that the 3rd order least squares approximation (3rd order lsa) exhibits the fastest relative convergent followed by 2nd order lsa then the 1st order lsa. it should also be noted that the convergence of the presented 3rd order lsa is very similar with that of Quo et al. (2012). order LSA, 2nd order LSA and 3rd order LSA, respectively, suggesting that relative to map of the interpolation theory of same order savings in CT rise with n. This has been shown to be true in (Ozoegwu, 2014) where 1st order LSA and 2nd order LSA, respectively, saved 13 and 21% of CT. It will be seen in what follows that the presented 3rd order LSA saves more CT than 2nd order LSA as expected.
The monodromy matrix must have eigenvalues existing on the unit circle centred at the origin of the complex plane for stability. The stability lobes are numerically tracked based on this stability criterion and plotted as critical depth of cut as a function of critical spindle speed. This is done by computing eigenvalues of with k = 40 on a 200 by 100 gridded plane of spindle speed and depth of cut. The generated stability diagrams for the down-milling mode (using Equation (22) for screen function) are presented in Figure 4. Figure 4 should be compared with those in the literature (Ding et al., 2010b;Insperger & Stépán, 2004) that are generated using identical set of parameters to see the close similarity and thus the applicability of the proposed 3rd order LSA map. The CT for each diagram using the currently proposed 3rd order LSA is contained in the caption. Using the form of milling models presented in this work in the method proposed by Quo et al. (2012), the equivalent stability diagrams are computed and given with CT's in Figure 5. The MATLAB programs written for both the presented 3rd order LSA map and the map of (Quo et al., 2012) are identical except that the relevant matrices are inserted at appropriate points to avoid bias. It is seen that the two methods produce similar stability boundaries, while 3rd order LSA computes much faster. Relative to the map of (Quo et al., 2012), the presented discrete map leads to about 30% gain in CT. It is seen that as expected the presented 3rd . it is seen that 3rd order lsa exhibits the fastest asymptotic behaviour thus the most convergent followed by 2nd order lsa then the 1st order lsa.  order LSA saves more CT than the 2nd order LSA that was seen in (Ozoegwu, 2014) to save 21% of CT. It is even more striking to note that against the general notion that an increase in the order of interpolation causes increase in CT (Ding et al., 2010b;Quo et al., 2012), the proposed 3rd order LSA map is faster than the first-order and secondorder interpolation-based maps proposed by Ding et al. (2010a,b), respectively. For example, Figure 4(a) and (c) took 442 and 413 s, respectively, but the equivalent figures are computed from the method of (Ding et al., 2010a) to take 488 and 468 s, respectively, and also computed from the method of (Ding et al., 2010b) to take 556 and 531 s, respectively. The presented method further consolidates the earlier notion (Ozoegwu, 2014) that least squares modelling of milling state term is powerful for minimizing CT and gets even more powerful in this respect with rise in order of approximation. The generated stability diagrams for the up-milling mode (using Equation (21) for screen function) are presented in Figure 6. Application of the presented map to stability analysis of 2DOF down-milling is demonstrated by generating Figure 7, while application to stability analysis of 2DOF up-milling is demonstrated in Figure 8. Wu, Luo, and Zhang (2013) provided an instructive screening procedure that identifies data points close to the stability limit (the so-called valid region) and used the fulldiscretization method of Ding et al. (2010a) in computing the monodromy operators at these data points. They reorganized computational itinerary of an existing method by eliminating seemly unnecessary computations with the aim of saving computational time. It is shown in the work of Ozoegwu, (2014) that the least squares approximated full-discretization methods save considerable time relative to the interpolation-based full-discretization methods presented in the works (Ding et al., 2010a,b;Insperger, 2010;Quo et al., 2012) thus anyone implementing the screening procedure of Wu et al. (2016) and uses the least squares approximated full-discretization methods instead of their interpolation-based counterparts in computations within the valid region will much further save computational time and thus further enhance speed of machining parameter selection in the workshop.

Conclusion
The general least squares model for milling process state term is presented. This kind of model is used in discrete mapping needed for stability analysis of milling process chatter. The effectiveness and computational time (CT) of the map generated using third-order case of the proposed general state term model are judged against that of Quo et al. (2012) that is based on the third-order interpolation theory. A gain of about 30% in CT is demonstrated by the proposed map. This beneficial result derives from the fact that the presented method results in monodromy matrix that is formed with less number of computations. It should be noted that the first-and second-order least squares approximation theory, respectively, saved 13 and 21% of CT thus the resulting 30% savings from the third-order theory are consistent with the conclusion in (Ozoegwu, 2014) that CT savings grow with order of approximation. Numerical rate of convergence analysis (NRCA) of both 1DOF and 2DOF systems is carried out. The NRCA shows that convergence of the 1DOF system using the presented third-order map is very similar with that of the map (Quo et al., 2012). The major Figure 7. stability lobes of 2dof down-milling mode of the system from the proposed third-order least squares map at (a) radial immersion ρ = 1 and (b) ρ = 0.05. these results should be compared with the equivalent stability diagrams in (Quo et al., 2012) to see close similarity. contribution of this work to the production and manufacturing industry is that a method of milling stability analysis that saves 30% of computational time of the state of art with full retention of accuracy is presented. This will allow faster stability characterization of milling process that will provide the machinist greater freedom in systematic selection of milling process parameters.

Disclosure statement
No potential conflict of interest was reported by the authors.