Issue 
Manufacturing Rev.
Volume 1, 2014



Article Number  14  
Number of page(s)  13  
DOI  https://doi.org/10.1051/mfreview/2014014  
Published online  25 September 2014 
Short Note
Reducing computational requirement of stability analysis of milling by partial averaging
Department of Mechanical Engineering, Nnamdi Azikiwe University, PMB 5025, Awka, Nigeria
^{*} email: chigbogug@yahoo.com
Received:
25
July
2014
Accepted:
27
August
2014
The trapezoidal rule and Taylor theorem are used to establish a novel reformulation technique called Partial Averaging (PA). PA is seen to reduce computational time (CT) of computer central processing unit (CPU) without loss of accuracy. Two types of PA are proposed, one involves use of trapezoidal rule to give a map called Full Trapezoidal Rule Map with Partial Averaging (FTRMPA) while the other is achieved through exact analysis to give a map called Partial Trapezoidal Rule Map with Partial Averaging (PTRMPA). The latter is demonstrated to be more accurate and more time conserving in stability analysis of milling process (both one and two degree of freedom systems were considered) and delayed damped Mathieu equation. With all other things being equal the results of stability analysis of fullyimmersed milling process using FTRMPA and PTRMPA are seen to be identical with those of fulldiscretization by Ding et al. [1] further validating the presented reformulation. PA is applied on the Second Order Least Squares Approximated Fulldiscretization method of [18] to illustrate its usefulness for reducing the cumbersome analysis and CT of the fulldiscretization method when analysis gets more complicated by higher order interpolation/approximation theory.
Key words: Chatter / Milling process / Delayed damped Mathieu equation / Stability / Simulation / Delay differential equation
© O.C. Godwin & O. Sam, Published by EDP Sciences, 2014
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
k_{x}: Stiffness in the feed direction, N/m
k_{y}: Stiffness in the feednormal direction, N/m
m, m_{x}: Modal mass in the feed direction, kg
m_{y}: Modal in the feednormal direction, kg
c_{x}: Damping coefficient in the feed direction, Ns/m
c_{y}: Damping coefficient in the feednormal direction, Ns/m
ω_{nx}: Tool natural frequency in feed direction, rad/s
ω_{ny}: Tool natural frequency in feednormal direction, rad/s
ξ_{x}: Tool damping ratio in feed direction, rad/s
ξ_{y}: Tool damping ratio in feednormal direction, rad/s
x(t): Tool total motion in feed direction, m
y(t): Tool total motion in feednormal direction, m
v: Feed speed in feed direction, m/s
f_{a}: Dynamic feed per tooth in feed direction, m
F_{norm, j}(t): Normal cutting force on jth tooth, N
F_{tan, j}(t): Tangential cutting force on jth tooth, N
X: The ratio ${F}_{\mathrm{norm},j}\left(t\right)/{F}_{\mathrm{tan},j}\left(t\right)$
F_{x}(t): Tool cutting force in feed direction, N
F_{y}(t): Tool cutting force in feednormal direction, N
z(t): Regenerative motion in feed direction, m
N: Number of teeth on the milling tool
θ_{j}(t): Angular displacement of jth tooth, rad
ρ = B/D: Radial immersion where
γ: Feed exponent in cutting force law
C_{tan}: Tangential cutting coefficients, Nm^{−1−γ}
1. Introduction
Regenerative chatter has been experimentally and mathematically modelled for over a century to reduce losses caused to the manufacturing industry to the barest minimum but the war is still far from won. Quintana and Ciurana [2] cited the report by Le Lan et al. [3] that Renault S.A.S., a typical automotive industry that produces around 3 million engines a year estimated the cost due to chatter on a cylinder block of a Renault 2.0 lDCi to be around 0.35€ per piece. This scale of loss cannot be overlooked. Mathematical modelling and stability analysis of regenerative chatter comes at a very heavy computational cost. This computational cost can be viewed from the standpoint of amount of analysis, computer resources and computational time needed. Providing new methods of analysis and reformulating the existing ones in a way that reduces computational cost is a step in the right direction. Among the early works that provided adequate mathematical description of the milling process considering multiple cutting edges, interrupted nature of cutting and changing direction of cutting forces is the work [19] which provided linear delay differential equation model with timevarying (periodic) coefficients. After about two decades Minis et al. [20] used Floquet’s theorem and the Fourier series for the formulation of a system for milling stability which they numerically solved using the Nyquist criterion. Shortly afterwards a new method of analysis which involves retaining only the first term in the Fourier series expansion of periodic coefficients appeared in the work [21]. This method of Zeroth Order Approximation was later used in experimental and numerical verification of fast computed stability lobe diagrams [22, 23]. The semidiscretization method [4–7] appeared and enjoyed wide acceptability because of its incisive applicability to many different kinds of delayed dynamical systems. For example, the semidiscretization method revealed some experimentally verified features of milling stability like flip chatter which are not revealed by method of Zeroth Order Approximation but revealed by the method of MultiFrequency Solution [24]. The method of MultiFrequency Solution is an extension of the Zeroth Order Approximation achieved by inclusion of higher order harmonics of the periodic coefficients. Methods to reduce CT and inaccuracy in semidiscretization were subsequently proposed [8]. The semidiscretization method basically involves discretization of the delayed term while leaving the undelayed terms undiscretized and approximating the periodic coefficient matrices as piecewise constant functions. Ding et al. [1] subsequently extended the discretization scheme to the periodic coefficient matrices and the associated state terms to develop the socalled fulldiscretization method. The fulldiscretization method turned out to be slightly less accurate than semidiscretization method of same order [9, 10] but it is as incisive in revealing various stability features of milling. The breakthrough with the advent of fulldiscretization method is that considerable amount of CT has been saved. Quo et al. [11] have improved the accuracy of the fulldiscretization method by developing a map based on the third order interpolation theory. An ingenious application of the Hermite interpolation theory by Liu et al. [12] has led to a fulldiscretization method with further improvements in accuracy and CT. Least squares approximated first and second order fulldiscretization methods have been introduced in [18] to save considerable amount of CT of the original fulldiscretization methods in the works [1, 10]. The basic interest in this work is to introduce a simplification called Partial Averaging (PA) that further reduces computational cost from the perspectives of amount of analysis involved and CT without loss of accuracy. The method PA presented in this work is potentially useful (as illustrated for the Second Order Least Squares Approximated Fulldiscretization method of [18]) for reducing the cumbersome analysis and CT of the fulldiscretization method especially when analysis gets more complicated by the use of higher order interpolation/approximation theory.
2. The basis for reformulation
Consider the application of the trapezoidal rule to the definite integral:$${\mathit{q}}_{1}={\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{\mathit{A}\left({t}_{i+1}s\right)}\mathit{f}\left(s\right)\mathrm{d}s,$$(1)to give$${\mathit{q}}_{1}=\frac{\u2206t}{2}\left({\mathrm{e}}^{\mathit{A}\u2206t}{\mathit{f}}_{i}+{\mathit{f}}_{i+1}\right).$$(2)
Without loss of generality f(s) could be viewed as an ndimensional state vector and A as an n × n constant matrix. The quantity ∆t = t_{i+1} − t_{i} is the interval of integration and ${\mathit{f}}_{i}=\mathit{f}\left({t}_{i}\right)$. Also consider the application of the trapezoidal rule to the definite integral:$${\mathit{q}}_{2}=\frac{1}{\u2206t}{\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{\mathit{A}\left({t}_{i+1}s\right)}\mathrm{d}s{\int}_{{t}_{i}}^{{t}_{i+1}}\mathit{f}\left(s\right)\mathrm{d}s,$$(3)to give$${\mathit{q}}_{2}=\frac{\u2206t}{4}\left({\mathrm{e}}^{\mathit{A}\u2206t}+\mathbf{I}\right)\left({\mathit{f}}_{i}+{\mathit{f}}_{i+1}\right).$$(4)
Equation (3) should be seen to mean definite integral of f(s) on interval $\left[{t}_{i},{t}_{i+1}\right]$ weighted with the mean of ${e}^{\mathit{A}\left({t}_{i+1}s\right)}$ and that constitutes the socalled technique of Partial Averaging (PA) of integration scheme. Making use of the matrix exponential expansion ${\mathrm{e}}^{\mathit{A}\u2206t}=\mathbf{I}+\sum _{p=1}^{\mathit{\infty}}\left({\u2206t}^{p}{\mathbf{A}}^{p}\right)/p!$ in equations (2) and (4) respectively gives:$${\mathit{q}}_{1}=\frac{\u2206t}{2}\left[\left(\sum _{p=1}^{\infty}\frac{{\u2206t}^{p}{\mathbf{A}}^{p}}{p!}\right){\mathit{f}}_{i}+{\mathit{f}}_{i}+{\mathit{f}}_{i+1}\right],$$(5) $${\mathit{q}}_{2}=\frac{\u2206t}{4}\left(\sum _{p=1}^{\mathit{\infty}}\frac{{\u2206t}^{p}{\mathbf{A}}^{p}}{p!}+2\mathbf{I}\right)\left({\mathit{f}}_{i}+{\mathit{f}}_{i+1}\right).$$(6)
The vector q_{2} is multiplied out and rewritten to become:$${\mathit{q}}_{2}=\frac{\u2206t}{2}\left[\left(\sum _{p=1}^{\mathit{\infty}}\frac{{\u2206t}^{p}{\mathbf{A}}^{p}}{p!}\right){\mathit{f}}_{i}+{\mathit{f}}_{i}+{\mathit{f}}_{i+1}+\frac{1}{2}\left(\sum _{p=1}^{\mathit{\infty}}\frac{{\u2206t}^{p}{\mathbf{A}}^{p}}{p!}\right)\left({\mathit{f}}_{i+1}{\mathit{f}}_{i}\right)\right].$$(7)
Equation (5) is substituted in equation (7) to give:$${\mathit{q}}_{2}{\mathit{q}}_{1}=\frac{1}{4}\left[\left(\sum _{p=1}^{\mathit{\infty}}\frac{{\u2206t}^{p+1}{\mathbf{A}}^{p}}{p!}\right)\left({\mathit{f}}_{i+1}{\mathit{f}}_{i}\right)\right].$$(8)
Equation (8) means that the error of PA is of the order of ∆t^{2} thus negligible when ∆t becomes very small due to shrink of the interval $\left[{t}_{i},{t}_{i+1}\right]$ and due to very small size of the difference f_{i+1} − f_{i}. This is symbolically written as ${\mathit{q}}_{2}{\mathit{q}}_{1}=\mathcal{O}{\u2206t}^{2}$ to mean that the error of PA over the interval $\left[{t}_{i},{t}_{i+1}\right]$ is proportional to the square of the size of the integration step ∆t. The next thing is to numerically verify this result for the milling process.
3. Application of PA in chatter stability analysis of milling process
The dynamical model for regenerative milling process is [1, 9–13, 18]:$${\stackrel{\u0307}{\mathit{\xi}}}^{\left(\mathrm{d}\right)}\left(t\right)={\mathbf{A}}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\right)+{\mathbf{B}}^{\left(\mathrm{d}\right)}\left(t\right){\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\right){\mathbf{B}}^{\left(\mathrm{d}\right)}\left(t\right){\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\tau \right),$$(9)where d = 1 or 2 denotes the degree of freedom of the milling process. If the tool is compliant only in the feed direction then d = 1 but if the tool is simultaneously compliant in the feed and feednormal directions then d = 2. These cases are illustrated in Figure 1.
Figure 1. Mechanical model of milling tool with indications of dynamic (F_{tan,j}(t), F_{norm,j}(t), F_{x}(t), θ_{j}(t)), process (v, B, Ω) and constant (k_{x}, k_{y}, c_{x}, c_{y}, m etc.) parameters. (a) 1DOF model and (b) 2DOF model. 
The vectors ${\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\right)$ and ${\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\tau \right)$ respectively mean the state and the delayed state of the system at time t. The constant matrix ${\mathbf{A}}^{\left(\mathrm{d}\right)}$ contains the timeinvariant parameters of the system while the periodic coefficient matrix ${\mathbf{B}}^{\left(\mathrm{d}\right)}\left(t\right)$ captures the periodicity of cutting force of the unperturbed milling process. Dividing the discrete delay τ of the system into k equal discrete time intervals $\left[{t}_{i},{t}_{i+1}\right]$ where i = 0, 1, 2, … (k − 1) and ${t}_{i}=i\frac{\tau}{k}=i\u2206t$ and approximating equation (9) in each discrete interval gives the local model:$${\stackrel{\u0307}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(t\right)={\mathbf{A}}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\right)+{\stackrel{\u0303}{\mathbf{B}}}^{\left(\mathrm{d}\right)}\left(t\right){\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(t\right){\stackrel{\u0303}{\mathbf{B}}}^{\left(\mathrm{d}\right)}\left(t\right){\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(t\tau \right),t\u03f5\left[{t}_{i},{t}_{i+1}\right].$$(10)
The basic idea of discretization employed in arriving at equation (10) is that ${\mathbf{B}}^{\left(\mathrm{d}\right)}\left(t\right){\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\right)$ and ${\mathbf{B}}^{\left(\mathrm{d}\right)}\left(t\right){\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\tau \right)$ are discretized while ${\stackrel{\u0307}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(t\right)$ and ${\mathbf{A}}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}^{\left(\mathrm{d}\right)}\left(t\right)$ are left undiscretized. This idea is based on the fact that the exact integral of the undiscretized terms can be determined while those of the discretized terms cannot. The solution of equation (10) in the interval $\left[{t}_{i},{t}_{i+1}\right]$ reads:$${\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}={\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+{\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\left({t}_{i+1}s\right)}\left[{\stackrel{\u0303}{\mathbf{B}}}^{\left(\mathrm{d}\right)}\left(s\right){\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(s\right){\stackrel{\u0303}{\mathbf{B}}}^{\left(\mathrm{d}\right)}\left(s\right){\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(s\tau \right)\right]\mathrm{d}s.$$(11)
The trapezoidal rule is used in solving the Duhamel term such that equation (11) becomes:$${\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}={\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+\frac{\u2206t}{2}\left({\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}\right)\frac{\u2206t}{2}\left({\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{ik}^{\left(\mathrm{d}\right)}+{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i+1k}^{\left(\mathrm{d}\right)}\right).$$(12)
It should be noted that the first published work on application of NewtonCortes method in chatter stability analysis of milling process is that of Ding et al. [13]. Equation (12) is rearranged and put in matrix form to read:$$\left\{\begin{array}{c}{\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}\\ {\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}\\ {\mathbf{\xi}}_{i1}^{\left(\mathrm{d}\right)}\\ \vdots \\ {\mathbf{\xi}}_{i+1k}^{\left(\mathrm{d}\right)}\end{array}\right\}=\left[\begin{array}{cccccc}{\mathbf{M}}_{11}^{i}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{M}}_{1k}^{i}& {\mathbf{M}}_{1,k+1}^{i}\\ \mathbf{I}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right]\left\{\begin{array}{c}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}\\ {\mathbf{\xi}}_{i1}^{\left(\mathrm{d}\right)}\\ {\mathbf{\xi}}_{i2}^{\left(\mathrm{d}\right)}\\ \vdots \\ {\mathbf{\xi}}_{ik}^{\left(\mathrm{d}\right)}\end{array}\right\},$$(13)where:$${\mathbf{G}}_{i+1}=\mathbf{I}\frac{\u2206t}{2}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)},$$(14) $$\begin{array}{c}{\mathbf{M}}_{11}^{i}={{\mathbf{G}}_{i+1}}^{1}{\mathbf{F}}_{0}\left(\mathbf{I}+\frac{\u2206t}{2}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}\right)\\ {\mathbf{M}}_{1k}^{i}=\frac{\u2206t}{2}{{\mathbf{G}}_{i+1}}^{1}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\\ {\mathbf{M}}_{1,k+1}^{i}=\frac{\u2206t}{2}{{\mathbf{G}}_{i+1}}^{1}{\mathbf{F}}_{0}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}\end{array}.$$(15)
The map for the system covering the whole discrete delay is constructed to become:$${\mathbf{x}}_{k+1}={\mathbf{M}}_{k1}\dots {\mathbf{M}}_{0}{\mathbf{x}}_{0}.$$(16)
This map is called Full Trapezoidal Rule Map (FTRM) in this work because equation (11) is completely handled using the Trapezoidal Rule. The Floquet transition matrix reads:$$\mathbf{\psi}={\mathbf{M}}_{k1}\dots {\mathbf{M}}_{0}.$$(17)
The concept of PA means that equation (11) can be represented as:$${\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}={\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+\frac{1}{\u2206t}{\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\left({t}_{i+1}s\right)}\mathrm{d}s{\int}_{{t}_{i}}^{{t}_{i+1}}\left[{\stackrel{\u0303}{\mathbf{B}}}^{\left(\mathrm{d}\right)}\left(s\right){\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(s\right){\stackrel{\u0303}{\mathbf{B}}}^{\left(\mathrm{d}\right)}\left(s\right){\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(s\tau \right)\right]\mathrm{d}s.$$(18)
Application of the trapezoidal rule to each of the two integrals of the second term of equation (18) gives:$${\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}={\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+\frac{\u2206t}{4}\left({\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}+\mathbf{I}\right)\left({\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}\right)\frac{\u2206t}{4}\left({\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}+\mathbf{I}\right)\left({\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{ik}^{\left(\mathrm{d}\right)}+{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i+1k}^{\left(\mathrm{d}\right)}\right).$$(19)
The monodromy matrix is constructed as M_{k−1}… M_{0} where:$${\mathbf{M}}_{i}=\left[\begin{array}{cccccc}{\mathbf{M}}_{11}^{i}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{M}}_{1k}^{i}& {\mathbf{M}}_{1,k+1}^{i}\\ \mathbf{I}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right],$$(20) $$\begin{array}{c}{\mathbf{M}}_{11}^{i}={\mathbf{P}}_{i+1}\left({\mathbf{F}}_{0}+\frac{\u2206t}{4}\left({\mathbf{F}}_{0}+\mathbf{I}\right){\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}\right)\\ {\mathbf{M}}_{1k}^{i}=\frac{\u2206t}{4}{\mathbf{P}}_{i+1}\left({\mathbf{F}}_{0}+\mathbf{I}\right){\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\\ {\mathbf{M}}_{1,k+1}^{i}=\frac{\u2206t}{4}{\mathbf{P}}_{i+1}\left({\mathbf{F}}_{0}+\mathbf{I}\right){\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}\end{array},$$(21) $${\mathbf{P}}_{i+1}={\left(\mathbf{I}\frac{\u2206t}{4}\left({\mathbf{F}}_{0}+\mathbf{I}\right){\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\right)}^{1}.$$(22)
Let this map be called Full Trapezoidal Rule Map with Partial Averaging (FTRMPA) in this work because equation (11) is completely handled using the Trapezoidal Rule after PA was introduced.
Application of the trapezoidal rule in equation (18) after exact evaluation of the first integral as:$$\mathit{E}=\frac{1}{\u2206t}{\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\left({t}_{i+1}s\right)}\mathrm{d}s=\frac{1}{\u2206t}{\int}_{0}^{\u2206t}{e}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\left(\u2206ts\right)}\mathrm{d}s=\frac{1}{\u2206t}\left({\mathbf{F}}_{0}\mathbf{I}\right){\left({\mathbf{A}}^{\left(\mathrm{d}\right)}\right)}^{1},{\mathbf{F}}_{0}={\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t},$$(23)gives$${\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}={\mathrm{e}}^{{\mathit{A}}^{\left(\mathrm{d}\right)}\u2206t}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+\frac{\u2206t}{2}\mathit{E}\left({\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}\right)\frac{\u2206t}{2}\mathit{E}\left({\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{ik}^{\left(\mathrm{d}\right)}+{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}{\mathbf{\xi}}_{i+1k}^{\left(\mathrm{d}\right)}\right).$$(24)
This leads to the construction of the monodromy matrix as M_{k−1}… M_{0} where:$${\mathbf{M}}_{i}=\left[\begin{array}{cccccc}{\mathbf{M}}_{11}^{i}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{M}}_{1k}^{i}& {\mathbf{M}}_{1,k+1}^{i}\\ \mathbf{I}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right],$$(25) $$\begin{array}{c}{\mathbf{M}}_{11}^{i}={\mathbf{P}}_{i+1}\left({\mathbf{F}}_{0}+\frac{\u2206t}{2}\mathit{E}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}\right)\\ {\mathbf{M}}_{1k}^{i}=\frac{\u2206t}{2}{\mathbf{P}}_{i+1}\mathit{E}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\\ {\mathbf{M}}_{1,k+1}^{i}=\frac{\u2206t}{2}{\mathbf{P}}_{i+1}\mathit{E}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}\end{array},$$(26) $${\mathbf{P}}_{i+1}={\left(\mathbf{I}\frac{\u2206t}{2}\mathit{E}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\right)}^{1}.$$(27)
Partial application of the trapezoidal rule is implemented in this case since the first integral in equation (18) is replaced with its exact value as given in equation (23) while the trapezoidal rule is used in handling the second integral. The arising map is then called the Partial Trapezoidal Rule Map with Partial Averaging (PTRMPA) because equation (18) is handled with the Trapezoidal Rule after exact Partial Averaging was carried out. It will be shown in what follows that the concept of PA is applicable to chatter stability analysis by showing in what follows that FTRM, FTRMPA and PTRMPA exhibit very similar rate of convergence and stability lobes.
4. Numerical simulation of milling stability
A numerical simulation requires complete description of equation (9). For completeness and independence of this work a brief derivation of equation (9) for the 1DOF system given in Figure 1a is as follows; the governing equation of motion of the 1DOF tool in terms of the total motion x(t) reads:$$m\stackrel{\u0308}{x}\left(t\right)+{c}_{x}\left[\stackrel{\u0307}{x}\left(t\right)v\right]+{k}_{x}\left[x\left(t\right)\mathrm{vt}\right]+{F}_{x}\left(t\right)=0.$$(28)
The xresultant of cutting force stemming from the chipforming interaction between the tool and workpiece is seen in Figure 1a to be:$${F}_{x}\left(t\right)=\sum _{j=1}^{N}{g}_{j}\left(t\right){F}_{\mathrm{tan},j}\left(t\right)\left[X\mathrm{sin}{\theta}_{j}\left(t\right)+\mathrm{cos}{\theta}_{j}\left(t\right)\right],$$(29)where N is the number of teeth on the milling tool indexed with the values j = 1, 2, 3… N and $X={F}_{\mathrm{norm},j}\left(t\right)/{F}_{\mathrm{tan},j}\left(t\right)$ is the ratio of normal to tangential cutting force on the jth tooth.
The instantaneous angular position of jth tooth ${\theta}_{j}\left(t\right)$ is given by:$${\theta}_{j}\left(t\right)=\left(\frac{\pi \mathrm{\Omega}}{30}\right)t+\left(j1\right)\frac{2\pi}{N}.$$(30)
The screen function ${g}_{j}\left(t\right)$ has values of either unity or nullity depending on whether the tool is cutting or not. Given start and end angles of cut θ_{s} and θ_{e}, respectively, ${g}_{j}\left(t\right)$ becomes [14]:$${g}_{j}\left(t\right)=\{\begin{array}{c}\begin{array}{cc}1& \mathrm{if}\hspace{1em}{\theta}_{s}<{\theta}_{j}\left(t\right){\theta}_{e}\end{array}\\ \begin{array}{cc}0& \mathrm{otherwise}\end{array}\end{array}.$$(31)
The precise mathematical form for ${g}_{j}\left(t\right)$ is given as [4]:$${g}_{j}\left(t\right)=\frac{1}{2}\left(1+\mathrm{sgn}\left\{\mathrm{sin}\left[{\theta}_{j}\left(t\right){\mathrm{tan}}^{1}\mathcal{P}\right]\mathrm{sin}\left[{\theta}_{s}{\mathrm{tan}}^{1}\mathcal{P}\right]\right\}\right),$$(32)where $\mathcal{P}=\left(\mathrm{sin}{\theta}_{s}\mathrm{sin}{\theta}_{e}\right)/\left(\mathrm{cos}{\theta}_{s}\mathrm{cos}{\theta}_{e}\right)$. The screen function ${g}_{j}\left(t\right)$ has been formulated in terms of radial immersion ρ = B/D and ${\theta}_{j}\left(t\right)$ for upmilling and downmilling modes respectively as [25]:$${g}_{j}\left(t\right)=\frac{1}{2}\left\{1+\mathrm{sgn}\left[\mathrm{sin}\left({\theta}_{j}\left(t\right)\mathrm{arctan}\left\{\frac{1}{2\rho}\mathrm{sin}\left[\mathrm{arccos}\left(12\rho \right)\right]\right\}\right)+\mathrm{sin}\left(\mathrm{arctan}\left\{\frac{1}{2\rho}\mathrm{sin}\left[\mathrm{arccos}\left(12\rho \right)\right]\right\}\right)\right]\right\},$$(33a) $${g}_{j}\left(t\right)=\frac{1}{2}\left\{1+\mathrm{sgn}\left[\mathrm{sin}\left({\theta}_{j}\left(t\right)\mathrm{arctan}\left\{\frac{1}{2\rho}\mathrm{sin}\left[\mathrm{arccos}\left(2\rho 1\right)\right]\right\}\right)\mathrm{sin}\left(\mathrm{arccos}\left(2\rho 1\right)\mathrm{arctan}\left\{\frac{1}{2\rho}\mathrm{sin}\left[\mathrm{arccos}\left(12\rho \right)\right]\right\}\right)\right]\right\}.$$(33b)
The tangential cutting force on jth tooth is generally given by the nonlinear law [4]:$${F}_{\mathrm{tan},j}\left(t\right)={C}_{\mathrm{tan}}w{\left[{f}_{a}\mathrm{sin}{\theta}_{j}\left(t\right)\right]}^{\gamma},$$(34)where w is the axial depth of cut, C_{tan} is the tangential cutting coefficient associated with the workpiece material and tool geometry, f_{a} is the actual feed per tooth (i.e. feed within single discrete delay τ = 60/(NΩ)) given as:$${f}_{a}=x\left(t\right)x\left(t\tau \right),$$(35)and γ is a positive exponent. Equations (29), (34) and (35) taken together give:$${F}_{x}\left(t\right)=\mathrm{wq}\left(t\right){\left[x\left(t\right)x\left(t\tau \right)\right]}^{\gamma},$$(36)where q(t) is a τ( = 60/NΩ) periodic function given by:$$q\left(t\right)=\sum _{j=1}^{N}{g}_{j}\left(t\right){C}_{\mathrm{tan}}{\mathrm{sin}}^{\gamma}{\theta}_{j}\left(t\right)\left[X\mathrm{sin}{\theta}_{j}\left(t\right)+\mathrm{cos}{\theta}_{j}\left(t\right)\right].$$(37)
The regenerative motion z(t) of the tool is given by:$$z\left(t\right)=x\left(t\right)\mathrm{vt}{x}_{\mathrm{t}}\left(t\right),$$(38)where x_{t}(t) is the τperiodic response due to periodic force of toolworkpiece interaction. Then making use of equations (38) and (36) in equation (1) when z(t) = 0 gives:$$m{\stackrel{\u0308}{x}}_{\mathrm{t}}\left(t\right)+{c}_{x}{\stackrel{\u0307}{x}}_{\mathrm{t}}\left(t\right)+{k}_{x}{x}_{\mathrm{t}}\left(t\right)=\mathrm{wq}\left(t\right){\left(\mathrm{v\tau}\right)}^{\gamma}.$$(39)
Inserting equations (39) and (38) in equation (28) gives:$$m\stackrel{\u0308}{z}\left(t\right)+{c}_{x}\stackrel{\u0307}{z}\left(t\right)+{k}_{x}z\left(t\right)=\mathrm{wq}\left(t\right){\left(\mathrm{v\tau}\right)}^{\gamma}\mathrm{wq}\left(t\right){\left\{\mathrm{v\tau}+\left[z\left(t\right)z\left(t\tau \right)\right]\right\}}^{\gamma}.$$(40)
This is linearized to become:$$m\stackrel{\u0308}{z}\left(t\right)+{c}_{x}\stackrel{\u0307}{z}\left(t\right)+{k}_{x}z\left(t\right)=wh\left(t\right)\left[z\left(t\right)z\left(t\tau \right)\right],$$(41)where:$$h\left(t\right)=\gamma {\left(\mathrm{v\tau}\right)}^{\gamma 1}q\left(t\right),$$(42)is the specific periodic cutting force variation. The modal form of governing model chatter of milling process becomes:$$\stackrel{\u0308}{z}\left(t\right)+2{\xi}_{x}{\omega}_{\mathrm{nx}}\stackrel{\u0307}{z}\left(t\right)+\left({\omega}_{\mathrm{nx}}^{2}+\frac{wh\left(t\right)}{m}\right)z\left(t\right)=\frac{wh\left(t\right)}{m}z\left(t\tau \right).$$(43)
The natural frequency and damping ratio of the tool are given as ${\omega}_{\mathrm{nx}}=\sqrt{{k}_{x}/m}$ and ${\xi}_{x}={c}_{x}/2\sqrt{m{k}_{x}}$. Equation (43) is easily put in state space form of equation (1) such that:$${\mathbf{A}}^{\left(1\right)}=\left[\begin{array}{cc}0& 1\\ {\omega}_{\mathrm{nx}}^{2}& 2{\xi}_{x}{\omega}_{\mathrm{nx}}\end{array}\right],\hspace{1em}{\mathbf{B}}^{\left(1\right)}\left(t\right)=\left[\begin{array}{cc}0& 0\\ \frac{wh\left(t\right)}{m}& 0\end{array}\right].$$(44)
In the two degree of freedom case (d = 2) the nonlinear force is used in a similar manner taking into account tool compliance in the feed and feednormal directions to give:$${\mathbf{A}}^{\left(2\right)}=\left[\begin{array}{cccc}0& 1& 0& 0\\ {\omega}_{\mathrm{nx}}^{2}& 2{\xi}_{x}{\omega}_{\mathrm{nx}}& 0& 0\\ 0& 0& 0& 1\\ 0& 0& {\omega}_{\mathrm{ny}}^{2}& 2{\xi}_{y}{\omega}_{\mathrm{ny}}\end{array}\right],$$(45) $${\mathbf{B}}^{\left(2\right)}\left(t\right)=\left[\begin{array}{cccc}0& 0& 0& 0\\ \frac{w{h}_{\mathrm{xx}}\left(t\right)}{{m}_{x}}& 0& \frac{w{h}_{\mathrm{xy}}\left(t\right)}{{m}_{x}}& 0\\ 0& 0& 0& 0\\ \frac{w{h}_{\mathrm{yx}}\left(t\right)}{{m}_{y}}& 0& \frac{w{h}_{\mathrm{yy}}\left(t\right)}{{m}_{y}}& 0\end{array}\right],$$(46) $${h}_{\mathrm{xx}}\left(t\right)={C}_{\mathrm{tan}}\gamma {\left(\mathrm{v\tau}\right)}^{\gamma 1}\sum _{j=1}^{N}{g}_{j}\left(t\right){\mathrm{sin}}^{\gamma}{\theta}_{j}\left(t\right)\left[X\mathrm{sin}{\theta}_{j}\left(t\right)+\mathrm{cos}{\theta}_{j}\left(t\right)\right],$$(47a) $${h}_{\mathrm{xy}}\left(t\right)={C}_{\mathrm{tan}}\gamma {\left(\mathrm{v\tau}\right)}^{\gamma 1}\sum _{j=1}^{N}{g}_{j}\left(t\right){\mathrm{sin}}^{\gamma 1}{\theta}_{j}\left(t\right)\mathrm{cos}{\theta}_{j}\left(t\right)\left[X\mathrm{sin}{\theta}_{j}\left(t\right)+\mathrm{cos}{\theta}_{j}\left(t\right)\right],$$(47b) $${h}_{\mathrm{yx}}\left(t\right)={C}_{\mathrm{tan}}\gamma {\left(\mathrm{v\tau}\right)}^{\gamma 1}\sum _{j=1}^{N}{g}_{j}\left(t\right){\mathrm{sin}}^{\gamma}{\theta}_{j}\left(t\right)\left[X\mathrm{cos}{\theta}_{j}\left(t\right)\mathrm{sin}{\theta}_{j}\left(t\right)\right],$$(47c) $${h}_{\mathrm{yy}}\left(t\right)={C}_{\mathrm{tan}}\gamma {\left(\mathrm{v\tau}\right)}^{\gamma 1}\sum _{j=1}^{N}{g}_{j}\left(t\right){\mathrm{sin}}^{\gamma 1}{\theta}_{j}\left(t\right)\mathrm{cos}{\theta}_{j}\left(t\right)\left[X\mathrm{cos}{\theta}_{j}\left(t\right)\mathrm{sin}{\theta}_{j}\left(t\right)\right].$$(47d)
The directional modal parameters are designated k_{x}, m_{x}, c_{x}, k_{y}, m_{y} and c_{y}. The directional natural frequencies and damping ratios are given as ${\omega}_{\mathrm{nx}}=\sqrt{{k}_{x}/{m}_{x}}$, ${\omega}_{\mathrm{ny}}=\sqrt{{k}_{y}/{m}_{y}}$, ${\xi}_{x}={c}_{x}/2\sqrt{{m}_{x}{k}_{x}}$ and ${\xi}_{y}={c}_{y}/2\sqrt{{m}_{y}{k}_{y}}$. The parameters used for numerical simulation are contained in Table 1.
4.1. Rate of convergence of milling discrete maps
Both rate of convergence and stability analysis will be used in assessing the accuracy of PA. The rate of convergence of the maps is judged from the asymptotic behaviour of $f\left(k\right)={\left\mu \left(k\right)\right}_{\mathrm{max}}$ with rise in k where f(k) is the absolute value of the maximummagnitude characteristic multiplier of the system at the approximation parameter k. A map is more convergent when f(k) approaches faster an asymptote defined by the exact value of ${\left\mu \left(k\right)\right}_{\mathrm{max}}$ at k = ∞. The asymptotic rate of convergence of fullyimmersed 1DOF milling process is presented in Figure 2. The milling process parameters used in generating Figure 2 are seen in the caption. It should be noted that PA only slightly affects rate of convergence when judged relative to FTRM. It should be noted that Figure 2c even suggests that PTRMPA could be more convergent than FTRM at some parameter points. The asymptotic rate of convergence of fullyimmersed 2DOF milling process is presented in Figure 3 for selected parameter points specified in the caption. Even though rate of convergence of FTRMPA and PTRMPA agrees very closely, it is seen that PTRMPA is consistently slightly more convergent and thus advocated as the preferred choice.
Figure 2. The asymptotic rate of convergence of fullyimmersed 1DOF milling process at the parameter points (a) Ω = 5000 and w = 0.0002 m (stable operation), (b) Ω = 20000 and w = 0.001 m (stable operation) and (c) Ω = 20000 and w = 0.002 m (unstable operation). It is seen that partial averaging does not diminish rate of convergence much. It is also seen that PTRMPA is always slightly more convergent than FTRMPA and thus advocated. 
Figure 3. The asymptotic rate of convergence of fullyimmersed 2DOF milling process at the parameter points (a) Ω = 5000 and w = 0.00001 m (stable operation), (b) Ω = 20,000 and w = 0.00001 m (stable operation) and (c) Ω = 20,000 and w = 0.0001 m (unstable operation). It is seen that partial averaging does not diminish rate of convergence much. It is also seen that PTRMPA is always slightly more convergent than FTRMPA and thus advocated. 
4.2. Computation of milling stability diagrams
Making use of the monodromy matrices of FTRM, FTRMPA and PTRMPA stability boundaries can be tracked under the condition that boundary characteristic multipliers are of highest magnitude and lie on the circumference of a unit circle centred at the origin of the complex plane. These boundary characteristic multipliers result from critical parameter combinations of spindle speed and depth of cut. Characteristic multipliers are computed with k = 40 on a 200 by 100 gridded plane of spindle speed and depth of cut using FTRM, FTRMPA and PTRMPA and the boundary characteristic multipliers connected to separate the plane into stable (upper plane) and unstable (lower) subregions as given in Figure 4 for 1DOF system. Only slight deviation is seen to be caused by PA at low spindle speed. The more striking thing to note is that results of FTRMPA and PTRMPA are almost identical with the result of Fulldiscretization by Ding et al. [1]. This suggests that PA is strongly applicable to the Fulldiscretization as verified in the next section. It should be noted that PTRMPA is less deviated from FTRM than FTRMPA. This means that the suggestion from rate of convergence analysis that PTRMPA is slightly more accurate than FTRMPA is true. The computational time (CT) as seen in the legend of each figure generally increases in the order PTRMPA, FTRMPA and FTRM. It is then concluded that PA is beneficial in terms of reducing computational difficulty and saving CT and that PTRMPA is superior to FTRMPA from both points of view of accuracy and CT.
Figure 4. Stability lobes of (a) fullyimmersed 1DOF milling process using FTRM (CT = 426 s), FTRMPA (CT = 434 s) and PTRMPA (CT = 417 s) placed on same axes to show that accuracy is marginally affected at low spindle speed, (b) halfimmersed 1DOF downmilling process using FTRM (CT = 433 s), FTRMPA (CT = 408 s) (c) PTRMPA (CT = 397 s) placed on same axes to reveal that PA is accurate relative to FTRM, (c) halfimmersed 1DOF upmilling process using FTRM (CT = 439 s), FTRMPA (CT = 404 s), PTRMPA (CT = 395 s) placed on same axes to reveal good accuracy of PA. 
The technique of PA is even more validated by the stability result of 2DOF fullyimmersed milling process since no deviation is seen amongst FTRM, FTRMPA and PTRMPA in Figure 5.
Figure 5. Stability lobes of fullyimmersed 2DOF milling process using FTRM, FTRMPA and PTRMPA placed on same axes to show that the presented two types PA retain accuracy. 
5. Application of PA in simplifying the second order fulldiscretization method
The second order least squares approximated polynomial for the state term ${\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(t\right)$ reads [18]:$${\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(t\right)=\frac{1}{2{\left(\u2206t\right)}^{2}}\left({t}^{2}t\u2206t\right){\mathbf{\xi}}_{i1}^{\left(\mathrm{d}\right)}+\frac{1}{{\left(\u2206t\right)}^{2}}\left[{\left(\u2206t\right)}^{2}{t}^{2}\right]{\mathbf{\xi}}_{i}^{\left(\mathrm{d}\right)}+\frac{1}{2{\left(\u2206t\right)}^{2}}\left({t}^{2}+t\u2206t\right){\mathbf{\xi}}_{i+1}^{\left(\mathrm{d}\right)}.$$(48)
Equation (48) together with the linear approximations:$${\stackrel{\u0303}{\mathbf{\xi}}}^{\left(\mathrm{d}\right)}\left(t\tau \right)=\frac{1}{\u2206t}\left(\u2206tt\right){\mathbf{\xi}}_{ik}^{\left(\mathrm{d}\right)}+\frac{1}{\u2206t}t{\mathbf{\xi}}_{i+1k}^{\left(\mathrm{d}\right)},$$(49) $${\stackrel{\u0303}{\mathbf{B}}}^{\left(\mathrm{d}\right)}\left(t\right)=\frac{1}{\u2206t}\left(\u2206tt\right){\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}+\frac{1}{\u2206t}t{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)},$$(50)is inserted in equation (18) and solved to construct the monodromy operator ψ = M_{k−1,}M_{k−2}… M_{0} such that:$${\mathbf{M}}_{i}=\left[\begin{array}{cccccc}{\mathbf{M}}_{11}^{i}& {\mathbf{M}}_{12}^{i}& \cdots & \mathbf{0}& {\mathbf{M}}_{1k}^{i}& {\mathbf{M}}_{1,k+1}^{i}\\ \mathbf{I}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right],$$(51)where:$${\mathbf{M}}_{11}^{i}={\mathbf{P}}_{i}\left({\mathbf{F}}_{0}+{\mathit{G}}_{22}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}+{\mathit{G}}_{25}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\right),$$(52a) $${\mathbf{M}}_{12}^{i}={\mathbf{P}}_{i}\left({\mathit{G}}_{21}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}+{\mathit{G}}_{24}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\right),$$(52b) $${\mathbf{M}}_{1k}^{i}={\mathbf{P}}_{i}\left({\mathit{G}}_{12}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}+{\mathit{G}}_{13}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\right),$$(52c) $${\mathbf{M}}_{1,k+1}^{i}={\mathbf{P}}_{i}\left({\mathit{G}}_{11}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}+{\mathit{G}}_{12}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\right),$$(52d) $${\mathbf{P}}_{i}={\left[\mathbf{I}{\mathit{G}}_{23}{\mathbf{B}}_{i}^{\left(\mathrm{d}\right)}{\mathit{G}}_{26}{\mathbf{B}}_{i+1}^{\left(\mathrm{d}\right)}\right]}^{1},$$(53) $${\mathit{G}}_{11}=\frac{\u2206t}{3}\mathit{E},\hspace{0.5em}{\mathit{G}}_{12}=\frac{\u2206t}{6}\mathit{E},\hspace{0.5em}{\mathit{G}}_{13}={\mathit{G}}_{11},\hspace{0.5em}{\mathit{G}}_{21}=\frac{\u2206t}{24}\mathit{E},\hspace{0.5em}{\mathit{G}}_{22}=\frac{5\u2206t}{12}\mathit{E},\hspace{0.5em}{\mathit{G}}_{23}=\frac{\u2206t}{8}\mathit{E},\hspace{0.5em}{\mathit{G}}_{24}=\frac{\u2206t}{24}\mathit{E},\hspace{0.5em}{\mathit{G}}_{25}=\frac{\u2206t}{4}\mathit{E},\hspace{0.5em}{\mathit{G}}_{26}=\frac{7\u2206t}{24}\mathit{E}.$$(54)
The generated map is called the Second Order Least Squares Approximated Map with Partial averaging (2nd OLSAMPA). Simplicity provided by PA is obvious when the six G matrices in equation (54) are compared with those of second order least squares approximated fulldiscretization method in [18]. The G matrices in equation (54) are each given in terms of single product of some fraction of linear power of ∆t and E while the G matrices in [18] are each given in terms of sum of products of various nonlinear powers of ∆t and F matrices. Making use of presented map with k = 40 on a 200 by 100 grid produces the stability diagram of fullyimmersed milling as shown in Figure 6a. The original Second Order Least Squares Approximated Fulldiscretization method presented in [18] is also used to generate stability diagram of same parameters as seen in Figure 6b. Agreement seen between the two results highlights the validity of applying PA in the fulldiscretization method in terms of conservation of accuracy and CT savings.
Figure 6. Stability lobes of fullyimmersed 1DOF milling process using (a) Second Order Least Squares Approximated Fulldiscretization method with Partial averaging (CT = 430 s) and (b) the original Second Order Least Squares Approximated Fulldiscretization method presented in [18] (CT = 442 s). Very good agreement between the two verifies reliability in terms of accuracy of using PA in the Fulldiscretization method. 
6. Application of PA in stability analysis of delayed damped Mathieu equation
The delayed damped Mathieu equation is normally given as:$$\stackrel{\u0308}{x}\left(t\right)+\kappa \stackrel{\u0307}{x}\left(t\right)+\left[\delta +\u03f5\mathrm{cos}\left(\omega t\right)\right]x\left(t\right)=\mathrm{bx}\left(t\tau \right).$$(55)
This is put in state space as:$$\stackrel{\u0307}{\mathit{x}}\left(t\right)=\mathbf{A}\mathit{x}\left(t\right)+\mathbf{B}\left(t\right)\mathit{x}\left(t\right)\mathbf{C}\mathit{x}\left(t\tau \right),$$(56)where $\mathit{x}\left(t\right)=\left\{\begin{array}{c}x\left(t\right)\\ \stackrel{\u0307}{x}\left(t\right)\end{array}\right\}$, $\mathit{x}\left(t\tau \right)=\left\{\begin{array}{c}x\left(t\tau \right)\\ \stackrel{\u0307}{x}\left(t\tau \right)\end{array}\right\}$, $\mathbf{A}=\left[\begin{array}{cc}0& 1\\ \delta & \kappa \end{array}\right]$, $\mathbf{B}\left(t\right)=\left[\begin{array}{cc}0& 0\\ \u03f5\mathrm{cos}\left(\omega t\right)& 0\end{array}\right]$ and $\mathbf{C}=\left[\begin{array}{cc}0& 0\\ b& 0\end{array}\right]$. The time delay is given as τ = 2π while the period of the matrix B(t) is given as T = 2π/ω. The weight of the delayed feedback state is captured in the gain parameter b. The symbol κ is the damping coefficient of the system. Equation (55) is divided into k equal discrete time intervals $\left[{t}_{i},{t}_{i+1}\right]$ as done in the case of milling process and integrated between the limits t_{i} and t_{i+1} to give:$${\mathbf{x}}_{i+1}={\mathrm{e}}^{\mathit{A}\u2206t}{\mathbf{x}}_{i}+{\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{\mathit{A}\left({t}_{i+1}s\right)}\mathbf{B}\left(s\right)\mathit{x}\left(s\right)\mathrm{d}s{\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{\mathit{A}\left({t}_{i+1}s\right)}\mathbf{C}\mathit{x}\left(s\tau \right)\mathrm{d}s.$$(57)
The trapezoidal rule is used in handling the integrals to give:$${\mathbf{x}}_{i+1}={\mathrm{e}}^{\mathit{A}\u2206t}{\mathbf{x}}_{i}+\frac{\u2206t}{2}\left({\mathrm{e}}^{\mathit{A}\u2206t}{\mathbf{B}}_{i}{\mathbf{x}}_{i}+{\mathbf{B}}_{i+1}{\mathbf{x}}_{i+1}\right)\frac{\u2206t}{2}\left({\mathrm{e}}^{\mathit{A}\u2206t}\mathbf{C}{\mathbf{x}}_{ik}+\mathbf{C}{\mathbf{x}}_{i+1k}\right).$$(58a)
Equation (58a) is rearranged to read:$${\mathbf{x}}_{i+1}={{\mathbf{G}}_{i+1}}^{1}{\mathrm{e}}^{\mathit{A}\u2206t}\left(\mathbf{I}+\frac{\u2206t}{2}{\mathbf{B}}_{i}\right){\mathbf{x}}_{i}\frac{\u2206t}{2}{{\mathbf{G}}_{i+1}}^{1}{\mathrm{e}}^{\mathit{A}\u2206t}\mathbf{C}{\mathbf{x}}_{ik}\frac{\u2206t}{2}{{\mathbf{G}}_{i+1}}^{1}\mathbf{C}{\mathbf{x}}_{i+1k},$$(58b)where ${\mathbf{G}}_{i+1}=\mathbf{I}\frac{\u2206t}{2}{\mathbf{B}}_{i+1}$. The matrix form of equation (58b) is:$$\left\{\begin{array}{c}{\mathbf{x}}_{i+1}\\ {\mathbf{x}}_{i}\\ {\mathbf{x}}_{i1}\\ \vdots \\ {\mathbf{x}}_{i+1k}\end{array}\right\}=\left[\begin{array}{cccccc}{\mathbf{M}}_{11}^{i}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{M}}_{1k}^{i}& {\mathbf{M}}_{1,k+1}^{i}\\ \mathbf{I}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right]\left\{\begin{array}{c}{\mathbf{x}}_{i}\\ {\mathbf{x}}_{i1}\\ {\mathbf{x}}_{i2}\\ \vdots \\ {\mathbf{x}}_{ik}\end{array}\right\},$$(59)with the submatrices of the first row reading:$$\begin{array}{c}{\mathbf{M}}_{11}^{i}={{\mathbf{G}}_{i+1}}^{1}{\mathrm{e}}^{\mathit{A}\u2206t}\left(\mathbf{I}+\frac{\u2206t}{2}{\mathbf{B}}_{i}\right)\\ {\mathbf{M}}_{1k}^{i}=\frac{\u2206t}{2}{{\mathbf{G}}_{i+1}}^{1}\mathbf{C}\\ {\mathbf{M}}_{1,k+1}^{i}=\frac{\u2206t}{2}{{\mathbf{G}}_{i+1}}^{1}{\mathrm{e}}^{\mathit{A}\u2206t}\mathbf{C}\end{array}.$$(60)
The FTRM for the delayed damped Mathieu equation is constructed to have the Floquet transition matrix M_{k−1}, M_{k−2}… M_{0} with:$${\mathbf{M}}_{i}=\left[\begin{array}{cccccc}{\mathbf{M}}_{11}^{i}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{M}}_{1k}^{i}& {\mathbf{M}}_{1,k+1}^{i}\\ \mathbf{I}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right].$$(61)
This FTRM is used in computation of stability charts of equation (55) as shown in Figure 7. The stability charts are computed on a 200 by 100 grid of the plane of δ and b using k = 40.
Figure 7. Stability charts of equation (55) from FTRM with ϵ = 1, ω = 1 when (a) damping coefficient κ = 0, (b) κ = 0.1, (c) κ = 0.2 and (d) κ = 0.3. The subdomains of stability are white while those of instability are dark. 
The concept of PA means that equation (57) can be represented as:$${\mathbf{x}}_{i+1}={\mathrm{e}}^{\mathit{A}\u2206t}{\mathbf{x}}_{i}+\frac{1}{\u2206t}{\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{\mathit{A}\left({t}_{i+1}s\right)}\mathrm{d}s{\int}_{{t}_{i}}^{{t}_{i+1}}\mathbf{B}\left(s\right)\mathit{x}\left(s\right)\mathrm{d}s\frac{1}{\u2206t}{\int}_{{t}_{i}}^{{t}_{i+1}}{e}^{\mathit{A}\left({t}_{i+1}s\right)}\mathrm{d}s{\int}_{{t}_{i}}^{{t}_{i+1}}\mathbf{C}\mathit{x}\left(s\tau \right)\mathrm{d}s.$$(62)
The FTRMPA for the delayed damped Mathieu equation is derived from equation (62) as follows:$${\mathbf{x}}_{i+1}={\mathrm{e}}^{\mathit{A}\u2206t}{\mathbf{x}}_{i}+\frac{\u2206t}{4}\left({\mathrm{e}}^{\mathit{A}\u2206t}+\mathbf{I}\right)\left({\mathbf{B}}_{i}{\mathbf{x}}_{i}+{\mathbf{B}}_{i+1}{\mathbf{x}}_{i+1}\right)\frac{\u2206t}{4}\left({\mathrm{e}}^{\mathit{A}\u2206t}+\mathbf{I}\right)\left(\mathbf{C}{\mathbf{x}}_{ik}+\mathbf{C}{\mathbf{x}}_{i+1k}\right).$$(63)
The monodromy matrix is constructed as M_{k−1}… M_{0} where:$${\mathbf{M}}_{i}=\left[\begin{array}{cccccc}{\mathbf{M}}_{11}^{i}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{M}}_{1k}^{i}& {\mathbf{M}}_{1,k+1}^{i}\\ \mathbf{I}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right],$$(64) $$\begin{array}{c}{\mathbf{M}}_{11}^{i}={{\mathbf{G}}_{i+1}}^{1}\left[{\mathrm{e}}^{\mathit{A}\u2206t}+\frac{\u2206t}{4}\left({\mathrm{e}}^{\mathit{A}\u2206t}+\mathbf{I}\right){\mathbf{B}}_{i}\right]\\ {\mathbf{M}}_{1k}^{i}=\frac{\u2206t}{4}{{\mathbf{G}}_{i+1}}^{1}\left({\mathrm{e}}^{\mathit{A}\u2206t}+\mathbf{I}\right)\mathbf{C}\\ {\mathbf{M}}_{1,k+1}^{i}={\mathbf{M}}_{1k}^{i}\end{array},$$(65) $${\mathbf{G}}_{i+1}=\mathbf{I}\frac{\u2206t}{4}\left({\mathrm{e}}^{\mathbf{A}\u2206t}+\mathbf{I}\right){\mathbf{B}}_{i+1}.$$(66)
The stability charts are computed using FTRMPA on a 200 by 100 grid on the plane of δ and b using k = 40 as shown in Figure 8. It is seen that PA has valid application in stability analysis of delayed damped Mathieu equation since results from FTRMPA (Figure 8) are almost the same as those from FTRM (Figure 7).
Figure 8. Stability charts of equation (55) from FTRMPA with ϵ = 1, ω = 1 when (a) damping coefficient κ = 0, (b) κ = 0.1, (c) κ = 0.2 and (d) κ = 0.3. The subdomains of stability are white while those of instability are dark. 
The PTRMPA for the delayed damped Mathieu equation is derived from equation (62) in what follows:$${\mathbf{x}}_{i+1}={\mathrm{e}}^{\mathit{A}\u2206t}{\mathbf{x}}_{i}+\frac{1}{2}\left({\mathrm{e}}^{\mathit{A}\u2206t}\mathbf{I}\right){\left(\mathbf{A}\right)}^{1}\left({\mathbf{B}}_{i}{\mathbf{x}}_{i}+{\mathbf{B}}_{i+1}{\mathbf{x}}_{i+1}\right)\frac{1}{2}\left({\mathrm{e}}^{\mathit{A}\u2206t}\mathbf{I}\right){\left(\mathbf{A}\right)}^{1}\left(\mathbf{C}{\mathbf{x}}_{ik}+\mathbf{C}{\mathbf{x}}_{i+1k}\right).$$(67)
The monodromy matrix for PTRMEA is constructed as M_{k−1}… M_{0} where:$${\mathbf{M}}_{i}=\left[\begin{array}{cccccc}{\mathbf{M}}_{11}^{i}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{M}}_{1k}^{i}& {\mathbf{M}}_{1,k+1}^{i}\\ \mathbf{I}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \cdots & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right],$$(68) $$\begin{array}{c}{\mathbf{M}}_{11}^{i}={{\mathbf{G}}_{i+1}}^{1}\left[{\mathrm{e}}^{\mathit{A}\u2206t}+\frac{1}{2}\left({\mathrm{e}}^{\mathit{A}\u2206t}\mathbf{I}\right){\left(\mathbf{A}\right)}^{1}{\mathbf{B}}_{i}\right]\\ {\mathbf{M}}_{1k}^{i}=\frac{1}{2}{{\mathbf{G}}_{i+1}}^{1}\left({\mathrm{e}}^{\mathit{A}\u2206t}\mathbf{I}\right){\left(\mathbf{A}\right)}^{1}\mathbf{C}\\ {\mathbf{M}}_{1,k+1}^{i}={\mathbf{M}}_{1k}^{i}\end{array},$$(69) $${\mathbf{G}}_{i+1}=\mathbf{I}\frac{1}{2}\left({\mathrm{e}}^{\mathit{A}\u2206t}\mathbf{I}\right){\left(\mathbf{A}\right)}^{1}{\mathbf{B}}_{i+1}.$$(70)
The stability charts are computed using PTRMPA on a 200 by 100 grid plane of δ and b using k = 40 as shown in Figure 9. It is also seen that PA has valid application in stability analysis of delayed damped Mathieu equation since results from PTRMPA (Figure 9) are almost the same as those from FTRM (Figure 7). A more critical look at Figures 7–9 gives that result from PTRMPA is slightly closer to result from FTRM than result from FTRMPA. This goes to say that PTRMPA is more accurate than FTRMPA in stability analysis of delayed damped Mathieu equation. This is in agreement with the conclusion drawn in use of PA in stability analysis of milling process. Agreement of stability charts in Figures 7–9 with those of same parameters in the literature [6, 16] should be noted. Since inverse of $\mathbf{A}=\left[\begin{array}{cc}0& 1\\ \delta & \kappa \end{array}\right]$ is needed in computations involving PTRMPA, δ should not be allowed to be zero. The problem is handled here by starting computation from a very small value of δ that is greater than zero. Both FTRM and FTRMPA are free from this problem. Alternatively the MoorePenrose generalized inverse [17] can be used in handling this problem.
Figure 9. Stability charts of equation (55) from PTRMPA with ϵ = 1, ω = 1 when (a) damping coefficient κ = 0, (b) κ = 0.1, (c) κ = 0.2 and (d) κ = 0.3. The subdomains of stability are white while those of instability are dark. 
7. Conclusion
Application of the trapezoidal rule and Taylor expansion theorem led to reformulated discrete maps for stability analysis of milling process and delayed damped Mathieu equation. The reformulation technique which is called Partial Averaging (PA) retains accuracy while simultaneously reducing computational time (CT). Two avenues of PA were unravelled. One involves use of trapezoidal rule in solving for the average of exponential matrix function to give a map called Full Trapezoidal Rule Map with Partial Averaging (FTRMPA). The other involves exact analysis of average of exponential matrix function to give a map called Partial Trapezoidal Rule Map with Partial Averaging (PTRMPA). By numerical simulation of stability and rate of convergence of milling process (both one and two degree of freedom systems) and delayed damped Mathieu equation it is demonstrated that PTRMPA is more accurate and more time conserving than FTRMPA. Further validation of the presented reformulations via PA stems from the results that stability analysis of fullyimmersed milling process using FTRMPA and PTRMPA is seen to be identical with those of fulldiscretization by Ding et al. [1]. The presented method of reformulation via PA is concluded to be specially useful for reducing the cumbersome analysis and CT of the fulldiscretization method of milling stability analysis when analysis gets more complicated by the use of higher order interpolation/approximation theory. This is numerically demonstrated by application to the second order least squares approximated map of [18].
References
 Y. Ding, L.M. Zhu, X.J. Zhang, H. Ding, A fulldiscretization method for prediction of milling stability, International Journal of Machine Tools and Manufacture 50 (2010) 502–509. [CrossRef] [Google Scholar]
 G. Quintana, J. Ciurana, Chatter in machining processes: a review, International Journal of Machine Tools and Manufacture 51 (2011) 363–376. [CrossRef] [Google Scholar]
 J.V. Le Lan, A. Marty, J.F. Debongnie, CIRP (Eds.), A stability diagram computation method for milling adapted to automotive industry, in: Proceedings of the CIRP Second International Conference on High Performance Cutting, Vancouver, British Columbia, Canada, 12–13 June, 2006. [Google Scholar]
 T. Insperger, Stability analysis of periodic delaydifferential equations modelling machine tool chatter, PhD dissertation, Budapest University of Technology and Economics, 2002. [Google Scholar]
 T. Insperger, G. Stepan, Semidiscretization method for delayed systems, International Journal for Numerical Methods in Engineering 55 (2002) 503–518. [CrossRef] [MathSciNet] [Google Scholar]
 T. Insperger, G. Stepan, Updated semidiscretization method for periodic delay differential with discrete delay, International Journal for Numerical Methods in Engineering 61 (2004) 117–141. [CrossRef] [MathSciNet] [Google Scholar]
 T. Insperger, G. Stepan, J. Turi, On the higherorder semidiscretizations for periodic delayed systems, Journal of Sound and Vibration 313 (2008) 334–341. [CrossRef] [Google Scholar]
 C. Henninger, P. Eberhard, Improving the computational efficiency and accuracy of the semidiscretization method for periodic delaydifferential equations, European Journal of Mechanics – A/Solids 27 (2008) 975–985. [CrossRef] [Google Scholar]
 T. Insperger, Fulldiscretization and semidiscretization for milling stability prediction: some comments, International Journal of Machine Tools and Manufacture 50 (2010) 658–662. [CrossRef] [Google Scholar]
 Y. Ding, L. Zhu, X. Zhang, H. Ding, Secondorder fulldiscretization method for milling stability prediction, International Journal of Machine Tools and Manufacture 50 (2010) 926–932. [CrossRef] [Google Scholar]
 Q. Quo, Y. Sun, Y. Jiang, On the accurate calculation of milling stability limits using thirdorder fulldiscretization method, International Journal of Machine Tools and Manufacture 62 (2012) 61–66. [CrossRef] [Google Scholar]
 Y. Liu, D. Zhang, W. Baohai, An efficient fulldiscretization method for prediction of milling stability, International Journal of Machine Tools and Manufacture 63 (2012) 44–48. [CrossRef] [Google Scholar]
 Y. Ding, L.M. Zhu, X.J. Zhang, H. Ding, Numerical integration method for prediction of milling stability, Journal of Manufacturing Science and Engineering 133 (2011) 031005–031009. [CrossRef] [Google Scholar]
 C.G. Ozoegwu, Chatter of plastic milling CNC machine: master of engineering thesis, Nnamdi Azikiwe University Awka, 2011. [Google Scholar]
 P.V. Bayly, T.L. Schmitz, G. Stepan, B.P. Mann, D.A. Peters, T. Insperger, Effects of radial immersion and cutting direction on chatter instability in endmilling, in: Proceedings of IMECE’02 2002 ASME International Mechanical Engineering Conference & Exhibition, New Orleans, Louisiana, November 17–22, 2002. [Google Scholar]
 E. Butcher, B. Mann, Stability analysis and control of linear periodic delayed systems using Chebyshev and Temporal Finite Element Methods, http://mae.nmsu.edu/faculty/eab/bookchapter_final.pdf [Google Scholar]
 C.D. Meyer, Matrix analysis and applied linear algebra, SIAM, Philadelphia, 2000. [CrossRef] [Google Scholar]
 C.G. Ozoegwu, Least squares approximated stability boundaries of milling process, International Journal of Machine Tools and Manufacture 79 (2014) 24–30. [CrossRef] [Google Scholar]
 R. Sridhar, R.E. Hohn, G.W. Long, A stability algorithm for the general milling process, Transactions of the ASME Journal of Engineering for Industry 90 (1968) 330–334. [CrossRef] [Google Scholar]
 I. Minis, T. Yanushevsky, R. Tembo, R. Hocken, Analysis of linear and nonlinear chatter in milling, Annals of the CIRP 39 (1990) 459–462. [CrossRef] [Google Scholar]
 E. Budak, The mechanics and dynamics of milling thinwalled structures, Ph.D Dissertation, University of British Columbia, 1994. [Google Scholar]
 Y. Altintas, E. Budak, Analytical prediction of stability lobes in milling, Annals of the CIRP 44 (1995) 357–362. [Google Scholar]
 E. Budak, Y. Altintas, Analytical prediction of chatter stability in milling – part I: general formulation, part II: application to common milling systems, Transactions of the ASME Journal of Dynamic Systems, Measurement and Control 20 (1998) 22–36. [Google Scholar]
 S.D. Merdol, Y. Altintas, Multi frequency solution of chatter stability for low immersion milling, Journal of Manufacturing Science and Engineering 126 (2004) 459–467. [CrossRef] [Google Scholar]
 C.G. Ozoegwu, S.N. Omenyi, S.M. Ofochebe, C.H. Achebe, Comparing up and down milling modes of endmilling using temporal finite element analysis, Applied Mathematics 3 (2013) 1–11. [Google Scholar]
Cite this article as: Godwin OC & Sam O: Reducing computational requirement of stability analysis of milling by partial averaging. Manufacturing Rev. 2014, 1, 14.
All Tables
All Figures
Figure 1. Mechanical model of milling tool with indications of dynamic (F_{tan,j}(t), F_{norm,j}(t), F_{x}(t), θ_{j}(t)), process (v, B, Ω) and constant (k_{x}, k_{y}, c_{x}, c_{y}, m etc.) parameters. (a) 1DOF model and (b) 2DOF model. 

In the text 
Figure 2. The asymptotic rate of convergence of fullyimmersed 1DOF milling process at the parameter points (a) Ω = 5000 and w = 0.0002 m (stable operation), (b) Ω = 20000 and w = 0.001 m (stable operation) and (c) Ω = 20000 and w = 0.002 m (unstable operation). It is seen that partial averaging does not diminish rate of convergence much. It is also seen that PTRMPA is always slightly more convergent than FTRMPA and thus advocated. 

In the text 
Figure 3. The asymptotic rate of convergence of fullyimmersed 2DOF milling process at the parameter points (a) Ω = 5000 and w = 0.00001 m (stable operation), (b) Ω = 20,000 and w = 0.00001 m (stable operation) and (c) Ω = 20,000 and w = 0.0001 m (unstable operation). It is seen that partial averaging does not diminish rate of convergence much. It is also seen that PTRMPA is always slightly more convergent than FTRMPA and thus advocated. 

In the text 
Figure 4. Stability lobes of (a) fullyimmersed 1DOF milling process using FTRM (CT = 426 s), FTRMPA (CT = 434 s) and PTRMPA (CT = 417 s) placed on same axes to show that accuracy is marginally affected at low spindle speed, (b) halfimmersed 1DOF downmilling process using FTRM (CT = 433 s), FTRMPA (CT = 408 s) (c) PTRMPA (CT = 397 s) placed on same axes to reveal that PA is accurate relative to FTRM, (c) halfimmersed 1DOF upmilling process using FTRM (CT = 439 s), FTRMPA (CT = 404 s), PTRMPA (CT = 395 s) placed on same axes to reveal good accuracy of PA. 

In the text 
Figure 5. Stability lobes of fullyimmersed 2DOF milling process using FTRM, FTRMPA and PTRMPA placed on same axes to show that the presented two types PA retain accuracy. 

In the text 
Figure 6. Stability lobes of fullyimmersed 1DOF milling process using (a) Second Order Least Squares Approximated Fulldiscretization method with Partial averaging (CT = 430 s) and (b) the original Second Order Least Squares Approximated Fulldiscretization method presented in [18] (CT = 442 s). Very good agreement between the two verifies reliability in terms of accuracy of using PA in the Fulldiscretization method. 

In the text 
Figure 7. Stability charts of equation (55) from FTRM with ϵ = 1, ω = 1 when (a) damping coefficient κ = 0, (b) κ = 0.1, (c) κ = 0.2 and (d) κ = 0.3. The subdomains of stability are white while those of instability are dark. 

In the text 
Figure 8. Stability charts of equation (55) from FTRMPA with ϵ = 1, ω = 1 when (a) damping coefficient κ = 0, (b) κ = 0.1, (c) κ = 0.2 and (d) κ = 0.3. The subdomains of stability are white while those of instability are dark. 

In the text 
Figure 9. Stability charts of equation (55) from PTRMPA with ϵ = 1, ω = 1 when (a) damping coefficient κ = 0, (b) κ = 0.1, (c) κ = 0.2 and (d) κ = 0.3. The subdomains of stability are white while those of instability are dark. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.