Stability Analysis of Milling Based on the Barycentric Rational Interpolation Differential Quadrature Method

: Chatter causes great damage to the machining process, and the selection of appropriate process parameters through chatter stability analysis is of great significance for achieving chatter-free machining. This article proposes a milling stability analysis method based on the barycentric rational interpolation differential quadrature method (DQM). The dynamics of the milling process considering the regeneration effect is first modelled as a time-delay differential equation (DDE). When adjacent pitch angles of the milling cutter are symmetric, the milling dynamic equation contains a single time delay. Otherwise, when adjacent pitch angles are asymmetric, the dynamic equation contains multiple time delays. The barycentric rational interpolation DQM is then used to approximate the differential and delay terms of the milling dynamics equation, and to construct a state transition matrix between adjacent milling periods. Finally, the chatter stability lobe diagram (SLD) is obtained based on the Floquet theory. According to the SLD, the appropriate spindle speed can be selected to obtain the maximum stable axial depth of cutting, thereby effectively improving the material removal rate. The accuracy and efficiency of the proposed method have been validated by two widely used milling models, and the results show that the proposed method has good accuracy and computational efficiency.


Introduction
Milling is a common machining process, which is widely used in the field of aerospace and other high-end equipment manufacturing because it can achieve high-precision machining of complex surfaces.However, the appearance of chatter on machine tools is disastrous [1].The chatter during the milling process seriously affects the stability of the machining process, thereby reducing the surface quality of the workpiece and shortening the service life of the tool.Therefore, achieving chatter-free milling has become one of the major concerns.Researches have shown that the selection of process parameters has a significant effect on the chatter stability of the milling process, and optimizing the process parameters to achieve chatter-free machining is a simple and effective strategy.
The SLD reveals the chatter stability limits under different process parameters, which can be utilized to optimize the process parameters effectively.Various methods have been proposed to construct the chatter SLD.The chatter stability limits under real working conditions can be obtained by performing cutting experiments.However, experimental-based methods [2][3][4] require extensive cutting experiments when constructing the SLD, which inevitably lead to significant time and economic costs.Therefore, it is necessary to seek some low-cost stability prediction methods.In the milling process, the regenerative effect Symmetry 2024, 16, 384 2 of 18 had been reported as the main reason for chatter [1].The milling process considering the regenerative effect is generally formulated as a DDE.Many stability analysis methods have been generated based on the DDE of the milling process.The time-domain simulation methods [5][6][7][8] can fully consider the nonlinear phenomena of the machining process and provide important parameters such as the cutting force and the vibration displacement while analyzing the stability of the milling process.However, the time-domain simulation still entails significant computational costs.To compensate for these shortcomings, many methods based on dynamic analysis have been proposed.Based on the frequency-domain stability theory, the single-frequency solution [9] and multi-frequency solution [10,11] methods were proposed by Altintas and Budak.According to the characteristics of the intermittent cutting during the milling process, Bayly et al. [12] divided the vibration state of the cutting tool into free vibration and forced vibration, and proposed a temporal finite element analysis (TFEA) method for predicting the stability of the milling process.Subsequently, Mann et al. [13] used the TFEA method to simultaneously obtain the stability of the milling process and the surface location error of the workpiece.Insperger et al. [14][15][16] put forward the semi-discretization method (SDM) to study the delayed system and obtained the SLD of the milling process based on the SDM.The SDM and TFEA were based on the differential equation theory and the variational method, respectively [17] and their main idea was first to approximate the infinite-dimensional delay system to a finite-dimensional one, then construct a transition matrix between adjacent periods, and finally apply the Floquet theory to determine the stability of the system.Similarly, Ding et al. developed the full-discretization method (FDM) [18] based on the direct integration scheme, and adopted the FDM to study the stability of milling processes.Butcher et al. introduced the Chebyshev polynomial method [19] and the Chebyshev collocation method [20] for milling stability analysis.Lin et al. [21] established a generalized regression neural network model to predict the limiting axial cutting depth of the milling process.
Based on the numerical solution technique of integral equations, Ding et al. [22] developed the numerical integration method for analyzing the chatter stability of the milling process.Within the framework of numerical integration, various other methods have also been derived.On the basis of Runge-Kutta methods, Niu et al. [23] proposed the classical fourth-order Runge-Kutta method and the generalized Runge-Kutta method to predict the stability of the milling process.Starting from the perspective of the initial value solution of the differential equations, Zhang et al. [24] derived a chatter stability prediction method based on the Simpson method.Mei et al. [25] developed a cutting chatter stability prediction method based on the linear multistep method.Later, Mei et al. [26] proposed an adaptive variable-step numerical integration method, and used it to construct the SLD of the milling process with multiple delays.Yang et al. [27] proposed a five-point Gaussian quadrature-based chatter prediction method of milling processes together with a transition matrix reduction scheme to improve the computational accuracy and efficiency.
As far as the solution of milling dynamics equations is concerned, in addition to numerical integration strategies, numerical differentiation techniques can also be used for milling stability analysis.Zhang et al. [28] proposed a stability prediction method based on the finite difference method and extrapolation method.Ding et al. [29] studied the stability of the milling process using the DQM.The DQM is essentially a numerical differentiation technique that approximates the derivatives at the nodes by a linear combination of function values at the nodes.Since its proposal by Bellman R. et al. [30,31], the DQM has been widely used to solve differential equations due to its advantages of simple concepts and low computational complexity.However, the classical DQM also has its shortcomings.For instance, it is sensitive to the distribution of discrete nodes, and is prone to generating an ill-conditioned coefficient matrix when the number of discrete nodes is large.To address these issues, Zong Z. et al. [32] introduced the localized differential quadrature method (LDQM) and successfully applied it to the problem of two-dimensional wave equations.Mei et al. [33] derived a chatter stability analysis method based on the LDQM.In their works, the derivative of the state vector in the milling dynamics equation is represented as Symmetry 2024, 16, 384 3 of 18 the weighted sum of state vector values at local nodes, which can generate a sparse weight matrix and improve the stability of numerical calculations.However, when there are too many local nodes, unstable results can still be generated.
Based on the above review, commonly used milling stability analysis methods can be classified and are listed in Table 1.
Table 1.The most commonly used milling stability analysis methods.
From the perspective of numerical analysis, the shortcomings of the classical DQM stem from its construction process of polynomial interpolation.It is well known that rational interpolation sometimes gives better approximations than polynomial interpolation, especially for large sequences of points.Floater and Hormann [34] proposed a family of barycentric rational interpolants that have no real poles and arbitrarily high approximation orders on any real interval, regardless of the distribution of the points.Inspired by their works, this article derives a milling chatter stability analysis method based on the barycentric rational interpolation DQM.The novelty of this work lies in the use of the barycentric rational interpolation to approximate the vibration response function in one cutting cycle and the construction of a DQM based on the barycentric rational interpolation by differentiating the interpolation function.The use of the barycentric rational interpolation DQM to process the milling dynamic equation effectively improves the shortcomings of the classical DQM in producing an ill-conditioned weight coefficient matrix when there is a large number of discrete nodes.
This paper is organized as follows: in Section 2, the dynamics model of the milling process is outlined, taking into account the regeneration effect.In Section 3, a chatter stability analysis method based on the barycentric rational interpolation DQM is presented.In Section 4, the effectiveness of the proposed method is verified, and the verification results are discussed.Finally, conclusions are derived in Section 5.

Dynamics Model of the Milling Process
The milling process considering the regeneration effect is usually modeled by a DDE [9,27], which can be written as follows: where M, C and K are the mass, damping and stiffness matrices of the milling system, respectively.y(t) = x(t) y(t) T and f(t) = F x (t) F y (t) T are the vibration displacement vector and milling force vector, respectively.The calculation of the milling force usually depends on different milling force models.Figure 1 shows a two-degree-of-freedom (DOF) flexible tool rigid workpiece milling model.In Figure 1, φ j (t, z) is the angular position of the cutting edge on the tooth (j) with a height z, ϕ j is the pitch angle between the tooth (j − 1) and the tooth (j), β is the helix angle of the milling cutter, z u,j and z l,j are the upper and lower bounds of the cutting edge participating in cutting on the tooth (j), φ st and φ ex are the start and exit angles of the tooth, ϕ lag (z) is the lag angle of the cutting edge at height z and ϕ s is the tooth sweep angle.
cutting on the tooth ( j ), st ϕ and ex ϕ are the start and exit angles of the tooth, ( ) φ is the lag angle of the cutting edge at height z and s φ is the tooth sweep angle.( ) where D is the diameter of the milling cutter.
( ) ( ) ( ) where ( ) 0, 0 j ϕ is the angular position of the tooth ( j ) at the initial time and Ω is the spindle speed.
The cutting force arises from the chip formation process, which consists of two parts: static thickness formed by the feed motion and dynamic thickness formed by the relative motion between the tool and the workpiece.These two parts of chip thickness form the static cutting force and dynamic cutting force, respectively.Therefore, the milling force f can be expressed as follows: where c f is the static component of the milling force, and d f is the dynamic component of the milling force.Since the stability of the cutting process is influenced only by the dynamic component of the cutting force [35], the static component of the cutting force will be omitted in the subsequent analysis process.Based on the linear cutting force model [36], the milling force vector can be calculated as follows: According to the geometric relationship described in Figure 1, ϕ lag (z) and φ j (t, z) can be calculated through the following formula: where D is the diameter of the milling cutter.
where φ j (0, 0) is the angular position of the tooth (j) at the initial time and Ω is the spindle speed.
The cutting force arises from the chip formation process, which consists of two parts: static thickness formed by the feed motion and dynamic thickness formed by the relative motion between the tool and the workpiece.These two parts of chip thickness form the static cutting force and dynamic cutting force, respectively.Therefore, the milling force f can be expressed as follows: where f c is the static component of the milling force, and f d is the dynamic component of the milling force.Since the stability of the cutting process is influenced only by the dynamic component of the cutting force [35], the static component of the cutting force will be omitted in the subsequent analysis process.Based on the linear cutting force model [36], the milling force vector can be calculated as follows: Equation ( 5) is a multiple-time-delay dynamic milling force calculation model.In Equation ( 5), N is the number of teeth, T j is the time delay of the tooth (j), y t − T j is the vibration displacement vector of the previous tooth-passing period of the tooth (j) and H j (t) is the dynamic milling force coefficient matrix corresponding to the tooth (j).Due to the periodic characteristics of the milling process, H j (t) satisfies the relationship of H j (t) = H j (t + T), where T is the spindle period.H j (t) can be written as follows: where In Equation ( 7) K t and K n are the tangential and radial cutting force coefficients, respectively, and can be determined through cutting experiments.
The stability analysis of the milling process is usually conducted in the state space.We can substitute Equation ( 5) into Equation ( 1) and transform the governing equation of the milling process into the state space form. .

y(t)
T is the state vector and A and B j (t) can be expressed as 8) represents a linear periodic system with time-delay terms.It is worth noting that when adjacent pitch angles are symmetric, the values of T j (j = 1, 2, • • • , N) are the same and the system degenerates into a single-time-delay system.Otherwise, the system is a multipletime-delay system.According to the Floquet theory, the linear periodic system is stable if and only if the spectral radius of its Floquet transition matrix is less than one.

DQM Based on the Barycentric Rational Interpolation
In the mathematical field of numerical analysis, interpolation is a type of estimation, a method of constructing new data points based on the range of a discrete set of known data points.In engineering and science, polynomials are very often used for interpolation because of their straightforward mathematical properties [37].The most commonly used methods for constructing interpolation polynomials include Lagrange interpolation and Newton interpolation.According to polynomial interpolation theory, if (x 0 , f (x 0 )), (x 1 , f (x 1 )) . . . (x n , f (x n )) are n + 1 points in the plane with distinct x i , then there exists one unique polynomial p of degree n or less that satisfies p(x i ) = f (x i ) for i = 0, • • • , n.According to the Lagrange interpolation, p(x) can be expressed as follows: where and the Lagrange interpolation basis function l i has the following property: The formula described in Equation ( 10) is often referred to as the classical form of Lagrange interpolation.Some improvements can be made to Equation (10) to improve the characteristics of the classical Lagrange interpolation.By defining the barycentric weight w i , Equation ( 10) can finally be expressed in the following form [38]: The formula described in Equation ( 14) is called the second (true) form of the barycentric interpolation formula.The use of the barycentric formula can effectively reduce the computational complexity of classical Lagrange interpolation and improve numerical stability.

Barycentric Rational Interpolation and Its Differentiation
The interpolation formula proposed by Floater and Hormann [34] can be expressed as follows: where d is any integer with 0 ≤ d ≤ n, and for each i = 0, • • • , n − d, p i (x) is the unique polynomial that interpolates f (x) at points x i , . . . x i+d .The blending function λ i is defined as follows: By introducing a new barycentric weight w j , Equation ( 15) can be written as the following barycentric form: where the barycentric weight w j is expressed as follows: where J j is an index set that satisfies J j := {i ∈ I : j − d ≤ i ≤ j} and The barycentric rational interpolation formula represented by Equation ( 17) can be used to evaluate the derivatives of r(x) based on the formulas proposed by Schneider and Werner [39].Equation ( 17) can be further written as follows: where Taking differentiation with respect to x on both sides of Equation ( 19) and substituting x = x i the following equation can be obtained: Equation ( 21) provides a DQM based on the barycentric rational interpolation.By letting A ij = q ′ j (x i ), the differentiation of r(x) at point x i (i = 0, • • • , n) can be expressed in the following matrix form: . . .
where A ij can be calculated using the following formula:

Stability Analysis of the Milling Process Based on the Barycentric Rational Interpolation DQM
According to the Floquet theory, when conducting stability analysis of the milling process, it is necessary to first construct the Floquet transition matrix between adjacent milling periods.To construct the transition matrix, one spindle period [0, T] is discretized into m small pieces using m + 1 different time nodes first.The discrete time node can be written as t i (i = 0, 1, • • • , m).Similarly, the previous spindle period [−T, 0] is discretized into m pieces using the same set of discrete points, and the discrete nodes are represented as t i−T (i = 0, 1, • • • , m).Then, the state Equation (8) at time node t i (i = 0, 1, • • • , m) can be expressed as follows: For the convenience of description, let x i denote x(t i ), x i−T denote x(t i − T), x i−T j denote x t i − T j and B j,i denote B j (t i ), the state Equation ( 8) at all the discrete time nodes can be expressed as follows: According to the DQM provided by Equation ( 21), the derivative of the state vector on the left side of Equation ( 25) can be expressed as follows: where W pq (p, q = 0, 1, • • • , m) is the corresponding weighted coefficient and can be calculated according to Equation (23).Next, the state vector x t i − T j can be represented by a linear combination of the value of According to the position of the time node t i and the length of the delay term T j , four different situations about the location of the time node t i − T j in adjacent spindle periods are shown in Figure 2.
In Figure 2a, the time node t i − T j belongs to [−T, 0] and does not coincide with discrete nodes t k−T (k = 0, 1, • • • , m); in Figure 2b, the time node t i − T j belongs to [−T, 0] and coincides with a certain discrete node t k−T (k = 0, 1, • • • , m); in Figure 2c, the time node t i − T j belongs to [0, T] and does not coincide with discrete nodes t k (k = 0, 1, • • • , m); in Figure 2d, the time node t i − T j belongs to [0, T] and coincides a certain discrete node t k (k = 0, 1, • • • , m).The calculation of state vector x t i − T j needs to be carried out according to the four different situations shown in Figure 2. If the time node t i − T j is located at the position described in Figure 2a or Figure 2c, the state vector x t i − T j can be interpolated using the barycentric rational interpolation formula given in Equation (19).If the time node t i − T j is located at the position described in Figure 2b or Figure 2d, the state vector x t i − T j is directly equal to the function value of the discrete point x(t k−T ) or x(t k ).Based on the above analysis, the value of x t i − T j can be expressed as follows: In Figure 2a, the time node i − is located at the position described in Figure 2a or Figure 2c, the state vector ( ) can be interpolated using the barycentric rational interpolation formula given in Equation (19).If the time node i j t T − is located at the position described in Figure 2b or Figure 2d, the state vector ( ) is directly equal to the function value of the discrete point ( ) . Based on the above analysis, the value of ( ) can be expressed as follows: Equation ( 27) can also be expressed in the form of vector multiplication as shown below.(28) For one specific time delay T j , the position of t i − T j will have two situations as the index i increases: t i − T j < 0 or t i − T j ≥ 0. Assuming that a certain index s satisfies t s − T j < 0 and t s+1 − T j ≥ 0, we have the following expression: Substituting Equations ( 26) and ( 29) into Equation ( 25) and carrying out a simple transformation yields the following expressions: . . .
Symmetry 2024, 16, 384 Simplifying and rearranging Equation ( 30) yields the following equation: Finally, the Floquet transition matrix Φ between adjacent milling periods can be obtained based on Equation (32) as follows: The stability of the milling process can be determined by the spectral radius of the Floquet transition matrix Φ according to the Floquet theory.

Numerical Validation and Discussion
The calculation accuracy and efficiency of the proposed method are verified on account of two benchmark milling models in this section.

Single-Time-Delay Milling Model
Without loss of generality, this section utilizes the two-DOF milling model same as that in References [15,22] for demonstration.The cutting tool is a 12.7 mm diameter, two-flute, carbide helical end mill with a 106.2 mm overhang held in an HSK 63A collet-type tool holder.This milling model has undergone extensive validation and is widely used as a benchmark example to validate algorithms.A two-flute cutter with symmetric pitch angles is used in this model, so its governing equation contains a single time delay and can be converted to the state space form as shown below. .
where T is equal to the tooth passing period, and where ω n is the nature angular frequency, ζ is the relative damping, w is the axial depth of cutting and m t is the modal mass.The values of the dynamics parameters of the milling model are shown in Table 2.The SLD of the milling model is obtained using the method derived in Section 3. It is worth noting that the method in Section 3 is derived based on a generalized multiple time-delay milling model and can be appropriately simplified when used for the singletime-delay milling model.Here, the first-order SDM (1st-SDM) is adopted as a benchmark to verify the effectiveness and efficiency of the proposed method, as the SDM has been experimentally validated and widely applied.In the simulation experiment, the cutting method is down-milling, the spindle speed ranges from 5000 to 10,000 rpm, the axial depth of cutting ranges from 0 to 6 mm, and the radial immersion ratio ranges from 0.2 to 1.0.In numerical simulation, the milling period is discretized using the uniform nodes.The SLD obtained is constructed on a 200 × 100-sized grid.The comparisons of the 1st-SDM, LDQM and the proposed method with different parameters are shown in Figures 3 and 4. The reference stability limits represented by the red line in Figures 3 and 4 3a.As mentioned earlier, when the number of discrete nodes is large, the classical DQM produces an illconditioned differential matrix, resulting in erroneous stability prediction results.Figure 3b shows the SLD obtained by the LDQM [33], from which it can be seen that when the parameter l is too large, the LDQM may still yield unstable results.In Figure 3c,d, the SLD is obtained by the method proposed in Section 3. As shown in Figure 3c,d, with the increase of the discrete parameter m, the SLD obtained by the method proposed in this work shows good consistency with the reference value, indicating that the method proposed in this work has good accuracy.Compared with the LDQM using l = 21 local nodes (Figure 3b), the method proposed in this paper can obtain accurate results when the discrete parameter m is 60 (Figure 3d).This indicates that the DQM based on barycentric rational interpolation effectively improves the shortcomings of the classical DQM that can easily produce an ill-conditioned matrix when the number of discrete nodes is large.Figure 4a-d show the SLD obtained by the proposed method with two other different radial immersions.In Figure 4a,b, the radial immersion is 0.6, and the discrete parameters (m) are 20 and 30, respectively.In Figure 4c,d, the discrete parameters (m) are 20 and 30 and the radial immersion is 0.2, which is a typical machining condition of the small radial depth of cutting.From Figure 4, it can also be concluded that the method proposed in this work can obtain the accurate SLD to predict the stability of the milling process.
From Figures 3 and 4, it can be seen that when the cutter with symmetric pitch angles is used, the maximum value of the stability limit appears in the right area (near the spindle speed of 9200 rpm) of the SLD.We can choose the optimal axial depth of cutting to achieve the maximum material removal rate.At the same time, it can also be concluded that the method proposed in this work is not affected by the machining conditions and can obtain accurate SLD under both large and small radial depths of cutting.
In order to further validate the efficiency of the method, the time required to calculate the SLD in numerical simulation is also compared with that of the 1st-SDM based on the above mentioned two-DOF milling model.The simulative calculations were performed in Matlab ® on a laptop computer (Intel ® Core TM i7-10870H CPU 2.20 GHz, 16.00 GB RAM).
The machining parameters used in numerical simulation are the same as those used in calculating the SLD in Figures 3 and 4. Here, the discrete parameter m is set as 60.The calculation time of the 1st-SDM and the proposed method are listed in Table 3. From Table 3, it can be seen that under three different radial immersion conditions, the 1st-SDM requires the longest calculational time.In addition, as the radial immersion decreases, the time consumption of the 1st-SDM also decreases accordingly.The reason for this phenomenon is that as the radial immersion decreases, the free vibration time of the cutting tool increases in one cutting period.In the free vibration range, the value of parameter B in the milling state equation is 0, which reduces the complexity of the numerical calculation and saves some calculation time.As for the method in this work, it can be seen from the derivation process that the proposed method approximates the differential and delay terms of the state equation when constructing the transition matrix, and ultimately approximates the state equation into an algebraic system of equations.The entire process has low computational complexity, thus saving computational time.From the calculation time listed in Table 3, it can be seen that compared with the semi discrete method, the method proposed in this paper can save up to 87% of the calculation time, indicating that the method proposed in this work has good computational efficiency.
Also, it can be seen From Table 3 that with the increase of the barycentric rational interpolation parameter d, the calculation time of the method proposed in this work slightly increases, but overall, the impact of the change in parameter d on the calculation time is negligible.
It is worth noting that although uniform nodes are used in the validation of the proposed method in this section, the barycentric rational interpolation itself has no restriction on the type of nodes, and other types of nodes can be just as effectively used in the method proposed in this paper.

Multiple-Time-Delay Milling Model
In this section, a multiple-time-delay milling model is used to verify the method derived in this work.For the convenience of comparison and verification, the milling model used in [40,41] is employed for numerical simulation.A 19.05 mm diameter variable pitch cutter with four flutes is adopted.The helix angle is 30 • and the asymmetric pitch angles are 70 • -110 • -70 • -110 • .Due to the use of a variable pitch cutter, the milling process is a typical multiple-time-delay milling model as described by Equation (8).The detailed values of the modal parameters and cutting force coefficients of this milling model are listed in Table 4.The SLD of this multiple-time-delay milling model is constructed over a 200 × 100-sized grid, the machining condition is down-milling and the simulation parameters are set as follows: the spindle speed Ω ranges from 2500 to 12,500 rpm, the axial depth of cutting w ranges from 0 to 10 mm and the radial immersion ratio ranges from 0.2 to 1.0.The milling period is also discretized using the uniform nodes.The results are shown in Figure 5; in the SLD, the reference stability limits demoed by the red line are obtained by the 1st-SDM with the discrete parameter m = 200.In Figure 5a-c, the radial immersions are 1, 0.6, and 0.2, respectively.It can be seen from the SLD displayed in Figure 5a-c that when the cutter with asymmetric pitch angles is used, the maximum value of the stability limits does not appear in the right area of the SLD, but in the middle area (near the spindle speed of 5400 rpm), which is different from using a cutter with symmetric pitch angles.Also, the SLD displayed in Figure 5a-c agrees well with the reference SLD.This shows that the milling stability analysis method proposed in this work is also suitable for multiple-time-delay milling models.

Conclusions
This work mainly studies the analysis of the chatter stability in the milling process.A method for analyzing the chatter stability of the milling process is derived based on the barycentric rational interpolation DQM.Firstly, the milling process considering the regeneration effect is modeled as a DDE.Then, the differential and delay terms of the milling state equation are approximated as a linear combination of discrete state vectors using the barycentric rational interpolation DQM.The state equation is approximated as a system of linear equations.Finally, the state transition matrix on adjacent milling periods is constructed, and the stability of the milling process is studied based on the Floquet theory.The accuracy and computational efficiency of the method are validated using two widely used milling models.The results of the simulation experiments show that the method proposed in this work has good accuracy and computational efficiency.The main characteristics of the method proposed in this article are listed as follows: 1. Using the barycentric rational interpolation DQM can effectively improve the shortcomings of the classical DQM, avoiding the generation of the ill-conditioned matrix when there are a large number of discrete nodes, thereby improving the stability and accuracy of numerical calculations.2. The proposed method approximates the state equation of the milling system to an algebraic system of equations through interpolation and numerical differentiation

Conclusions
This work mainly studies the analysis of the chatter stability in the milling process.A method for analyzing the chatter stability of the milling process is derived based on the barycentric rational interpolation DQM.Firstly, the milling process considering the regeneration effect is modeled as a DDE.Then, the differential and delay terms of the milling state equation are approximated as a linear combination of discrete state vectors using the barycentric rational interpolation DQM.The state equation is approximated as a system of linear equations.Finally, the state transition matrix on adjacent milling periods is constructed, and the stability of the milling process is studied based on the Floquet theory.The accuracy and computational efficiency of the method are validated using two widely used milling models.The results of the simulation experiments show that the method proposed in this work has good accuracy and computational efficiency.The main characteristics of the method proposed in this article are listed as follows: 1.
Using the barycentric rational interpolation DQM can effectively improve the shortcomings of the classical DQM, avoiding the generation of the ill-conditioned matrix when there are a large number of discrete nodes, thereby improving the stability and accuracy of numerical calculations.

2.
The proposed method approximates the state equation of the milling system to an algebraic system of equations through interpolation and numerical differentiation techniques, which can quickly obtain the state transition matrix, and thus obtain the SLD of the system.

3.
The proposed milling stability analysis method is applicable to single-time-delay and multiple-time-delay milling systems, and is suitable for the machining conditions of both large and small radial depths of cutting.

Figure 1 .
Figure 1.Schematic of two-DOF milling model, (a) schematic of the milling model; (b) A-A direction view; (c) distribution of the cutter teeth; (d) the lag angle and tooth sweep angle.

Figure 1 .
Figure 1.Schematic of two-DOF milling model, (a) schematic of the milling model; (b) A-A direction view; (c) distribution of the cutter teeth; (d) the lag angle and tooth sweep angle.
.1.Polynomial Interpolation and the Barycentric Formula

−
Figure 2d, the time node i

Figure 2 .
Figure 2. Different locations of the time nodet i − T j .(a) t i − T j ∈ [−T, 0] and t i − T j ̸ = t k−T ; (b) t i − T j ∈ [−T, 0] and t i − T j = t k−T ; (c) t i − T j ∈ [0, T] and t i − T j ̸ = t k ; (d) t i − T j ∈ [0, T] and t i − T j = t k .

Figure 3a -
Figure3a-d display the SLD of the single-time-delay milling model with the radial immersion 1.The SLD obtained by the classical DQM is shown in Figure3a.As mentioned earlier, when the number of discrete nodes is large, the classical DQM produces an illconditioned differential matrix, resulting in erroneous stability prediction results.Figure3bshows the SLD obtained by the LDQM[33], from which it can be seen that when the parameter l is too large, the LDQM may still yield unstable results.In Figure3c,d, the SLD is obtained by the method proposed in Section 3. As shown in Figure3c,d, with the increase of the discrete parameter m, the SLD obtained by the method proposed in this work shows good consistency with the reference value, indicating that the method proposed in this work has good accuracy.Compared with the LDQM using l = 21 local nodes (Figure3b), the method proposed in this paper can obtain accurate results when the discrete parameter m is 60 (Figure3d).This indicates that the DQM based on barycentric rational interpolation effectively improves the shortcomings of the classical DQM that can easily produce an ill-conditioned matrix when the number of discrete nodes is large.Figure4a-dshow the SLD obtained by the proposed method with two other different radial immersions.In Figure4a,b, the radial immersion is 0.6, and the discrete parameters (m) are 20 and 30, respectively.In Figure4c,d, the discrete parameters (m) are 20 and 30 and the radial immersion is 0.2, which is a typical machining condition of the small radial depth of cutting.From Figure4, it can also be concluded that the method proposed in this work can obtain the accurate SLD to predict the stability of the milling process.From Figures3 and 4, it can be seen that when the cutter with symmetric pitch angles is used, the maximum value of the stability limit appears in the right area (near the spindle speed of 9200 rpm) of the SLD.We can choose the optimal axial depth of cutting to achieve the maximum material removal rate.At the same time, it can also be concluded that the method proposed in this work is not affected by the machining conditions and can obtain accurate SLD under both large and small radial depths of cutting.In order to further validate the efficiency of the method, the time required to calculate the SLD in numerical simulation is also compared with that of the 1st-SDM based on the above mentioned two-DOF milling model.The simulative calculations were performed in Matlab ® on a laptop computer (Intel ® Core TM i7-10870H CPU 2.20 GHz, 16.00 GB RAM).

Table 2 .
Values of the dynamics parameters of the single-time-delay milling model.

Table 3 .
Calculation time required for the 1st-SDM and the proposed method (unit: s).

Table 4 .
Values of the dynamics parameters of the multiple-time-delay milling model.