A general order full-discretization algorithm for chatter avoidance in milling

Abstract

Based on the full-discretization method, this work presents a generalized monodromy matrix as an exact function of the order of polynomial approximation of the milling state for chatter avoidance algorithm. In other words, the computational process is made smarter since the usual derivation of the monodromy matrices on order-by-order basis – a huge analytical involvement that rapidly gets heavier with a rise in the order of approximation – is bypassed. This is the highest possible level of generalization that seems to be the first of its kind among the time-domain methods as the known generalizations are limited to the interpolating/approximating polynomial of the milling state. It then became convenient in this work to study the stability of milling process up to the tenth order. More reliable methods of the rate of convergence analysis were suggested and utilized in consolidating the known result that the best accuracy of the full-discretization method lies with the third and fourth order. It is seen from numerical convergence analyses that, although accuracy most often decreases with rising order beyond the third-order methods, the trend did not persist with a continued rise in order.

Introduction

Regenerative vibration is a violent and noisy chatter of machine tools that limits material removal rate. It causes scraping of parts due to its inscription of vibration marks on the parts. Such compromise of integrity leads to a rise in economic loses. Regenerative chatter is self-excited and it occurs when an unfavourable phase difference occurs between two successive cuts thus it is governed by delay differential equation models. The regenerative delay differential equation models for turning process are autonomous^1–9 while those of milling process are periodically non-autonomous.^10–13 The stability analyses of milling processes are either in the frequency domain or time domain. Frequency-domain approaches are based on truncated Fourier series expansion of the periodic coefficients and utilization of the Nyquist stability criterion.^14–18 A frequency domain approach that is based on utilizing only the first term in the Fourier series expansion of the periodic coefficients is called the zeroth-order approximation (ZOA) method. ¹⁶ The ZOA method fails to reveal the flip bifurcation that appears in low radial immersion milling.^19–21 Another frequency domain approach that overcame the shortcomings of ZOA is based on utilizing additional higher order harmonics in the Fourier series expansion of the periodic coefficients. It is called the multi-frequency solution method. ¹³ Insperger and colleagues^22–24 introduced a time-domain approach called the semi-discretization method which was later improved for computational efficiency and accuracy by Henninger and Eberhard. ²⁵ The semi-discretization method basically involves discretizing the delayed term while leaving the undelayed terms undiscretized and approximating the periodic coefficient matrices as piecewise constant functions thus allowing the formation of ordinary differential equations from which a monodromy matrix is constructed. Further applications of the semi-discretization method can be seen.^26–31 The semi-discretization method has been recently reconstructed with Shannon standard orthogonal basis by Dong et al. ³² to generate an equally accurate but much faster method. Temporal finite element analysis was introduced for milling process stability analysis by Bayly et al. ³³ The method has been applied in the stability analysis of milling process in other works.^34–37 A cluster treatment of characteristic roots for stability analysis of autonomous delayed systems and regenerative machining was developed by Olgac and Sipahi.^38–40 Chebyshev collocation method was introduced for stability analysis of milling by Butcher et al.^41,42 Shifted Chebyshev polynomials were very recently applied in the prediction of stability of milling thin-walled workpiece by Yan et al. ⁴³ Ding et al. ⁴⁴ introduced the method of full-discretization which discretized the delayed state like the semi-discretization method and further discretized the current state. The structure of the full-discretization method (FDM) allows the reduction of computational time relative to the semi-discretization method. The FDM, as applied in the stability analysis of milling, can be classified into two: the methods based on the interpolation theory^44–49 together with the more-time-saving methods based on the approximation (from the least squares sense) theory^50,51 and the methods based numerical integration which are further sub-divided into those based on Newton–Cotes method^52–54 and those based on spectral method. ⁵⁵ The FDM has been applied in the study of different aspects of milling, for example, the least squares version has been applied^56,57 in minimizing pocketing time. Complete discretization method for prediction of milling stability was introduced by Li et al. ⁵⁸ The method not only discretizes the delayed state and the current state like the FDM but further discretizes other terms. An improved complete discretization method has been proposed very recently. ⁵⁹ A generalized Runge–Kutta method based on Volterra integral equation of the second kind has been proposed by Niu et al. ⁶⁰ A wavelet-based approach has been proposed for the stability analysis of periodic delay differential equations with discrete delay and application of the method demonstrated with the milling process by Ding et al. ⁶¹

Detailed review of the development of the FDM, as seen in the works,^44–51 shows that each of the available first- to fifth-order methods were derived from the first principles. The derivation process requires a computation of a polynomial of needed order, an evaluation of a vector/matrix integration problem and a manipulation of the arising results towards constructing a monodromy matrix. This process rapidly gets too complicated and unwieldy with rising order. Increase in order was necessitated by the initial recognition that it causes a rise in the accuracy of the method. The aim of this work is to derive the FDM of general-order p such that a monodromy matrix of general-order results, preempting the difficult derivation process in future use of the method. Such generalized monodromy matrix could then be written down at the specification of the order of approximation. The proposed generalization is carried out on the first class of the FDMs that are based on the interpolation/approximation theories. Specifically, the generalization is carried out on the methods based on the approximation theory because they are of same accuracy as those based on the interpolation theory but considerably saves more computational time. ⁵⁰

A mathematical model of milling process

The milling models presented in this section are adaptations from existing works such as Ozoegwu et al. ⁵¹ The model that conforms with the stability analysis based on the FDM is

$$\stackrel{.}{y}\left(t\right)=Ay\left(t\right)+B\left(t\right)y\left(t\right)-B\left(t\right)y\left(t-\tau \right)$$

(1)

where the matrices and the state vector for the 1 degree of freedom (1DOF) model has the form

$$A=\left[\begin{array}{cc}0& 1\\ -{\omega }_n^2& -2\zeta {\omega }_n\end{array}\right]$$

(2)

$$B\left(t\right)=\left[\begin{array}{cc}0& 0\\ -\frac{wh\left(t\right)}{m}& 0\end{array}\right]$$

(3)

$$y\left(t\right)={\left\{\begin{array}{cc}z\left(t\right)& \stackrel{.}{z}\left(t\right)\end{array}\right\}}^T$$

(4)

The symbols $z(t)$ and $\stackrel{·}{z}(t)$ respectively represent the regenerative displacement and regenerative velocity of the tool system in the feed direction. The specific force variation $h(t)$ is given as

$$h\left(t\right)=C_t\gamma {\left(v\tau \right)}^{\gamma -1}{\displaystyle \sum _{j=1}^N{g}_j}\left(t\right){\sin }^{\gamma }{\theta }_j\left(t\right)\left[X\sin {\theta }_j\left(t\right)+\cos {\theta }_j\left(t\right)\right]$$

(5)

The symbols ${\omega }_n$ and $\zeta$ respectively represent the natural frequency and damping ratio of the tool system in the feed direction. The meaning of the other symbols can be seen in Table 1.

Table 1.

The parameters used in all numerical computations.

Parameter	Value
Natural frequency of 1DOF tool ${\omega }_n$	$5793\mathrm{rad}{\mathrm{s}}^{-1}$
Mass of 1DOF tool m	$0.03993\mathrm{kg}$
Damping ratio of 1DOF tool $\zeta$	$0.011$
Feed exponent in the cutting force law $\gamma$	$1$
Tangential cutting coefficient $C_t$	$6\times {10}^8\mathrm{N}{\mathrm{m}}^{-1-\gamma }$
Tangential to normal cutting force coefficient ratio X	0.33333
Number of teeth N	$2$

DOF: degree of freedom.

These parameters are identical with the experimentally determined parameters of Bayly et al. ³³

The symbol w, as seen in the non-zero element of $\mathbf{B}(t)$ represents axial depth of cut. The angular position of $j\mathrm{th}$ tooth ${\theta }_j(t)$ is given by

$${\theta }_j\left(t\right)=\left(\frac{\pi \Omega }{30}\right)t+\left(j-1\right)\frac{2\pi }{N}$$

(6)

where $N$ is the number of flutes and $\Omega$ is the spindle speed (r/min). The screen function $g_j(t)$ has values of either one or zero depending on whether the tool is cutting or not. The screen function of $j\mathrm{th}$ tooth is given for up-milling and down-milling in terms of ${\theta }_j(t)$ and $\rho$ (radial immersion) as

$$g_j\left(t\right)=\frac{1}{2}\{1+\mathrm{sgn}[\sin ({\theta }_j\left(t\right)-\arctan \left\{\frac{-1}{2\rho }\sin \left[\arccos (1-2\rho )\right]\right\})+\sin \left(\arctan \left\{\frac{-1}{2\rho }\sin \left[\arccos (1-2\rho )\right]\right\}\right)]\}$$

(7)

$$\begin{array}{l}{g}_j\left(t\right)=\frac{1}{2}\\ \left\{1+\mathrm{sgn}\left[\sin \left({\theta }_j\left(t\right)-\arctan \left\{\frac{1}{2\rho }\sin \left[\arccos \left(2\rho -1\right)\right]\right\}\right)-\sin \left(\arccos \left(2\rho -1\right)-\arctan \left\{\frac{1}{2\rho }\sin \left[\arccos \left(1-2\rho \right)\right]\right\}\right)\right]\right\}\end{array}$$

(8)

A generalized order FDM

The discrete delay $\tau$ of equation (1) is discretized into equal k time intervals $\Delta t$ . The discrete interval is given as $\Delta t=\tau /k=t_{i+1}-t_i$ where $i=0,1,2,\dots ,(k-1)$ , $t_i=i\tau /k=i\Delta t=i(t_{i+1}-t_i)$ . Equation (1) is solved using direct integration scheme to give

$$y_{i+1}=e^{A\Delta t}{y}_i+{\displaystyle \underset{t_i}{\overset{t_{i+1}}{\int }}{e}^{A\left(t_{i+1}-s\right)}}\left[B\left(s\right)y\left(s\right)-B\left(s\right)y\left(s-\tau \right)\right]\mathrm{d}s$$

(9)

The general p-order least squares approximation of $\mathbf{y}(t)$ is seen in Ozoegwu and colleagues^51,62 to be

$$\mathbf{y}\left(t\right)={\left[\mathbf{a}\left(t\right)\right]}^T{\left\{\sum _{l=i-p+1}^{i+1}\mathbf{a}\left(t_l\right){\left[\mathbf{a}\left(t_l\right)\right]}^T\right\}}^{-1}\sum _{l=i-p+1}^{i+1}\mathbf{a}\left(t_l\right){\mathbf{y}}_l$$

(10)

The matrix $\mathbf{a}(t_l)[\mathbf{a}(t_l){]}^T$ and vector $\mathbf{a}(t_l){\mathbf{y}}_l$ are expanded and the substitutions, $t_{i-r}=-r\Delta t$ , $r=p-1,p-2,\dots ,0,-1$ , introduced to give

$$\begin{array}{c}\sum _{l=i-p+1}^{i+1}\mathbf{a}\left(t_l\right){\left[\mathbf{a}\left(t_l\right)\right]}^T\hfill \\ =\left[\begin{array}{ccccc}p\hfill & \sum _{\alpha =0}^{\alpha =p}(1+\alpha -p)\Delta t\hfill & \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^2\Delta t^2\hfill & \dots \hfill & \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^p\Delta t^p\hfill \\ \multicolumn{1}{c}{\sum _{\alpha =0}^{\alpha =p}(1+\alpha -p)\Delta t}& \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^2\Delta t^2\hfill & \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^3\Delta t^3\hfill & \dots \hfill & \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^{p+1}\Delta t^{p+1}\hfill \\ \multicolumn{1}{c}{\sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^2\Delta t^2}& \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^3\Delta t^3\hfill & \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^4\Delta t^4\hfill & \dots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & \dots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^p\Delta t^p}& \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^{p+1}\Delta t^{p+1}\hfill & \dots \hfill & \dots \hfill & \sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^{2p}\Delta t^{2p}\hfill \end{array}\right]\hfill \end{array}$$

(11a)

$$\sum _{l=i-p+1}^{i+1}\mathbf{a}\left(t_l\right){\mathbf{y}}_l=\left\{\begin{array}{c}\sum _{\alpha =0}^{\alpha =p}{\mathbf{y}}_{i+1+\alpha -p}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}\sum _{\alpha =0}^{\alpha =p}(1+\alpha -p)\Delta t{\mathbf{y}}_{i+1+\alpha -p}\hfill \\ \multicolumn{1}{c}{\sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^2\Delta t^2{\mathbf{y}}_{i+1+\alpha -p}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{\alpha =0}^{\alpha =p}{(1+\alpha -p)}^p\Delta t^p{\mathbf{y}}_{i+1+\alpha -p}}\end{array}}\end{array}\right\}$$

(11b)

$$\mathbf{a}\left(t\right)=\left\{\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{t}\\ \multicolumn{1}{c}{\begin{array}{c}{t}^2\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{t^p}\end{array}}\end{array}\right\}$$

(11c)

Equations (11a)-(11c) are introduced in equation (10) and evaluated to give the general form

$$\mathbf{y}\left(t\right)=\sum _{q=1-p}^{q=1}\left[\frac{-1^{1-q}}{(-q+1)!(p+q-1)!\Delta t^p}\underset{j=-1,j\ne -q}{\overset{p-1}{\Pi }}(t+j\Delta t)\right]{\mathbf{y}}_{i+q}$$

(12)

where $q=1-p,2-p,\dots ,-1,0,1$ . The simple structure of any least squares approximated state term identified in Ozoegwu and colleagues^50,51,62 is obeyed by equation (12). The simple structure allows the state term assume the value of the discrete state that corresponds to the substituted discrete time. This criterion is used to guarantee accuracy of the analysis leading to equation (12). The delayed state $\mathbf{y}(t-\tau )$ and the periodic coefficient matrix $\mathbf{B}(t)$ are basically modelled linear with $p=1$ in equation (12). In adapting equation (12) to the delayed state and the periodic coefficient matrix, the discrete states ${\mathbf{y}}_{i+q}$ are replaced with ${\mathbf{y}}_{i+q-k}$ and ${\mathbf{B}}_{i+q}$ , respectively to give

$$y\left(t-\tau \right)={\displaystyle \sum _{q=0}^{q=1}\left[\frac{-1^{1-q}}{\left(-q+1\right)!\left(q\right)!\Delta t}{\displaystyle \prod _{j=-1,j\ne -q}^0\left(t+j\Delta t\right)}\right]}{y}_{i+q-k}$$

(13)

$$\mathbf{B}\left(t\right)=\sum _{q=0}^{q=1}\left[\frac{-1^{1-q}}{(-q+1)!\left(q\right)!\Delta t}\underset{j=-1,j\ne -q}{\overset{0}{\Pi }}(t+j\Delta t)\right]{\mathbf{B}}_{i+q}$$

(14)

Equations (12)–(14) are inserted into equation (9) and the integration carried out between the limits $t_i$ and $t_{i+1}$ to give

$${\mathbf{y}}_{i+1}={\mathbf{P}}_i({\mathbf{F}}_0+{\mathbf{G}}_{p,p}{\mathbf{B}}_i+{\mathbf{G}}_{p,2p+1}{\mathbf{B}}_{i+1}){\mathbf{y}}_i+{\mathbf{P}}_i\sum _{q=-1}^{q=1-p}({\mathbf{G}}_{p,p+q}{\mathbf{B}}_i+{\mathbf{G}}_{p,2p+1+q}{\mathbf{B}}_{i+1}){\mathbf{y}}_{i+q}-{\mathbf{P}}_i({\mathbf{G}}_{1,2}{\mathbf{B}}_i+{\mathbf{G}}_{1,4}{\mathbf{B}}_{i+1}){\mathbf{y}}_{i+1-k}-{\mathbf{P}}_i({\mathbf{G}}_{1,1}{\mathbf{B}}_i+{\mathbf{G}}_{1,3}{\mathbf{B}}_{i+1}){\mathbf{y}}_{i-k}$$

(15)

where for $q=0$

$$G_{p,p+q}=\frac{{\displaystyle {\sum }_{l=0}^{l=p+1}{\Delta }_{p,p+q}^{\left(l\right)}}{F}_{p+2-l}}{\left(-q+1\right)!\left(p+q-1\right)!\Delta t^{p+1}}=\frac{{\displaystyle {\sum }_{l=0}^{l=p+1}{\Delta }_{p,p+q}^{\left(l\right)}}{F}_{p+2-l}}{\left(p-1\right)!\Delta t^{p+1}}$$

(16a)

$$G_{p,2p+1+q}=\frac{{\displaystyle {\sum }_{l=0}^{l=p}{\Delta }_{p,2p+1+q}^{\left(l\right)}}{F}_{p+2-l}}{\left(-q+1\right)!\left(p+q-1\right)!\Delta t^{p+1}}=\frac{{\displaystyle {\sum }_{l=0}^{l=p}{\Delta }_{p,2p+1+q}^{\left(l\right)}}{F}_{p+2-l}}{\left(p-1\right)!\Delta t^{p+1}}$$

(16b)

and for $q\ne 0$

$${\mathbf{G}}_{p,p+q}=\frac{{\sum }_{l=0}^{l=p}{\Delta }_{p,p+q}^{\left(l\right)}{\mathbf{F}}_{p+2-l}}{(-q+1)!(p+q-1)!\Delta t^{p+1}}$$

(16c)

$${\mathbf{G}}_{p,2p+1+q}=\frac{{\sum }_{l=0}^{l=p-1}{\Delta }_{p,2p+1+q}^{\left(l\right)}{\mathbf{F}}_{p+2-l}}{(-q+1)!(p+q-1)!\Delta t^{p+1}}$$

(16d)

The symbol ${\Delta }_{p,2p+1+q}^{\left(l\right)}$ is used to represent the $l\mathrm{th}$ term of any of the quantities

$${\Delta }_{p,2p+1+q}=\sum _{j=1}^{p-1}(1+j\Delta t),q=1$$

(17a)

$${\Delta }_{p,2p+1+q}=(-1{)}^{-q}(-1+\Delta t)\sum _{j=1,j\ne -q}^{p-1}(1+j\Delta t),q=0,-1,-2,\dots ,1-p$$

(17b)

$${\Delta }_{p,p+q}=\left(-1+\Delta t\right){\Delta }_{p,2p+1+q},q=1,0,-1,-2,\dots ,1-p$$

(17c)

$${\mathbf{F}}_{p+2-l(=p+2)}={\mathbf{F}}_0=e^{\mathbf{A}\Delta t}$$

(18a)

$${\mathbf{F}}_{p+2-l(=p+1)}={\mathbf{F}}_1=({\mathbf{F}}_0-\mathbf{I}){\mathbf{A}}^{-1}$$

(18b)

$${\mathbf{F}}_{p+2-l}=[(p+1-l){\mathbf{F}}_{p+1-l}-\Delta t^{p+1-l}\mathbf{I}]{\mathbf{A}}^{-1},l=p,p-1,\dots ,0$$

(18c)

$${\mathbf{P}}_i^{\left(p\right)}={[\mathbf{I}-{\mathbf{G}}_{p,p+1}{\mathbf{B}}_i-{\mathbf{G}}_{p,2p+2}{\mathbf{B}}_{i+1}]}^{-1}$$

(19)

The monodromy matrix for the system is constructed as

$${\mathbf{ѱ}}^{\left(p\right)}={\mathbf{M}}_{k-1}^{\left(p\right)}{\mathbf{M}}_{k-2}^{\left(p\right)}\dots {\mathbf{M}}_0^{\left(p\right)}$$

(20)

where

$$M_{i\left(=k-1,k-2,\dots ,1,0\right)}^{\left(p\right)}=\left[\begin{array}{cccccccccc}{M}_{1,1}^{\left(i,p\right)}& M_{1,2}^{\left(i,p\right)}& \cdots & M_{1,p-1}^{\left(i,p\right)}& M_{1,p}^{\left(i,p\right)}& 0& \cdots & 0& M_{1k}^{\left(i,p\right)}& M_{1,k+1}^{\left(i,p\right)}\\ I& 0& \cdots & 0& 0& 0& \cdots & 0& 0& 0\\ 0& I& \cdots & 0& 0& 0& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & 0& 0& 0& 0& 0& I& 0\end{array}\right]$$

(21)

and

$${\mathbf{M}}_{1,1}^{(i,p)}={\mathbf{P}}_i^{\left(p\right)}({\mathbf{F}}_0+{\mathbf{G}}_{p,p}{\mathbf{B}}_i+{\mathbf{G}}_{p,2p+1}{\mathbf{B}}_{i+1})$$

(22a)

$${\mathbf{M}}_{1,1-q}^{(i,p)}={\mathbf{P}}_i^{\left(p\right)}({\mathbf{G}}_{p,p+q}{\mathbf{B}}_i+{\mathbf{G}}_{p,2p+1+q}{\mathbf{B}}_{i+1}),q=-1,-2,\dots ,1-p$$

(22b)

$${\mathbf{M}}_{1k}^{(i,p)}=-{\mathbf{P}}_i^{\left(p\right)}({\mathbf{G}}_{1,2}{\mathbf{B}}_i+{\mathbf{G}}_{1,4}{\mathbf{B}}_{i+1})$$

(22c)

$${\mathbf{M}}_{1,k+1}^{(i,p)}=-{\mathbf{P}}_i^{\left(p\right)}({\mathbf{G}}_{1,1}{\mathbf{B}}_i+{\mathbf{G}}_{1,3}{\mathbf{B}}_{i+1})$$

(22d)

The summary of generating the general p-order least squares approximated full-discretization method (LSAFDM) simply involves

1. Making use of equation (17) and equations (18a)-(18c) to generate the ( $p+2$ ) different $\mathbf{G}$ matrices in equations (16a)-(16d).

2. The $\mathbf{G}$ matrices are used to generate ${\mathbf{P}}_i^{(p)}$ according to equation (19).

3. The first-row sub-matrices of the local transition matrices ${\mathbf{M}}_{i(=k-1,k-2,\dots ,1,0)}^{(p)}$ are then computed using equations (22a)-(22d) where ${\mathbf{M}}_{i(=k-1,k-2,\dots ,1,0)}^{(p)}\equiv {\mathbf{M}}_i^{(p)}$ for $i=k-1,k-2,\dots ,1,0$ .

4. Finally the monodromy matrix ${\mathbf{ѱ}}^{(p)}$ is formed.

Application in deriving known first- to fifth-order FDMs

Generating LSAFDM at a given numerical value of p previously required challenging analytical processes of deriving the p-order polynomial and evaluating the symbolic integration problem needed to get to the discrete transition. On the traditional order-by-order basis, FDM requires discretization, modelling of state variation between mesh points through approximation/interpolation, solution of complex integration problem and manipulation of results to get the G matrices used to build the ${\mathbf{M}}_i^{(p)}$ matrices that are in turn needed to form the monodromy matrix ${\mathbf{ѱ}}^{(p)}$ . The approximation/interpolation and solution of the complex integration problem are unarguably the most challenging repetitive steps that get rapidly more complicating and discouraging (in fact, the rise in complexity is more than geometric since the number of G matrices increases from 3 to 2n + 5 when the order of approximation increases from 1 to n) with rise in order. Therefore, in terms of symbolic analysis, the generic p-order formulation is much more compact and more efficient to handle. By virtue of this generalization, the monodromy matrix of LSAFDM of general p-order can directly be written down without encountering the challenging analytical processes. This work marks a strong point of improvement by deriving the general p-order monodromy. This generalization removes the discouraging heavy analytical requirement that prevented analysis of very high-order FDMs. In order to demonstrate the use of the presented generalization, the known first- to fifth-order FDMs are directly written down following the four steps of the generalization given above; details are given for the second- and fourth-order LSAFDM as presented in Appendices 1 and 2. These monodromy matrices of the second- and fourth-order LSAFDM from the generalized order LSAFDM are seen to be identical with that derived from the first principle in Ozoegwu ⁵⁰ and Ozoegwu et al., ⁵¹ respectively. Attempt to derive the fourth-order LSAFDM from the first principle, as done in Ozoegwu et al., ⁵¹ will highlight the considerable analytical simplification resulting from the presented generalization.

Application in deriving hyper-fifth-order FDMs

Investigation of FDM beyond fifth order is carried out with the presented generalization. The process of generalization, though rigorous, has truncated the most difficult analytical steps (generation of relevant polynomial and handling of a convoluted integration problem) in the generation of FDM of any given order such that any needed FDM can directly be written down. Making use of the presented general p-order FDM, the sixth- to tenth-order LSAFDMs are directly written down (as given in Table 4 in Appendix 3) and numerically studied in the following sections.

Error and convergence analysis

The 1DOF case

The estimate of numerical rate of convergence with rising approximation parameter (ENRCAP) of maps has been determined by considering how fast the function $f(k)=||{\mu }_{\mathrm{sr}}(k)|-\left|{\mu }_{\mathrm{rsr}}\right||$ converges to zero^45–48,51 where $f(k)$ is size of deviation of modulus of spectral radius ${\mu }_{\mathrm{sr}}(k)$ at approximation parameter k from the modulus of the reference spectral radius ${\mu }_{\mathrm{rsr}}$ of high accuracy at high value (reference) of k designated $k_R$ . The more accurate method is expected to exhibit $f(k)$ that vanishes faster. That is, the more accurate method has higher rate of convergence with the rise of approximation parameter. It is slightly misleading to consider $f(k)=||{\mu }_{\mathrm{sr}}(k)|-|{\mu }_{\mathrm{rsr}}(k_R)||$ as size of deviation of spectral radius from the reference. This is illustrated in Figure 1 where a pair of spectral radii ${\mu }_{\mathrm{sr}1}(k)=p_{1x}+j{p}_{1y}$ and ${\mu }_{\mathrm{sr}2}(k)=p_{2x}+j{p}_{2y}$ , located at the point of intersection of two circles and a straight line, give $f(k)=0$ relative to a reference spectral radius ${\mu }_{\mathrm{rsr}}=a+\mathrm{jb}$ at the centre of the smaller circle while any other spectral radii ${\mu }_{\mathrm{sr}}(k)=x_1+j{y}_1$ at the circumference of the smaller circle gives $f(k)>0$ . A more accurate representation of size of deviation of spectral radius from the reference spectral radius is $r(k)=|{\mu }_{\mathrm{sr}}(k)-{\mu }_{\mathrm{rsr}}(k_R)|$ which is seen to give same deviation for all points on the circumference of the smaller circle. It is thus suggested, based on the ongoing argument, that $r(k)=|{\mu }_{\mathrm{sr}}(k)-{\mu }_{\mathrm{rsr}}|$ be used in determining ENRCAP instead of the hitherto used $f(k)=||{\mu }_{\mathrm{sr}}(k)|-|{\mu }_{\mathrm{rsr}}(k_R)||$ . An alternative, that could be in the form of normalized error, is $r_n(k)=|({\mu }_{\mathrm{sr}}(k)-{\mu }_{\mathrm{rsr}})/{\mu }_{\mathrm{rsr}}|$ or, better still, the percentage error $100r_n(k)$ . Using the parameters in Table 1, the $r(k)$ is plotted in Figure 2(a) at Ω = 5000 r/min and w = 0.0002 for the k range of $[k_s,k_e]=[20,100]$ . The reference spectral radius ${\mu }_{\mathrm{rsr}}$ is computed using $k_R=500$ in the third-order FDM. The judgement of convergence could be said to be clear below k = 25 where the fourth-order FDM seems most convergent but not above k = 25 except for the first-order FDM clearly seen to be of least convergence for k above 25. A plot in the k range $[k_s,k_e]=[40,100]$ as shown in Figure 2(b) basically gives a different judgement, tending to suggest the third-order FDM as the most convergent. It is seen from Figure 2(c) and (d) that same results are gotten when based on $100r_n(k)$ rather than $r(k)$ .

Figure 1.

Illustrating the superiority of $r(k)$ over $f(k)$ .

Figure 2.

(a) Plot of $r(k)$ in the k range of $[k_s,k_e]=[20,100]$ , (b) plot of $r(k)$ in the k range of $[k_s,k_e]=[40,100]$ , (c) plot of $100r_n(k)$ in the k range of $[k_s,k_e]=[20,100]$ and (d) plot of $100r_n(k)$ in the k range of $[k_s,k_e]=[40,100]$ .

The ENRCAP based on $f(k)$ has been judged visually in the works^45–48,51 as done with $r(k)$ or $100r_n(k)$ here. Visual judgements could lack clarity especially at higher values of k where $f(k)$ , $r(k)$ or $100r_n(k)$ get packed together. It is reasoned here that it is better to compare various orders of approximation with summed deviation given as $S_r^{(k_R)}(p)=\sum _{k=k_s}^{k_e}r(k)$ where $k_s$ is the start value of k and $k_e$ is the end value. The $S_r^{(k_R)}(p)$ , which could be called estimate of summed deviations from the reference spectral radii (ESDRSR), will give an idea of performance of the FDMs of different orders p over a given range $[k_s,k_e]$ of k. The most converged method will exhibit least value for $S_r^{(k_R)}(p)$ . An alternative ESDRSR which should also be noted is $S_r^{(k_e)}(p)=\sum _{k=k_s}^{k_e}|{\mu }_{\mathrm{sr}}(k)-{\mu }_{\mathrm{rsr}}(k_e)|$ . The difference between $S_r^{(k_R)}(p)$ and $S_r^{(k_e)}(p)$ is clear; the deviation for $S_r^{(k_R)}(p)$ is calculated relative to $k_R$ and the deviation for $S_r^{(k_e)}(p)$ is calculated relative to $k_e$ . The $S_r^{(k_R)}(p)$ is computed at Ω = 5000 r/min and w = 0.0002 for the k range of $[k_s,k_e]=[20,100]$ and the result presented as a plot of $S_r^{(k_R)}(p)$ against the order p of FDM as shown in Figure 3(a). It is seen that the minimum value of $S_r^{(k_R)}(p)$ occurs at the fourth order meaning that the fourth-order FDM is expected to compete for the most accurate method. In fact, the trend of Figure 3(a) agrees with the work ⁵¹ that accuracy improves in the sequence: first, second, third or fifth and fourth. Additional deduction is that the decrease of accuracy with rising order beyond the third- and fourth-order methods may not necessarily persist being that there is a partial recovery of accuracy before the tenth-order. This extent of clarity would not have been possible from Figure 2. Plots of $S_r^{(k_e)}(p)$ at the process parameter combination Ω = 5000 r/min and w = 0.0002 are generated for different k ranges as shown in Figure 3(b). This figure also suggests best convergence for the fourth-order FDM. A similar plot is generated as shown in Figure 3(c) for the k ranges of $[k_s,k_e]=[40,100]$ and $[k_s,k_e]=[40,120]$ . This figure suggests that the third-order FDM could be the most converged. It must be noted that Figure 3(a) and (b), with common starting value $k_s=20$ , agrees in suggesting that the fourth-order FDM is likely to be the most convergent. Also, both curves in Figure 3(c), with common starting value $k_s=40$ , are consistent in suggesting that the third-order FDM is likely to be most convergent. The benefit of using the $S_r^{(k_R)}(p)$ or $S_r^{(k_e)}(p)$ over $r(k)$ is the clarity with which the former indicates the potentially most accurate methods.

Figure 3.

(a) A plot of $S_r^{(k_R)}(p)$ at the process parameter point Ω = 5000 r/min and w = 0.0002 for the k range of $[k_s,k_e]=[20,100]$ with $k_R=500$ , (b) a plot of $S_r^{(k_e)}(p)$ at the process parameter point Ω = 5000 r/min and w = 0.0002 for the k ranges $[20,k_e]$ and (c) a plot of $S_r^{(k_e)}(p)$ at the process parameter point Ω = 5000 r/min and w = 0.0002 for the k ranges $[40,k_e]$ .

Further calculations of ENRCAP and calculations of computational time (CT) for ENRCAP at additional depth-speed pairs, chosen on the basis of central composite design, indicate that accuracy most likely peaked at the third order, and that CT rises with a rise of order as shown in Figure 4. The natural and coded levels of the depths and speeds used in the simulation are given in Table 2. The coded levels are used to title each element of Figure 4. Among the nine numerical experiments shown in Figure 4, the third-order method is the most accurate in six cases and the first, second and fourth are each the most accurate in a single case.

Figure 4.

Calculations of estimate of summed deviations from the reference spectral radii (ESDRSR) and calculations of computational time (CT) of ESDRSR at depth–speed pairs chosen on the basis of central composite design. The plots are titled with the coded depth–speed pairs.

Table 2.

The natural and coded levels of the depths and speeds used in the simulations of Figure 4.

Factors	Levels
	$-\sqrt{2}$	$-1$	$0$	$1$	$\sqrt{2}$
Ω (r/min)	3863.6	4166.7	5000	6000	6470.7
w (m)	$1.1272\times {10}^{-4}$	$1.3333\times {10}^{-4}$	$2.000\times {10}^{-4}$	$3.000\times {10}^{-4}$	$3.5486\times {10}^{-4}$

The 2DOF case

Equation (1) also governs the 2DOF model but with the matrices and the specific force variation changed as seen in Ozoegwu et al. ⁵¹ Therefore, the generalization does not discriminate between 1DOF and 2DOF systems. To compare the accuracy of various orders of the generalized monodromy for 2DOF model, the $100r_n(k)$ , $\log [r_n(k)]$ and $S_r^{(k_R)}(p)$ (with $k_R=300$ ) are used. In computing and plotting $100r_n(k)$ , $\log [r_n(k)]$ and $S_r^{(k_R)}(p)$ , the discretization integer k is varied in the range of $[k_s,k_e]=[30,100]$ . Three different milling systems are used for simulations so that global conclusions can be drawn. The milling systems include both symmetric and non-symmetric 2DOF tools as shown in Table 3. System 3 in Table 3 is drawn from Tekeli and Budak. ⁶³ The results of convergence and computational time analyses are given in Figures 5–7. The ESDRSR, as plotted in Figure 5, indicates that the third-order method is the most accurate for System 1. Based on ESDRSR, all the orders except the tenth-order seem to be of same accuracy for System 2 as shown in Figure 6. All the methods from third- to tenth-order seem to have same accuracy for System 3 as shown in Figure 7. The most important result is that the methods are very accurate for Systems 2 and 3, reducing percentage error below 0.01% within the plot range $[k_s,k_e]=[30,100]$ .

Table 3.

The parameters used in all numerical computations for 2DOF systems.

Parameter	System 1	System 2	System 3
	Bayly et al. 33	Bobrenkov et al. 36	Tekeli and Budak 63
${\omega }_{\mathrm{nx}}$	$5793\mathrm{rad}{\mathrm{s}}^{-1}$	920.49 rad s−1	3769.9 rad s−1
${\omega }_{\mathrm{ny}}$	$5793\mathrm{rad}{\mathrm{s}}^{-1}$	920.49 rad s−1	4146.9 rad s−1
$m_x$	$0.03993\mathrm{kg}$	2.5729 kg	0.3940 kg
$m_y$	$0.03993\mathrm{kg}$	2.5729 kg	0.3256 kg
${\zeta }_x$	$0.011$	0.0032	0.035
${\zeta }_y$	$0.011$	0.0032	0.035
$\gamma$	$1$	1	1
$C_t$	$6\times {10}^8\mathrm{N}{\mathrm{m}}^{-1-\gamma }$	$5.5\times {10}^8\mathrm{N}{\mathrm{m}}^{-1-\gamma }$	$6\times {10}^8\mathrm{N}{\mathrm{m}}^{-1-\gamma }$
X	$0.33333$	0.3636	$1.1667\times {10}^{-4}$
N	$2$	4	3

DOF: degree of freedom

These parameters are identical with the experimentally determined parameters of cited references.

Figure 5.

The System 1 plots of $100r_n(k)$ , $\log [r_n(k)]$ , ESDRSR and CT at the coded depth–speed pairs (a) −1, −1; (b) 1, 1; (c) 0, 0 and (d) 0, 1.4142.

Figure 6.

The System 2 plots of $100r_n(k)$ , $\log [r_n(k)]$ , ESDRSR and CT at the coded depth–speed pairs (a) −1, −1; (b) 1, 1; (c) 0, 0 and (d) 0, 1.4142.

Figure 7.

The System 3 plots of $100r_n(k)$ , $\log [r_n(k)]$ , ESDRSR and CT at the coded depth–speed pairs (a) −1, −1; (b) 1, 1; (c) 0,0 and (d) 0, 1.4142.

By comparing the exact Taylor series solution of milling process and the approximate solution as given in equation (15) over the first discrete interval $[0,\Delta t]$ , an estimate of the local discretization error of the general-order FDM can be determined. Determination of the local discretization error of the general-order FDM in this manner is extremely difficult; thus, such analysis is done on order-by-order basis to give a generalized error of form ${\mathbf{E}}_p=\mathcal{O}(\Delta t{)}^{p+1}$ . This means that a FDM of order p is expected to have local discretization error that is proportional to the $(p+1)\mathrm{th}$ power of the discrete time. This means that accuracy should increase progressively with order but that is not the case as seen in the numerical convergence analysis. This discrepancy could be due to the Runge phenomenon. Insperger et al. ⁶⁴ explained the problem as due to adoption of first-order Magus series for $\mathbf{B}(t)$ rather the higher order series. Recent works^65–67 have shown that higher order interpolation of the delayed state $\mathbf{y}(s-\tau )$ rather than the matrix $\mathbf{B}(t)$ leads to improved accuracy.

The condition for asymptotic stability is that all the $(2k+2)d$ eigen-values of the monodromy matrix ${\mathbf{\psi }}^{(p)}$ , where d is the degree of milling model, have modulus less than unity. The stability lobes are curves of critical process parameter points at which spectral radii of ${\mathbf{ѱ}}^{(p)}$ are of unit magnitude. The 1DOF milling process with parameters in Table 1 and the 2DOF milling process (System 1) with parameters in Table 3 are considered for illustration. The above results suggest that the FDM of third order is the most accurate for stability analyses of the 1DOF and 2DOF systems. The stability analysis of the 1DOF system using the FDM of third- order is already conducted ⁵¹ while the stability analysis of the 2DOF system using the FDM of third-order is also already conducted. ⁶²

Conclusion

This work generalizes the FDM such that the monodromy matrices are presented as exact function of the order of state approximation. In other words, the usual derivation of monodromy matrices of various FDMs on order-by-order basis (which entails heavy analytical involvement that rapidly gets heavier with the rise in the order of approximation) is bypassed. This is the highest possible level of generalization that seems to be the first of its kind among the time-domain methods. The known generalizations are limited to the interpolating/approximating polynomial of milling state. A potential application is complete computerization of stability analysis such that order of approximation is input with the stability diagram as output rather than the usual time-consuming approach of programming on order-by-order basis. This generalization has made it convenient in this work to study the stability of milling process with the FDMs of first- to tenth-order. By considering the actual geometric nature of characteristic multipliers as complex numbers, more reliable methods for rate of convergence analysis are suggested and utilized in consolidating the known result that best accuracy of the FDM lies with the third- and fourth-order methods. It is seen from numerical convergence analyses that, although accuracy decreases with rising order beyond the third- and fourth-order methods, the trend of decreasing accuracy did not persist with rise of order. The generalized method proved to be very accurate in the stability analysis of 2DOF milling process; less than 0.01% error was attained within $k<100$ in majority of the studied 2DOF systems. It is hoped that the generalization will be extended to the more accurate cases of nonlinear interpolation of the delayed state.