Next Article in Journal
On q-Hermite Polynomials with Three Variables: Recurrence Relations, q-Differential Equations, Summation and Operational Formulas
Previous Article in Journal
Load Calculation Method for Deep-Buried Layered Soft Rock Tunnel Based on Back-Analysis of Structural Deformation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
School of Construction Machinery, Chang’an University, Xi’an 710064, China
2
College of Mechanical & Electrical Engineering, Shaanxi University of Science & Technology, Xi’an 710021, China
3
School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China
4
School of Mechanics and Safety Engineering, Zhengzhou University, Zhengzhou 450001, China
*
Authors to whom correspondence should be addressed.
Symmetry 2024, 16(4), 384; https://doi.org/10.3390/sym16040384
Submission received: 26 February 2024 / Revised: 19 March 2024 / Accepted: 22 March 2024 / Published: 24 March 2024

Abstract

:
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.

1. 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 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 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.
From Table 1, it can be seen that the number of studies based on dynamic analysis methods is the highest, indicating that such methods have good application prospects.
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.

2. 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:
M y ¨ t + C y ˙ t + K y t = f t
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 ), φ s t and φ e x are the start and exit angles of the tooth, ϕ l a g z is the lag angle of the cutting edge at height z and ϕ s is the tooth sweep angle.
According to the geometric relationship described in Figure 1, ϕ l a g z and φ j t , z can be calculated through the following formula:
ϕ l a g z = 2 tan β D z
where D is the diameter of the milling cutter.
φ j t , z = φ j 0 , 0 + 2 π Ω 60 t ϕ l a g z
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:
f t = f c t + f d t
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:
f t = j = 1 N H j t y t y t T j
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:
H j t = h x x , j t h x y , j t h y x , j t h y y , j t
where
h x x , j t = z l , j t z u , j t K t cos φ j t , z sin φ j t , z K n sin φ j t , z sin φ j t , z d z h x y , j t = z l , j t z u , j t K t cos φ j t , z cos φ j t , z K n sin φ j t , z cos φ j t , z d z h y x , j t = z l , j t z u , j t + K t sin φ j t , z sin φ j t , z K n cos φ j t , z sin φ j t , z d z h y y , j t = z l , j t z u , j t + K t sin φ j t , z cos φ j t , z K n cos φ j t , z cos φ j t , z d z
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.
x ˙ t = A x t + j = 1 N B j t x t x t T j
where x t = y t y ˙ t T is the state vector and A and B j ( t ) can be expressed as
A = 0 I M 1 K M 1 C B j t = 0 0 M 1 H j t 0
Since H j t = H j t + T , it can be inferred that B j t = B j t + T . Hence, Equation (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 multiple-time-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.

3. Algorithm Derivation

3.1. 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.

3.1.1. Polynomial Interpolation and the Barycentric Formula

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:
p x = i = 0 n l i x f x i
where
l i x = k = 0 , k i n x x k k = 0 , k i n x i x k
and the Lagrange interpolation basis function l i has the following property:
l i x k = 1 , k = i 0 , otherwise k , i = 0 , , n
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 ,
w i = 1 k i x i x k , i = 0 , , n
Equation (10) can finally be expressed in the following form [38]:
p x = i = 0 n w i x x i f x i k = 0 n w k x x k
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.

3.1.2. Barycentric Rational Interpolation and Its Differentiation

The interpolation formula proposed by Floater and Hormann [34] can be expressed as follows:
r x = i = 0 n d λ i x p i x i = 0 n d λ i x
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:
λ i x = 1 i x x i x x i + 1 x x i + d
By introducing a new barycentric weight w j , Equation (15) can be written as the following barycentric form:
r x = j = 0 n w j x x j f x j k = 0 n w k x x k
where the barycentric weight w j is expressed as follows:
w j = i J j 1 i k = i , k j i + d 1 x j x k
where J j is an index set that satisfies J j : = i I : j d i j and I : = 0 , 1 , , n d .
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:
r x = j = 0 n q j x f x j
where
q j x = w j x x j k = 0 n w k x x k
Taking differentiation with respect to x on both sides of Equation (19) and substituting x = x i the following equation can be obtained:
r x i = j = 0 n q j x i f x j
Equation (21) provides a DQM based on the barycentric rational interpolation. By letting A i j = q j x i , the differentiation of r x at point x i i = 0 , , n can be expressed in the following matrix form:
r x 0 r x 1 r x n = A 00 A 01 A 0 n A 10 A 11 A 1 n A n 0 A n 1 A n n f x 0 f x 1 f x n
where A i j can be calculated using the following formula:
A i j = w j / w i x i x j i j k = 0 , k i n w k / w i x k x i i = j

3.2. 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:
x ˙ t i = A + j = 1 N B j t i x t i j = 1 N B j t i x t i T j
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:
x ˙ 0 x ˙ 1 x ˙ m = A + j = 1 N B j , 0 A + j = 1 N B j , 1 A + j = 1 N B j , m x 0 x 1 x m j = 1 N B j , 0 B j , 1 B j , m x 0 T j x 1 T j x m T j
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:
x ˙ 0 x ˙ 1 x ˙ m = W 00 W 01 W 0 m W 10 W 11 W 1 m W n 0 W n 1 W m m x 0 x 1 x m
where W p q 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 x t i i = 0 , 1 , , m or x t i T i = 0 , 1 , , m based on the location of t i T j . 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:
x t i T j = l = 0 m q l t i T j x t l T situation a x t k T situation b r = 0 m q r t i T j x t r situation c x t k situation d
Equation (27) can also be expressed in the form of vector multiplication as shown below.
x i T j = q 0 , i T j q 1 , i T j q m , i T j x 0 T x 1 T x m T T situation a 0 I 0 x 0 T x k T x m T T situation b q 0 , i T j q 1 , i T j q m , i T j x 0 x 1 x m T situation c 0 I 0 x 0 x k x m T situation d
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:
x 0 T j x s T j x s + 1 T j x m T j = q 0 , 0 T j q 1 , 0 T j q m , 0 T j q 0 , s T j q 1 , s T j q m , s T j 0 0 0 0 0 0 x 0 T x s T x s + 1 T x m T + 0 0 0 0 0 0 q 0 , s + 1 T j q 1 , s + 1 T j q m , s + 1 T j q 0 , m T j q 1 , m T j q m , m T j x 0 x s x s + 1 x m
Substituting Equations (26) and (29) into Equation (25) and carrying out a simple transformation yields the following expressions:
W 00 W 01 W 0 m W 10 W 11 W 1 m W m 0 W m 1 W m m x 0 x 1 x m = A + j = 1 N B j , 0 A + j = 1 N B j , 1 A + j = 1 N B j , m x 0 x 1 x m j = 1 N B j , 0 B j , 1 B j , m q 0 , 0 T j q 1 , 0 T j q m , 0 T j q 0 , s T j q 1 , s T j q m , s T j 0 0 0 0 0 0 x 0 T x s T x s + 1 T x m T j = 1 N B j , 0 B j , 1 B j , m 0 0 0 0 0 0 q 0 , s + 1 T j q 1 , s + 1 T j q m , s + 1 T j q 0 , m T j q 1 , m T j q m , m T j x 0 x s x s + 1 x m
Letting
W = W 00 W 01 W 0 m W 10 W 11 W 1 m W m 0 W m 1 W m m , T d , j = B j , 0 B j , 1 B j , m X T = x 0 x 1 x m T , X T = x 0 T x 1 T x m T T T c = A + j = 1 N B j , 0 A + j = 1 N B j , 1 A + j = 1 N B j , m T L , j = q 0 , 0 T j q 1 , 0 T j q m , 0 T j q 0 , s T j q 1 , s T j q m , s T j 0 0 0 0 0 0 , T R , j = 0 0 0 0 0 0 q 0 , s + 1 T j q 1 , s + 1 T j q m , s + 1 T j q 0 , m T j q 1 , m T j q m , m T j
Simplifying and rearranging Equation (30) yields the following equation:
X T = T c W j = 1 N T d , j T R , j 1 j = 1 N T d , j T L , j X T
Finally, the Floquet transition matrix Φ between adjacent milling periods can be obtained based on Equation (32) as follows:
Φ = T c W j = 1 N T d , j T R , j 1 j = 1 N T d , j T L , j
The stability of the milling process can be determined by the spectral radius of the Floquet transition matrix Φ according to the Floquet theory.

4. 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.

4.1. 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.
x ˙ t = A x t + B t x t x t T
where T is equal to the tooth passing period, and
A = 0 0 1 0 0 0 0 1 ω n 2 0 2 ζ ω n 0 0 ω n 2 0 2 ζ ω n ,   B t = 0 0 0 0 0 0 0 0 w h x x t m t w h x y t m t 0 0 w h y x t m t w h y y t m t 0 0
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 single-time-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 Figure 3 and Figure 4. The reference stability limits represented by the red line in Figure 3 and Figure 4 are obtained by the 1st-SDM with m = 200 .
Figure 3a–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 Figure 3a. As mentioned earlier, when the number of discrete nodes is large, the classical DQM produces an ill-conditioned 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 Figure 3 and Figure 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® CoreTM 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 Figure 3 and Figure 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.

4.2. 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.

5. 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:
  • 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.
  • 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.
  • 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.

Author Contributions

Conceptualization, Y.M.; Methodology, Y.M. and B.H.; Supervision, B.H. and S.H.; Validation, Y.M., X.R. and Z.Z.; Writing—original draft, Y.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number 12102239 and the Fundamental Research Funds for the Central Universities, CHD, grant number 300102252107 and the Natural Science Basic Research Program of Shaanxi Province of China, grant number 2021JQ-521.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Nomenclature

M Mass matrix of the milling system
K Stiffness matrix of the milling system
x t Vibration displacement of the tool in the x direction
f t Milling force vector
F y t The component of the milling force in the y direction
ϕ j The pitch angle between the tooth ( j 1 ) and the tooth ( j )
z u , j The upper bound of the cutting edge participating in cutting on the tooth ( j )
φ s t The start angle of the tooth
ϕ l a g z The lag angle of the cutting edge at height z
D The diameter of the milling cutter
f c The static component of the milling force
N The number of teeth
y t T j The vibration displacement vector of the previous tooth-passing period of the tooth ( j )
T The spindle period
K n The radial cutting force coefficient
y ˙ t Vibration velocity vector
B j ( t ) Time-varying parameter matrix corresponding to the tooth ( j ) in the state equation
l i ( x ) The Lagrange interpolation basis function
r x The barycentric rational interpolation function
w j The barycentric weight of the barycentric rational interpolation
q j x The barycentric rational interpolation basis function
W p q The corresponding weighted coefficient of the milling state vector
C Damping matrix of the milling system
y t Vibration displacement vector
y t Vibration displacement of the tool in the y direction
F x t The component of the milling force in the x direction
φ j t , z The angular position of the cutting edge on the tooth ( j ) with a height z
β The helix angle of the milling cutter
z l , j The lower bound of the cutting edge participating in cutting on the tooth ( j )
φ e x The exit angle of the tooth
ϕ s The tooth sweep angle
Ω The spindle speed
f d The dynamic component of the milling force
T j The time delay of the tooth ( j )
H j t The dynamic milling force coefficient matrix corresponding to the tooth ( j )
K t The tangential cutting force coefficient
x t The state vector
A Time-invariant parameter matrix in the state equation
p x The polynomial interpolation function
w i The barycentric weight of the barycentric Lagrange interpolation
λ i Blending function
J j The index set
A i j Weight coefficient of the barycentric rational interpolation differential quadrature method
Φ The Floquet transition matrix
DQMDifferential quadrature method
SLDStability lobe diagram
SDMSemi-discretization method
LDQMLocalized differential quadrature method
DDEDelay differential equation
TFEATemporal finite element analysis
FDMFull-discretization method
DOFDegree of freedom

References

  1. Munoa, J.; Beudaert, X.; Dombovari, Z.; Altintas, Y.; Budak, E.; Brecher, C.; Stepan, G. Chatter suppression techniques in metal cutting. CIRP Ann. 2016, 65, 785–808. [Google Scholar] [CrossRef]
  2. Ismail, F.; Soliman, E. A new method for the identification of stability lobes in machining. Int. J. Mach. Tools Manuf. 1997, 37, 763–774. [Google Scholar] [CrossRef]
  3. Peng, C.; Wang, L.; Liao, T.W. A new method for the prediction of chatter stability lobes based on dynamic cutting force simulation model and support vector machine. J. Sound Vib. 2015, 354, 118–131. [Google Scholar] [CrossRef]
  4. Quintana, G.; Ciurana, J.; Ferrer, I.; Rodriguez, C.A. Sound mapping for identification of stability lobe diagrams in milling processes. Int. J. Mach. Tools Manuf. 2009, 49, 203–211. [Google Scholar] [CrossRef]
  5. Campomanes, M.L.; Altintas, Y. An improved time domain simulation for dynamic milling at small radial immersions. J. Manuf. Sci. Eng. 2003, 125, 416–422. [Google Scholar] [CrossRef]
  6. Li, H.; Li, X.; Chen, X. A novel chatter stability criterion for the modelling and simulation of the dynamic milling process in the time domain. Int. J. Adv. Manuf. Technol. 2003, 22, 619–625. [Google Scholar] [CrossRef]
  7. Sims, N.D. The self-excitation damping ratio: A chatter criterion for time-domain milling simulations. J. Manuf. Sci. Eng. 2005, 127, 433–445. [Google Scholar] [CrossRef]
  8. Zhongqun, L.; Qiang, L. Solution and analysis of chatter stability for end milling in the time-domain. Chin. J. Aeronaut. 2008, 21, 169–178. [Google Scholar] [CrossRef]
  9. Altintaş, Y.; Budak, E. Analytical prediction of stability lobes in milling. CIRP Ann.-Manuf. Technol. 1995, 44, 357–362. [Google Scholar] [CrossRef]
  10. Budak, E.; Altintas, Y. Analytical prediction of chatter stability in milling—Part I: General formulation. J. Dyn. Syst. Meas. Control 1998, 120, 22–30. [Google Scholar] [CrossRef]
  11. Budak, E.; Altintas, Y. Analytical prediction of chatter stability in milling—Part II: Application of the general formulation to common milling systems. J. Dyn. Syst. Meas. Control 1998, 120, 31–36. [Google Scholar] [CrossRef]
  12. Bayly, P.; Halley, J.; Mann, B.P.; Davies, M. Stability of interrupted cutting by temporal finite element analysis. J. Manuf. Sci. Eng. 2003, 125, 220–225. [Google Scholar] [CrossRef]
  13. Mann, B.P.; Young, K.A.; Schmitz, T.L.; Dilley, D.N. Simultaneous stability and surface location error predictions in milling. J. Manuf. Sci. Eng. 2005, 127, 446–453. [Google Scholar] [CrossRef]
  14. Insperger, T.; Stépán, G. Semi-discretization method for delayed systems. Int. J. Numer. Methods Eng. 2002, 55, 503–518. [Google Scholar] [CrossRef]
  15. Insperger, T.; Stépán, G. Updated semi-discretization method for periodic delay-differential equations with discrete delay. Int. J. Numer. Methods Eng. 2004, 61, 117–141. [Google Scholar] [CrossRef]
  16. Insperger, T.; Stépán, G.; Turi, J. On the higher-order semi-discretizations for periodic delayed systems. J. Sound Vib. 2008, 313, 334–341. [Google Scholar] [CrossRef]
  17. Ding, H.; Ding, Y.; Zhu, L. On time-domain methods for milling stability analysis. Chin. Sci. Bull. 2012, 57, 4336–4345. [Google Scholar] [CrossRef]
  18. Ding, Y.; Zhu, L.; Zhang, X.; Ding, H. A full-discretization method for prediction of milling stability. Int. J. Mach. Tools Manuf. 2010, 50, 502–509. [Google Scholar] [CrossRef]
  19. Butcher, E.A.; Ma, H.; Bueler, E.; Averina, V.; Szabo, Z. Stability of linear time-periodic delay-differential equations via Chebyshev polynomials. Int. J. Numer. Methods Eng. 2004, 59, 895–922. [Google Scholar] [CrossRef]
  20. Butcher, E.A.; Bobrenkov, O.A.; Bueler, E.; Nindujarla, P. Analysis of milling stability by the Chebyshev collocation method: Algorithm and optimal stable immersion levels. J. Comput. Nonlinear Dyn. 2009, 4, 031003. [Google Scholar] [CrossRef]
  21. Lin, L.; He, M.; Wang, Q.; Deng, C. Chatter Stability Prediction and Process Parameters’ Optimization of Milling Considering Uncertain Tool Information. Symmetry 2021, 13, 1071. [Google Scholar] [CrossRef]
  22. Ding, Y.; Zhu, L.; Zhang, X.; Ding, H. Numerical integration method for prediction of milling stability. J. Manuf. Sci. Eng. 2011, 133, 031005. [Google Scholar] [CrossRef]
  23. Niu, J.; Ding, Y.; Zhu, L.; Ding, H. Runge–Kutta methods for a semi-analytical prediction of milling stability. Nonlinear Dyn. 2014, 76, 289–304. [Google Scholar] [CrossRef]
  24. Zhang, Z.; Li, H.; Meng, G.; Liu, C. A novel approach for the prediction of the milling stability based on the Simpson method. Int. J. Mach. Tools Manuf. 2015, 99, 43–47. [Google Scholar] [CrossRef]
  25. Mei, Y.; Mo, R.; Sun, H.; He, B.; Wan, N. Stability prediction in milling based on linear multistep method. Int. J. Adv. Manuf. Technol. 2019, 105, 2677–2688. [Google Scholar] [CrossRef]
  26. Mei, Y.; Mo, R.; Sun, H.; He, B.; Bu, K. Stability analysis of milling process with multiple delays. Appl. Sci. 2020, 10, 3646. [Google Scholar] [CrossRef]
  27. Yang, Y.; Yuan, J.-W.; Tie, D.; Wan, M.; Zhang, W.-H. An efficient and accurate chatter prediction method of milling processes with a transition matrix reduction scheme. Mech. Syst. Signal Process. 2023, 182, 109535. [Google Scholar] [CrossRef]
  28. Zhang, X.; Xiong, C.; Ding, Y.; Ding, H. Prediction of chatter stability in high speed milling using the numerical differentiation method. Int. J. Adv. Manuf. Technol. 2017, 89, 2535–2544. [Google Scholar] [CrossRef]
  29. Ding, Y.; Zhu, L.; Zhang, X.; Ding, H. Stability analysis of milling via the differential quadrature method. J. Manuf. Sci. Eng. 2013, 135, 044502. [Google Scholar] [CrossRef]
  30. Bellman, R.; Kashef, B.; Casti, J. Differential quadrature: A technique for the rapid solution of nonlinear partial differential equations. J. Comput. Phys. 1972, 10, 40–52. [Google Scholar] [CrossRef]
  31. Bellman, R.; Casti, J. Differential quadrature and long-term integration. J. Math. Anal. Appl. 1971, 34, 235–238. [Google Scholar] [CrossRef]
  32. Zong, Z.; Lam, K. A localized differential quadrature (LDQ) method and its application to the 2D wave equation. Comput. Mech. 2002, 29, 382–391. [Google Scholar] [CrossRef]
  33. Mei, Y.; He, B.; He, S.; Ren, X. Stability Analysis in Milling Based on the Localized Differential Quadrature Method. Micromachines 2023, 15, 54. [Google Scholar] [CrossRef]
  34. Floater, M.S.; Hormann, K. Barycentric rational interpolation with no poles and high rates of approximation. Numer. Math. 2007, 107, 315–331. [Google Scholar] [CrossRef]
  35. Gradišek, J.; Kalveram, M.; Insperger, T.; Weinert, K.; Stépán, G.; Govekar, E.; Grabec, I. On stability prediction for milling. Int. J. Mach. Tools Manuf. 2005, 45, 769–781. [Google Scholar] [CrossRef]
  36. Altintas, Y. Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design; Cambridge University Press: Cambridge, UK, 2012. [Google Scholar]
  37. Sauer, T. Numerical Analysis, 2nd ed.; Pearson Education, Inc.: Boston, MA, USA, 2012. [Google Scholar]
  38. Berrut, J.-P.; Trefethen, L.N. Barycentric lagrange interpolation. SIAM Rev. 2004, 46, 501–517. [Google Scholar] [CrossRef]
  39. Schneider, C.; Werner, W. Some new aspects of rational interpolation. Math. Comput. 1986, 47, 285–299. [Google Scholar] [CrossRef]
  40. Altıntas, Y.; Engin, S.; Budak, E. Analytical stability prediction and design of variable pitch cutters. J. Manuf. Sci. Eng. 1999, 121, 173–178. [Google Scholar] [CrossRef]
  41. Tao, J.; Qin, C.; Liu, C. Milling stability prediction with multiple delays via the extended Adams-Moulton-based method. Math. Probl. Eng. 2017, 2017, 7898369. [Google Scholar] [CrossRef]
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. 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.
Symmetry 16 00384 g001
Figure 2. Different locations of the time node t 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 2. Different locations of the time node t 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 .
Symmetry 16 00384 g002
Figure 3. SLD of the single-time-delay milling model obtained by the classical DQM, the LDQM, and the method proposed in this work. (a) classical DQM (m = 60); (b) LDQM (m = 60, l = 21); (c) the proposed method (m = 40, d = 4); (d) the proposed method (m = 60, d = 4).
Figure 3. SLD of the single-time-delay milling model obtained by the classical DQM, the LDQM, and the method proposed in this work. (a) classical DQM (m = 60); (b) LDQM (m = 60, l = 21); (c) the proposed method (m = 40, d = 4); (d) the proposed method (m = 60, d = 4).
Symmetry 16 00384 g003
Figure 4. SLD of the single-time-delay milling model obtained by the method proposed in this work. (a) radial immersion 0.6 (m = 20, d = 4); (b) radial immersion 0.6 (m = 30, d = 4); (c) radial immersion 0.2 (m = 20, d = 4); (d) radial immersion 0.2 (m = 30, d = 4).
Figure 4. SLD of the single-time-delay milling model obtained by the method proposed in this work. (a) radial immersion 0.6 (m = 20, d = 4); (b) radial immersion 0.6 (m = 30, d = 4); (c) radial immersion 0.2 (m = 20, d = 4); (d) radial immersion 0.2 (m = 30, d = 4).
Symmetry 16 00384 g004
Figure 5. SLD of the multiple-time-delay milling model obtained by the method proposed in this work. (a) radial immersion 1.0 (m = 100, d = 5); (b) radial immersion 0.6 (m = 100, d = 5); (c) radial immersion 0.2 (m = 100, d = 5).
Figure 5. SLD of the multiple-time-delay milling model obtained by the method proposed in this work. (a) radial immersion 1.0 (m = 100, d = 5); (b) radial immersion 0.6 (m = 100, d = 5); (c) radial immersion 0.2 (m = 100, d = 5).
Symmetry 16 00384 g005
Table 1. The most commonly used milling stability analysis methods.
Table 1. The most commonly used milling stability analysis methods.
Stability Analysis MethodsReferences
Experimental-based methods[2,3,4]
Time-domain simulation methods[5,6,7,8]
Dynamic analysis methods[9,10,11,12,13,14,15,16,18,19,20,22,23,24,25,26,27,28,29,33]
Table 2. Values of the dynamics parameters of the single-time-delay milling model.
Table 2. Values of the dynamics parameters of the single-time-delay milling model.
mt
(kg)
ζωn
(rad/s)
Kt
(N/m2)
Kn
(N/m2)
N
0.033930.01157936 × 1082 × 1082
Table 3. Calculation time required for the 1st-SDM and the proposed method (unit: s).
Table 3. Calculation time required for the 1st-SDM and the proposed method (unit: s).
m = 60 Radial Immersion: 1Radial Immersion: 0.6Radial Immersion: 0.2
1st-SDM460.98321.87280.20
The   proposed   method   ( d = 3 )56.9654.2652.92
The   proposed   method   ( d = 4 )56.4555.1753.10
Table 4. Values of the dynamics parameters of the multiple-time-delay milling model.
Table 4. Values of the dynamics parameters of the multiple-time-delay milling model.
mt
(kg)
ζωn
(rad/s)
Kt
(N/m2)
Kn
(N/m2)
m t x = 1.4986 ζ x = 0.0558 ω n x = 563.55 × 2 π 6.79 × 10 8 2.56 × 10 8
m t y = 1.199 ζ y = 0.025 ω n y = 516.27 × 2 π
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Mei, Y.; He, B.; He, S.; Ren, X.; Zhang, Z. Stability Analysis of Milling Based on the Barycentric Rational Interpolation Differential Quadrature Method. Symmetry 2024, 16, 384. https://doi.org/10.3390/sym16040384

AMA Style

Mei Y, He B, He S, Ren X, Zhang Z. Stability Analysis of Milling Based on the Barycentric Rational Interpolation Differential Quadrature Method. Symmetry. 2024; 16(4):384. https://doi.org/10.3390/sym16040384

Chicago/Turabian Style

Mei, Yonggang, Bingbing He, Shangwen He, Xin Ren, and Zeqi Zhang. 2024. "Stability Analysis of Milling Based on the Barycentric Rational Interpolation Differential Quadrature Method" Symmetry 16, no. 4: 384. https://doi.org/10.3390/sym16040384

APA Style

Mei, Y., He, B., He, S., Ren, X., & Zhang, Z. (2024). Stability Analysis of Milling Based on the Barycentric Rational Interpolation Differential Quadrature Method. Symmetry, 16(4), 384. https://doi.org/10.3390/sym16040384

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop