Stability Analysis of Multi-Discrete Delay Milling with Helix Effects Using a General Order Full-Discretization Method Updated with a Generalized Integral Quadrature

: A tensor-based general order full-discretization method is enhanced with the capacity to handle multiple discrete delays and helix effects leading to a unique automated algorithm in the stability analysis of milling process chatter. The automated algorithm is then exploited in investigating the effects of interpolation order of chatter states and helix-induced terms on the convergence of milling stability lobes. The enhanced capacity to handle the distributed helix effects is based on a general order formulation of the Newton-Coates integral quadrature method. Application to benchmark milling models showed that high order methods are necessary for convergence of the low speed domain of stability lobes while all the numerically stable orders converge in the high speed domain where the ultra-high order methods are prone to numerical instability. Also, composite numerical integration of the helix-induced integrand beyond the usual zero-th order method leads to higher accuracy of stability lobes especially in the low speed domain.


Introduction
The flutes of technically used milling tools have helix shapes which help to minimize impacts with the workpiece hence, at face value, mitigates chatter and improves machined surface quality. Helix shape also plays positive roles in chip evacuation by inducing axial chip motions. In addition to helix shape, other shape/geometric effects have been introduced in milling tools to possibly induce destructive interference that passively mitigates chatter. Such techniques of chatter suppression which have attracted research attention for several decades include: inter-flute non-uniformity of pitch [1], inter-flute non-uniformity of both pitch and helix angles [2], cutting edge serration of helix tools [3], inter-flute non-uniformity of helix angle [4] and on-flute helix angle variation [5]. In order to highlight the aim and novelty of the current work, the rest of this section surveys the computational techniques applied so far in the stability analysis of milling dynamics of non-uniform pitch helix tools (helix tools with multiple discrete delays). An analytical frequency domain approach based on zero-th order Fourier coefficients was used for identification of the stability lobes of milling with non-uniform pitch helix tools and for establishment of an optimal method for pitch angles selection [1]. The method is an update of the zero-th order method proposed in [6]. The same method was used in [7] for establishing pitch angles that maximize stability limits for given chatter frequency, spindle speed, and number of cutting edges. The cluster treatment of characteristic roots was applied to the stability analysis of milling with non-uniform pitch helix tools in [8], and the arising results reflected the same level of accuracy as the method in [1]. Adaptations of the first order semi-discretization method [9] and temporal-finite element method [10] in stability prediction of non-uniform pitch tools at low radial immersion have been compared and validated with time-domain simulations [2]. The semi-discretization method predicted stability lobes adequately across all radial immersions ranges while the temporal finite element method performed better only in the low radial immersion range. In [11], an adaptation of the first order semi-discretization method was used to identify the stability lobes of non-uniform pitch helix tools, and comparison to the method of [1,7] showed some predictive improvement. The first order full-discretization method [12] was updated for stability prediction of non-uniform pitch helix tools considering runout in [13]. The results showed a noticeable region of improved prediction relative to the method in [1].
Zero-th order and first order semi-discretization methods, which were based on Ackermann's approach to control systems with multiple discrete delays, were used for stability lobes identification of helix tools with non-uniform pitch [14]. While the first method showed marginal predictive accuracy relative to [1], the second method showed a noticeable region of improved prediction that was validated with time-domain simulations. A variable-step trapezoidal rule was used to predict the stability lobes of non-uniform pitch helix tools in [15]. Using time-domain simulations, a comparison of the predictive capacity of the numerical integration method with the frequency domain method in [1] and the full-discretization method in [13] showed improved predictive accuracy. Based on the frequency method in [1], a design procedure for optimizing pitch non-uniformity towards minimizing a novel chatter index called "regeneration factor" was present in [16]. The approach in [17] was an adaptation of the first order semi-discretization for non-uniform pitch helix tools which showed a noticeable region of improved prediction relative to the method in [1]. A combined implicit subspace iteration and zero-th order semi-discretization has been proposed for stability analysis of non-uniform pitch helix tools [18]. The results compared well with those in [14,19]. Semi-third order full-discretization method ("semi" is used here because one and third order interpolations of the current and delayed states were used) has been adapted for stability identification of multiple discrete delay tools due to pitch non-uniformity [20]. Consistent with most of the earlier-mentioned semiand full-discretization methods, this semi-third order full-discretization method revealed a noticeable region of improved prediction relative to [1]. In [21], the combined effect of inter-flute non-uniformities of both pitch and helix angles was investigated with an updated first order semi-discretization method. Their results agreed closely with the method of the work [2]. A novel methodology for the design of variable pitch tool geometry that enhances stable sub-domains was established in [22]. The arising enhancements due to the designs were verified with the semi-discretization and the frequency domain methods.
From the computational viewpoint, the following gaps regarding the time domain approaches, which form the foundation of the novelty of this work, are identified; • Low order interpolation methods have been applied in handling milling states thus accuracy optimization has not been investigated holistically. In this regard, only [20] so far ventured further by exploiting a single case-involving third order interpolation of the delayed state-in the domain of hyper-first order methods. • Stability analyses are mostly done manually thus needing re-analysis when interpolation parameters change. • The distributed cutting force on helix tools has been mostly handled as zero-th order variations on discrete depth segments of the tools. Here again, re-analysis is needed when interpolation order of the variations change.
Time-domain methods of milling stability analysis require the protocol of temporal discretization, chatter states interpolation, monodromy matrix construction, and critical eigenvalue tracking. The first three steps are symbolic analyses which are traditionally done manually while the last step is usually a computerized numerical analysis. Any change on the interpolation order will require re-analysis of the four steps. This protocol is hard enough for the cases of low order interpolation of the cutting states of the simplest milling models let alone high order interpolation of the complex milling models with multiple discrete delays and helix effects. It is proposed to generalize and automatize the four steps for chatter stability analysis of non-uniform pitch tools subjected to helix effects in order to avoid case by case manual analysis/re-analysis, and to allow easy investigation of the effect of interpolation orders of the states and helix-induced integrand on stability lobes accuracy. The methods in [23,24] provided the framework for the first least squares based attempt in [25] on partial generalization of monodromy matrix of the full-discretization methods for arbitrary order approximation of the current state (linear approximation was maintained on the delayed state thus the reference to the method as partial generalization). A complete generalization for the single delay milling model without helix effects was attained in [26] using polynomial tensor representation of cutting states. A major adaptation is needed to extend the method of [26] to take care of the more advanced cases having multiple discrete delays due to non-uniformity of pitch and distributed force due to helix angle. It is aimed here to develop the upgraded generalization and automate it to enable for the first time an investigation of the effects of interpolation parameters of milling states and helix force variation on computational accuracy of stability lobes.
This paper is organized as follows. The current section is the introduction which surveys the related works and specifies the novelty of this work. Section 2 presents a modelling of the regenerative dynamics of milling subjected to multi-delay and helix effects. Section 3 presents a novel general order full-discretization algorithm for the stability lobes identification of the developed model. A novel representation of the helix-induced integral with tensor-based general order Newton-Coates composite quadrature approach is formulated in Section 4 for implementation in the programming of the algorithm of Section 3. Section 5 presents numerical simulation of error and stability results, and discusses the effects of interpolation parameters on computational convergence of stability lobes. The conclusions in Section 6 summarize the key results and their implications, and point out the future potential of the research.

Regenerative Dynamics of Milling with Multi-Delay and Helix Effects
Real milling tools are usually compliant in both the feed (x) and feed-normal (y) directions, see Figure 1. Thus, the directional total motions of a compliant node are x(t) = vt + x t (t) + z x (t) and y (t) = y t (t) + z y (t) where z x (t) and z y (t) are the directional regenerative components and x t (t) and y t (t) are the directional periodic motions. Another characteristic of real milling is the helix shape of the cutting edge which helps in chip removal and shock mitigation. This means that the axial depth of cut w is the integral of differential depths; w = w 0 dz where z is the axial coordinate of a general point on the cutting edge. Therefore, x and y components of the cutting force on the tool are integrated from the resolved differential normal (thrust) dF n,j (t, z) and tangential dF t,j (t, z) forces. The scenarios of the force system and geometry for the uniform pitch (single delay) milling tool and the non-uniform pitch (multiple delays) milling tool are shown in planar view in Figure 1. The x and y components of cutting force on the tool become (a) Uniform pitch (single delay) tool (b) Non-uniform pitch (multiple delays) tool The effect of the helix profile of the cutting edge is to induce a lag angle which is seen from Figure 2 to be given over a differential depth dz by dα = 2 tan β D dz. In Figure 2, ACEG represents a cylindrical segment of the surface traced by the tool which is occupied by the cutting edge segment AE. The angular displacement of the j-th cutting edge for j = 1, 2, . . . , N is therefore given as where β is the helix angle. Though the directional forces given in Equations (1) and (2) drive the responses of the tool only in x and y directions, the differential forces from which they are integrated are functions of axial coordinate z due to the lag angle introduced via the angular displacement θ j (t, z). The screen function g j (t, z) which switches on and off the cutting action of the j-th cutting edge according to whether it is cutting or not is given by where ℘ = (sin θ s − sin θ e )/(cos θ s − cos θ e ). The start angle θ s and end angle θ e of the tooth action are given as θ s = 0 and θ e = arccos(1 − 2ρ) for up-milling and as θ s = arccos(2ρ − 1) and θ e = π for down-milling, where the radial immersion is given as ρ = B/D for a tool of diameter D subjected to radial depth of cut B. The non-linear cutting force law reads In the equations, C t is the tangential cutting coefficient associated with the workpiece material properties and tool shape, χ is the ratio of thrust to tangential force and γ is the feed exponent in the cutting force law. The directional feeds f a,x and f a,y are given as f a,x = (x (t) − x t − τ j ) and f a,y = (y (t) − y t − τ j ) where τ j = 30ψ j /πΩ is the time interval between the j-th and (j − 1)-th cutting edges and ψ j is the j-th pitch angle or angle between the j-th and (j − 1)-th cutting edges. Inserting Equations (6) and (7) in Equations (1) and (2) gives where the total feed function Q T (t, z) is given by The abbreviations used in Equation (10) are as follows; The non-regenerative feed function Q P (t, z) is obtained by setting ∆z x,τ j and ∆z y,τ j in Equation (11) to zero giving Q P (t, z) = vτ j s j γ + γ ∆x t,τ j s j + ∆y t,τ j c j vτ j s j γ−1 . The regenerative feed function Replacing Q T (t, z) in Equations (8) and (9) with Q R (t, z) gives the directional regenerative forces By Newton's third law, these regenerative forces equally excite the tool and the workpiece, meaning that the linear second-order model of either the tool or the workpiece, depending on which is flexible, can be constructed as where and z(t) ∈ R n d is the state, M ∈ R n d ×n d is the mass matrix, C ∈ R n d ×n d is the damping matrix, K ∈ R n d ×n d is the stiffness matrix. The matrices B H ∈ R n d ×2 and C H ∈ R 2×n d serve to project H(t) which is 2-dimensional to the n d -dimensional response. For the case of the tool or rigid workpiece mounted on flexible clamps, only point responses in the feed and feed-normal directions are possible as k x , m x and c x for the feed direction and as k y , m y and c y for the feed normal-direction. For the case of a flexible workpiece, the matrices M, C and K are the assemblages of the elemental global mass matrices. Equation (15) is transformed to the first-order form with the substitutions x 1 (t) = z(t) and where Equation (18) is a periodic multiple discrete delay model with period T = 60/Ω.

A Generalized Stability Analysis Considering Multiple Discrete Delays and the Helix Tool Angle
The period T of the system is divided into This means that (18) is integrated between the limits t i and t i+1 to become The state x(t) is replaced with the polynomial tensor a T (t)T −1 Sv of order (degree) p c , where a(t) is the vector of polynomial basis. The numerical matrices T ∈ R (p c +1)×(p c +1) and S ∈ R (p c +1)×(p c +1) , and the vector of discrete milling states v ∈ R 2n d ×(p c +1) are given as In a similar vein, the delayed states Adapting the method of generalized solution for the non-helix single discrete delay system presented in [26] to the current case of multiple discrete delay, the generalized solution of Equation (22) upon insertion of Equations (23) to (28) will read where and The matrices P and D p d ,p d +2+q are still exactly the same as for the single delay case in [26,27] but additional complication arises from the fact that the dimension of the discrete map is now determined by k m = max(k) which is the largest element of the vector k = {k 1 k 2 · · · k N } T having the discretization integers k j as the elements. The dimension is 2(k m + 1) such that a transition from state It must be noted that while M shows that there are p c current discrete states. Therefore, the current state transition matrix M (p c ) i is given as where, according to [26], the top-most p c submatrices are given as Each of I d,j in Equation (29) has p d + 1 delayed discrete states. The coefficients of the states occupy the top-most p d + 1 locations from k j + 1 − p d to k j + 1 in the j-th delayed state transition matrix M (k j ,p d ) given as If m cn = 1, 2, . . . , k m + 1 denotes the column number in M (k j ,p d ) i , then to ensure proper positioning of the top-most submatricies associated with the delayed sate in the 2(k m + 1)-dimensional discrete transition matrix, the screening parameter f j introduced in Equation (36) is defined as At the columns where f j = 1, the top-most sub-matrices in M (k j ,p d ) i are given as Since the discrete transition is carried out on every of the k equal discrete time intervals [t i , t i+1 ] to cover the spindle period T of the system, the monodromy matrix for the system is constructed as It must be noted that non-uniform pitch lengthens the span of discrete transition and is thus extending the computational time. The models and discrete map established above can still apply to real helix tool milling with uniform pitch having single discrete delay τ = 60/NΩ when the simplification τ j = τ is made for all j, but unnecessary computational time is incurred. The model for such tools isẋ(t) = Ax(t) + B(t)x(t) − B(t)x(t − τ) where A and B(t) are as given in Equations (19) and (20). Thus, the unified algorithm for stability analysis of single delay milling presented in [26] is applicable once the integrated specific force variation matrix H j (t) as given in Equation (16) is used in place of the corresponding quantity in [26].

General Order Tensor-Based Interpolation and Integration of Helix-Induced h j (t, z)
The level of accuracy in handling the helix-induced integrand h j (t, z) will affect the overall accuracy of the stability lobes of helix tools. It has been argued in [15] that the use of numerical integration in handling the integral is contributory to the higher accuracy recorded relative to the semiand full-discretization methods. Most of the earlier works such as described in [2] approximate the corresponding integral as a constant (zero-th order variation) within every depth element of the tool. Here, a general order polynomial tensor variation of the helix-induced integrand is applied such as to investigate the effect of variation of order on the accuracy of stability lobes. In other words, the Newton-Coates integral quadrature approach is generalized for use in automatic integration of h j (t, z) for formation of B(t) and B j (t) at every discrete time step of the generalized discrete mapping presented above for stability analysis.
The numerical integration requires discretization of the axial depth of cut w into K depth intervals [z q , z q+1 ], where q = 0, 1, 2, . . . , K − 1 and ∆z = w K = z q+1 − z q . Then, p h -steps numerical integration rules are generated The helix-induced integrand h j (t, z) is then interpolated with general order polynomial tensors according to Equations (23) to (25). The integration is implemented to give the general order Newton-Coates integral quadrature method The K/p h steps of numerical integration rules are summed to generate the composite numerical integration methods which are represented generally as and for p h ≥ 2 where C r = C p h −r for r = 0, 1, 2, . . . , p h , and these coefficients are exactly the coefficients in the interpolation results from Equations (40) and (41). For illustration, the row vectors of the coefficients C p h = {C 0 C 1 C 2 · · · C p h } for the zero-th to the sixth order integration methods are given as These coefficients, which derive from zero-th to hexic order polynomial tensor interpolation of h j (t, z), are compared for accuracy in the following numerical simulations.

Numerical Simulation and Discussions
The application of the presented generalized algorithm for identifying the stability boundaries of chatter of a 1DOF tool and a 2DOF tool is demonstrated here. The test cases are drawn from literature [14]. More complicated models of the tool are possible but are not in the focus of this algorithmically oriented contribution. In addition to achieving a fully automated algorithm for stability analysis of the milling model considering the complications of helix and multiple delay effects, attention is additionally focused on exploiting the developed algorithm in assessing the sensitivity of computational precision on interpolation orders of states, and also the sensitivity of precision on the method of handling the distributed regenerative force.

Two Degree of Freedom Tool
The values in [14] for the parameters of a 2DOF tool are adopted for numerical simulation here. The parameters of the 2DOF tool are close to that originally established experimentally in [1] and are summarized in Table 1. The parameters are used to study the effect of interpolation order of chatter states and helix-induced integrand on the convergence rate and computational precision of the spectral radii Ψ(p c , p d ) and, thus, the precision of the stability lobes of the adopted milling process. By virtue of the choice of modal parameters provided in Table 1

Convergence with Interpolation Orders
The effects of the interpolation orders of the current state p c and the delayed state p d on the spectral radii computation error (SRCE) is investigated at the process parameter coordinate of spindle speed and axial depth of cut (1000 rpm, 6 mm) by treating the error as a bivariate function of p c and p d . The percentage error is given on a logarithmic scale to base ten as where the computed spectral radius µ SR (k) at a low value of k = 160 is compared to the so-called exact spectral radius µ ESR computed at a high value of k = 400 designated k R . The error surface for the process parameter point (1000 rpm, 6 mm) is presented in Figure 3a which shows a gradual rise with approach of either the current state order p c or the delayed state order p d to zero and sudden rise, when either one of the orders crosses a threshold large value, numerical instability being more sensitive to p c than p d . A planar view of the error surfaces is shown in Figure 3b with a colour bar to show that the domain of minimal error for the simulated point is defined by the order ranges p c ≤ 13 and p d ≤ 16. A magnified view of the undulating base of the SRCE surface which is seen to have several local minima is shown in Figure 3c. The global minimum for the considered process parameter point (1000 rpm, 6 mm) occurs at (p c , p d ) = (1, 3), and the other six local minima in ascending order are (p c , p d ) = (7, 1), (1,10), (1,9), (7, 2), (1, 2) and (7,8). The local minima of other process parameter points occur at somewhat different order combinations but on statistical basis, p c in the neighborhood of 7 and p d in the neighborhood of 2 give the best results for the studied milling system.

Stability Lobes Identification
Eigenvalues are computed on a 200 by 200 gridded plane of the spindle speed range [1000, 10,000] rpm and appropriate axial depth of cut ranges, and used for tracking the stability boundary curves under the criterion of neutral stability. The computational flow for tracking of stability lobes using the established algorithm is explained using a flowchart in the Appendix A. For simplicity of the flowchart, the loop for numerical integration needed for computing H j (t) and the loop needed for summing B(t) from all the cutting edges are not shown. A reference stability lobe is computed using p = 3 and k = 320 as a benchmark for the rest of the stability lobes which are computed using k = 160. The stability lobes of the square unification orders p = 1 to 14 (that is, the cases for which p = p c = p d ) are plotted in Figure 4. Figure 4a shows a full scale plot of the stability lobes for the lower orders p = 1 to 7 while Figure 4b shows the results for the higher orders p = 8 to 14. It can be seen that at p = 14, the method fails at the high speed. The failure continues to extend to the lower speed at higher p values. A magnification of a low speed portion of the stability lobes in Figure 4c shows that high order interpolation is needed for accuracy in the low speed domain since the stability lobes for p = 5 to 14 converged more to the reference compared to the methods from p = 1 to 4. The orders p = 5 and 6 are the most accurate for this case. This result agrees with Figure 3a,b for convergence of spectral radius of an isolated process parameter coordinate with order p. Figures 5  and 6 are generated to study the effects of independent variations of p c and p d on the convergence of the stability lobes of the studied multi-delay and helix milling model. Figure 5 shows the results for a variation of p c from 0 to 14 (Figure 5a shows results for p c = 0 to 7 while Figure 5b shows results for p c = 8 to 14), when p d is fixed at 3. For a study of the effects of variation of p d , Figure 6 shows the stability lobes from the variation of p d from 0 to 14 (Figure 6a shows results for p d = 0 to 7 while Figure 6b shows results for p d = 8 to 14), when p c is fixed at 3. From comparing Figures 5 and 6, it can be seen that while p c = 0 is too erroneous to appear on the full computational plane, p d = 0 is accurate enough to appear on the computational plane. These results agree with the trends of spectral radii computation error in Figure 3a,b. Comparison of Figures 5b and 6b (especially the cases of p c = 14 and p d = 14), shows that the high speed numerical instability stems more from high p d than from high p c . Also, a comparison of Figures 5c and 6c shows that p c plays a higher role in convergence of stability lobes than p d . In fact, the best low speed convergence during variation of p d occurred here when p d = 2 meaning that high order interpolation of the delayed state is not advisable. This is obverse to the conclusion deducible from Figure 5c that high p c is needed for convergence (here, p c ≥ 4). The enlarged view of a high speed portion of each of stability lobes in Figures 4-6 are given in Figure 7. The figure shows that at high speed, stability lobes of all the numerically stable maps for orders greater than one converge. It is noteworthy that in a part of the high speed lobes, p c = 1 resulted in the highest convergence, highlighting that the most accurate stability lobe is multi-order in the sense that it should be built from different sub-ranges computed from different combinations of interpolation orders that are most accurate for the sub-ranges.        The effect of interpolation order p h of the helix-induced integrand h j (t, z) on convergence of the stability lobes of the studied system is investigated by generating the stability lobes for p h = 0 to 6 using the parameters p c = p d = 3, K = 24 and k = 160. The reference stability diagram is computed using p c = p d = 3, K = 24 and k = 320. As seen in Figure 8a for the lower speed range, the most accurate stability lobes occurred at p h = 1, 2, 4 to 6 while the worst results occurred at p h = 0 and 3. Therefore, against the usual practice in earlier works of using zero-th order composite numerical integration to handle the helix-induced integral, the generalized composite numerical integration presented here, which allows automatic higher order treatment of the integral, leads to higher accuracy. The effect of p h is the same in the high speed range as seen in Figure 8b. The results in Figure 9 show that the number of discretization intervals K of axial depth of cut does not have a significant effect on accuracy of the stability lobes. Therefore, in light of the need for reduced computational time without increasing inaccuracy, the computational parameters recommended for the studied system are p c = 5 or 6, p d = 2, p h = 1 and K = 6. (b) A higher speed range Figure 8. Stability lobes of the system for different integration orders of the helix-induced integrand p h = 0 to 6. (b) A higher speed range Figure 9. Stability lobes of the system for different discretization intervals K of axial depth of cut.

Single Degree of Freedom Tool
Here, single degree of freedom tool chatter is considered in the feed direction. This simply implies that in Equations (19) to (21), M = m x , K = k x , C = c x , 0 = 0, I = 1, B H = 1, C H = 1, H j (t) = H j (t) = − C t γ(vτ j ) (γ−1) w 0 h j (t, z)dz and h j (t, z) = g j (t, z)s γ−1 j (χs 2 j + s j c j ). A benchmark single degree of freedom tool with parameters summarized in Table 2 is adopted from [14] for numerical investigation here. A computational grid of 200 by 200 for the process parameter plane of spindle speed range [500, 5000] rpm and axial depth range of [0, 80] mm is adopted. Stability boundary curves are tracked on the plane at the discretization parameter k = 200. The curves for the cases of p = p c = p d = 1 to 13, where are plotted in Figure 10a on the full axes. Computational accuracy is judged relative to the reference curve (red) computed using p = 3 at discretization level of k = 400. The matlab command, trapz, for numerical integration is used here to handle the helix induced integral. It can be seen that both the profile and the parameter range of the stability boundary curves agree with the results for the same system in [14,18]. A low speed range and a higher speed range of the curves are given in Figures 10b,c. It is seen that the first order method clearly stands out from the rest in terms of error. It can be seen from Figure 10a that the ultra-high orders p > 10 are not appropriate in the high speed range while it can be seen from Figures 10b,c that the first order method is not appropriate in the low speed range. A stable island in Figure 10a (which is also seen in the equivalent result in [14]) is selected for validation with time domain simulations. The process parameter coordinate (1000 rpm, 55 mm) which is labeled B lies within the stable island thus expected to be stable while the coordinate (1000 rpm, 4 mm) which is labeled C and the coordinate (1000 rpm, 70 mm) which is labeled A lie below (and also below the main stability curve) and above the B, and are expected to be stable and unstable. The results in Figure 11 confirm the expectations thus validating the stability curve. The matlab integrator dde23 was used to generate the time domain regenerative responses. The parameters used in the simulation are k = 500, K = 24, x(t = 0) = {10 −6 0} T and p h = 2. By the virtue of the computerized generalization, the recommendations which may be different for other systems and process parameter ranges can be automatically revised. (c) A high speed range Figure 10. Stability lobes of the 1DOF system under square unification maps for p = 1 to 13.

Conclusions
In this work, an algorithm based on the arbitrary order full-discretization method and arbitrary order Newton-Coates integral quadrature is developed and implemented for fully automated identification of the stability lobes of multiple delay milling processes with helix effects. Because of the flexibility provided by the generalized algorithm, the study is able to asses the influence of interpolation orders of the cutting states and the helix-induced variation on the identification accuracy of the stability lobes of case studies of 1DOF and 2DOF milling systems. It is found that high order interpolation (p = 5 to 14) is needed for accuracy in the low speed domain. With respect to the separate effects of p d and p c on the stability lobes precision, it is found that high speed numerical instability stems more from high p d than high p c while p c is more influential in the overall convergence of stability lobes than p d . High order interpolation of the delayed state is not advisable for the studied 2DOF case as the best low speed convergence during the variation of p d occurred when p d = 2. Stability lobes of all the numerically stable interpolation orders greater than one converge at high speed. It is concluded that the most accurate stability lobe is multi-order in the sense that it should be built from different sub-ranges computed from different combinations of interpolation orders that are most accurate for the sub-ranges. Increasing the order of interpolation p h of the helix induced function can lead to rise and fall of accuracy in the range of the variation of p h . The number of discretization intervals K of axial depth of cut does not have a significant effect on accuracy of the stability lobes. Recalling that the model presented equally applies to the chatter vibration of flexible workpiece, this work lays a solid mathematical and numerical ground for the stability analysis of flexible workpiece chatter subjected to multiple discrete delays and helix effects.
Author Contributions: All authors have worked on the conceptualization, methodology, investigation, writing, reviewing, editing. The software was done by C.O. All authors have read and agreed to the published version of the manuscript.