Reducing computational requirement of stability analysis of milling by partial averaging

– 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 lose 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 fully-immersed milling process using FTRMPA and PTRMPA are seen to be identical with those of full-discretization by Ding et al. [1] further validating the presented reformulation. PA is applied on the Second Order Least Squares Approximated Full-discretization method of [18] to illustrate its usefulness for reducing the cumbersome analysis and CT of the full-discretization method when analysis gets more complicated by higher order interpolation/ approximation theory


Nomenclature k x
Stiffness in the feed direction, N/m k y Stiffness in the feed-normal direction, N/m m, m x Modal mass in the feed direction, kg m y Modal in the feed-normal direction, kg c x Damping coefficient in the feed direction, Ns/m c y Damping coefficient in the feed-normal direction, Ns/m x nx Tool natural frequency in feed direction, rad/s x ny Tool natural frequency in feed-normal direction, rad/s n x Tool damping ratio in feed direction, rad/s n y Tool damping ratio in feed-normal direction, rad/s x(t) Tool total motion in feed direction, m y(t) Tool total motion in feed-normal 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 norm;j t ð Þ=F tan;j t ð Þ F x (t) Tool cutting force in feed direction, N F y (t) Tool cutting force in feed-normal direction, N z(t) Regenerative motion in feed direction, m

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 time-varying (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 semi-discretization method [4][5][6][7] appeared and enjoyed wide acceptability because of its incisive applicability to many different kinds of delayed dynamical systems.For example, the semi-discretization 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 Multi-Frequency Solution [24].The method of Multi-Frequency 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 semi-discretization were subsequently proposed [8].The semi-discretization 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 so-called full-discretization method.The full-discretization method turned out to be slightly less accurate than semi-discretization method of same order [9,10] but it is as incisive in revealing various stability features of milling.The breakthrough with the advent of full-discretization method is that considerable amount of CT has been saved.Quo et al. [11] have improved the accuracy of the full-discretization 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 full-discretization method with further improvements in accuracy and CT.Least squares approximated first and second order full-discretization methods have been introduced in [18] to save considerable amount of CT of the original full-discretization 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 Full-discretization method of [18]) for reducing the cumbersome analysis and CT of the full-discretization method especially when analysis gets more complicated by the use of higher order interpolation/approximation theory.

The basis for reformulation
Consider the application of the trapezoidal rule to the definite integral: to give Without loss of generality f (s) could be viewed as an n-dimensional state vector and A as an n • n constant matrix.The quantity Dt = t i+1 À t i is the interval of integration and Also consider the application of the trapezoidal rule to the definite integral: to give Equation (3) should be seen to mean definite integral of f (s) on interval t i ; t iþ1 ½ weighted with the mean of e A t iþ1 Às ð Þ and that constitutes the so-called technique of Partial Averaging (PA) of integration scheme.Making use of the matrix exponential expansion e AÁt ¼ I þ P ' p¼1 Át p A p ð Þ=p! in equations ( 2) and (4) respectively gives: The vector q 2 is multiplied out and re-written to become: Equation ( 5) is substituted in equation (7) to give:  8) means that the error of PA is of the order of Dt 2 thus negligible when Dt becomes very small due to shrink of the interval t i ; t iþ1 ½ and due to very small size of the difference f i+1 À f i .This is symbolically written as q 2 À q 1 ¼ OÁt 2 to mean that the error of PA over the interval t i ; t iþ1 ½ is proportional to the square of the size of the integration step Dt.The next thing is to numerically verify this result for the milling process.

Application of PA in chatter stability analysis of milling process
The dynamical model for regenerative milling process is [1,[9][10][11][12][13]18]: 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 feed-normal directions then d = 2.These cases are illustrated in Figure 1.
The vectors n ðdÞ ðtÞ and n ðdÞ t À s ð Þ respectively mean the state and the delayed state of the system at time t.The constant matrix A ðdÞ contains the time-invariant parameters of the system while the periodic coefficient matrix B ðdÞ t ð Þ captures the periodicity of cutting force of the unperturbed milling process.Dividing the discrete delay s of the system into k equal discrete time intervals t i ; t iþ1 ½ where i = 0, 1, 2, . . .(k À 1) and t i ¼ i s k ¼ iÁt and approximating equation (9) in each discrete interval gives the local model: The basic idea of discretization employed in arriving at equation (10) is that B ðdÞ t ð Þn ðdÞ ðtÞ and B ðdÞ ðtÞn ðdÞ t À s ð Þ are discretized while _ n ðdÞ t ð Þ and A ðdÞ n ðdÞ ðtÞ 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 t i ; t iþ1 ½ reads: The trapezoidal rule is used in solving the Duhamel term such that equation (11) becomes: It should be noted that the first published work on application of Newton-Cortes 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: . . . where: The map for the system covering the whole discrete delay is constructed to become: 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: The concept of PA means that equation ( 11) can be represented as : Application of the trapezoidal rule to each of the two integrals of the second term of equation (18) gives: The monodromy matrix is constructed as M kÀ1 . . .M 0 where : 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: gives This leads to the construction of the monodromy matrix as M kÀ1 . . .M 0 where : 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.

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: The x-resultant of cutting force stemming from the chipforming interaction between the tool and workpiece is seen in Figure 1a to be: where N is the number of teeth on the milling tool indexed with the values j = 1, 2, 3. . .N and X ¼ F norm;j t ð Þ=F tan;j t ð Þ is the ratio of normal to tangential cutting force on the jth tooth.
The instantaneous angular position of jth tooth h j t ð Þ is given by: The screen function g j ðtÞ has values of either unity or nullity depending on whether the tool is cutting or not.Given start and end angles of cut h s and h e , respectively, g j ðtÞ becomes [14]: The precise mathematical form for g j t ð Þ is given as [4]: The screen function g j ðtÞ has been formulated in terms of radial immersion q = B/D and h j t ð Þ for up-milling and down-milling modes respectively as [25]: The tangential cutting force on jth tooth is generally given by the non-linear law [4]: 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 s = 60/(NX)) given as: and c is a positive exponent.Equations ( 29), (34) and ( 35) taken together give: where q(t) is a s( = 60/NX) periodic function given by: The regenerative motion z(t) of the tool is given by: where x t (t) is the s-periodic response due to periodic force of tool-workpiece interaction.Then making use of equations ( 38) and (36) in equation ( 1) when z(t) = 0 gives: Inserting equations (39) and (38) in equation (28) gives: This is linearized to become: where: is the specific periodic cutting force variation.The modal form of governing model chatter of milling process becomes: The natural frequency and damping ratio of the tool are given as . Equation ( 43) is easily put in state space form of equation ( 1) such that: O.C. Godwin and O. Sam: Manufacturing Rev. 2014, 1, 14 In the two degree of freedom case (d = 2) the non-linear force is used in a similar manner taking into account tool compliance in the feed and feed-normal directions to give: 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 and n y ¼ c y =2 ffiffiffiffiffiffiffiffiffi ffi m y k y p .The parameters used for numerical simulation are contained in Table 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 k ð Þ ¼ l k ð Þ j j max with rise in k where f (k) is the absolute value of the maximum-magnitude 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 l k ð Þ j j max at k = 1.The asymptotic rate of convergence of fully-immersed 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 fully-immersed 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.

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) sub-regions 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 Full-discretization by Ding et al. [1].This suggests that PA is strongly applicable to the Full-discretization 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.The technique of PA is even more validated by the stability result of 2DOF fully-immersed milling process since no deviation is seen amongst FTRM, FTRMPA and PTRMPA in Figure 5.

Application of PA in simplifying the second order full-discretization method
The second order least squares approximated polynomial for the state term ñðdÞ t ð Þ reads [18]: Equation (48) together with the linear approximations: is inserted in equation ( 18) and solved to construct the monodromy operator w = M kÀ1, M kÀ2 . . .M 0 such that:   where: 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 full-discretization method in [18].The G matrices in equation ( 54) are each given in terms of single product of some fraction of linear power of Dt and E while the G matrices in [18] are each given in terms of sum of products of various non-linear powers of Dt and F matrices.Making use of presented map with k = 40 on a 200 by 100 grid produces the stability diagram of fully-immersed milling as shown in Figure 6a.The original Second Order Least Squares Approximated Full-discretization 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.

Application of PA in stability analysis of delayed damped Mathieu equation
The delayed damped Mathieu equation is normally given as: This is put in state space as: The time delay is given as s = 2p while the period of the matrix B(t) is given as T = 2p/x.The weight of the delayed feedback state is captured in the gain parameter b.The symbol j is the damping coefficient of the system.Equation ( 55) is divided into k equal discrete time intervals t i ; t iþ1 ½ as done in the case of milling process and integrated between the limits t i and t i+1 to give: The trapezoidal rule is used in handling the integrals to give: Equation (58a) is re-arranged to read: where G iþ1 ¼ I À Át 2 B iþ1 .The matrix form of equation (58b) is: x i x iÀ1 x iÀ2 . . .
with the sub-matrices of the first row reading: 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: 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 d and b using k = 40.[18] (CT = 442 s).Very good agreement between the two verifies reliability in terms of accuracy of using PA in the Fulldiscretization method.
The concept of PA means that equation (57) can be represented as: The FTRMPA for the delayed damped Mathieu equation is derived from equation (62) as follows: The monodromy matrix is constructed as M kÀ1 . . .M 0 where: The stability charts are computed using FTRMPA on a 200 by 100 grid on the plane of d 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).
The PTRMPA for the delayed damped Mathieu equation is derived from equation ( 62) in what follows: The monodromy matrix for PTRMEA is constructed as M kÀ1 . . .M 0 where: The stability charts are computed using PTRMPA on a 200 by 100 grid plane of d 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 A ¼ 0 1 Àd Àj ! is needed in computations involving PTRMPA, d should not be allowed to be zero.The problem is handled here by starting computation from a very small value of d that is greater than zero.Both FTRM and FTRMPA are free from this problem.Alternatively the Moore-Penrose generalized inverse [17] can be used in handling this problem.

Conclusion
Application of the trapezoidal rule and Taylor expansion theorem led to reformulated discrete maps for stability analysis .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 fully-immersed milling process using FTRMPA and PTRMPA is seen to be identical with those of full-discretization 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].

Figure 1 .
Figure 1.Mechanical model of milling tool with indications of dynamic (F tan,j (t), F norm,j (t), F x (t), h j (t)), process (v, B, X) and constant (k x , k y , c x , c y , m etc.) parameters.(a) 1DOF model and (b) 2DOF model.

Figure 2 .
Figure2.The asymptotic rate of convergence of fully-immersed 1DOF milling process at the parameter points (a) X = 5000 and w = 0.0002 m (stable operation), (b) X = 20000 and w = 0.001 m (stable operation) and (c) X = 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 .
Figure3.The asymptotic rate of convergence of fully-immersed 2DOF milling process at the parameter points (a) X = 5000 and w = 0.00001 m (stable operation), (b) X = 20,000 and w = 0.00001 m (stable operation) and (c) X = 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.

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

Figure 6 .
Figure 6.Stability lobes of fully-immersed 1DOF milling process using (a) Second Order Least Squares Approximated Full-discretization method with Partial averaging (CT = 430 s) and (b) the original Second Order Least Squares Approximated Full-discretization 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.