Next Article in Journal
New Era of Electroceuticals: Clinically Driven Smart Implantable Electronic Devices Moving towards Precision Therapy
Next Article in Special Issue
Numerical Modeling of the Effect of Cutting-Edge Radius on Cutting Force and Stress Concentration during Machining
Previous Article in Journal
Mixing Improvement in a T-Shaped Micro-Junction through Small Rectangular Cavities
Previous Article in Special Issue
Influence of Some Microchanges Generated by Different Processing Methods on Selected Tribological Characteristics
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Novel Updated Full-Discretization Method for Prediction of Milling Stability

1
School of Mechanical and Power Engineering, Henan Polytechnic University, Jiaozuo 454000, China
2
Key Laboratory of Contemporary Design and Integrated Manufacturing Technology, Ministry of Education, Northwestern Polytechnical University, Xi’an 710072, China
3
College of Computer Science and Technology, Henan Polytechnic University, Jiaozuo 454000, China
*
Author to whom correspondence should be addressed.
Micromachines 2022, 13(2), 160; https://doi.org/10.3390/mi13020160
Submission received: 30 December 2021 / Revised: 13 January 2022 / Accepted: 19 January 2022 / Published: 21 January 2022
(This article belongs to the Special Issue Advances in Ultra-Precision Machining Technology and Applications)

Abstract

:
This paper presents an updated full-discretization method for milling stability prediction based on cubic spline interpolation. First, the mathematical model of the time-delay milling system considering regenerative chatter is represented by a dynamic delay differential equation. Then, in a single tooth passing period, the time is divided into a finite time intervals, the state item and the time-delay item are approximated in each time interval by cubic spline interpolation and third-order Newton interpolation, respectively. Afterward, a transition matrix is constructed to represent the transfer relationship of the teeth in a period. Finally, based on Floquet theory, the milling stability lobes can be obtained. Meanwhile, in order to improve computational efficiency, an optimized method is proposed based on the traditional algorithm and the proposed method has high precision without losing high efficiency. Finally, several milling experiments are conducted to verify the accuracy of the proposed method, and the results show that the predicted results agree well with the experimental results.

1. Introduction

Chatter is a serious problem in the milling of a thin-walled workpiece such as aeroengine blades, casings, impellers, blisks etc., which not only reduces the surface quality and production efficiency but also shortens the life of machine spindles and cutters. Therefore, it is necessary to investigate chatter in the milling of thin-walled workpieces for obtaining chatter-free operations.
Up to now, chatter has been studied by considering the time-varying milling process system by several researchers, including T. Insperger [1,2,3], Chandiramani [4], and Eksioglu [5] et al. From these studies, we found that in the milling process, considerable valuable results on chatter are investigated by a mathematical model of a time-periodic delay differential equation (DDE), and stability lobes diagram can be used to obtain the chatter-free process parameters more accurately.
To investigate machine stability considering regeneration chatter, a considerable number of methods, including analytical methods, numerical methods, and experimental methods, have been proposed to predict the stability lobes diagram (SLD), e.g., the analytical methods [6], the temporal finite element methods [7], the semi-discretization and full-discretization methods [8,9], the time domain numerical simulation methods [10], and the experimental methods [11]. Among the methods mentioned above, the semi-discretization methods and full-discretization methods are widely used due to their efficiency and accuracy.
In recent years, Altintas [12] proposed a zero-order analytical (ZOA) method by the Fourier series average method of dynamic milling force coefficient. Then, Merdol [13] transformed the low radial immersion milling dynamics into an eigenvalue problem by considering the tooth spacing angle and tooth passing frequencies for accurately predicting the stability lobes diagram. The semi-discrete cycloid method based on the nonlinear cutting force milling model was presented by Faassen [14] for investigating the stability of a milling system, which was proved effectively.
Furthermore, to obtain a high-precision stability lobe diagram, the numerical methods, including the numerical integration method [15,16], Runge–Kutta-based discretization method [17], and precise integration method [18] were developed. Shortly after that, a semi-discrete method (SDM) was proposed by Insperger [19] and then, a full-discretization method (FDM) was presented by Ding et al. [20]. Later on, a second-order FDM method [21], third-order FDM method [22], high-order FDM method [23], Hermite interpolation FDM [24], and the update FDM [25] are successively proposed. It is proved that the accuracy and efficiency of these methods were improved to some extent. However, these extended methods have more complex algorithm structures, which will take more calculation time and obtain the low convergence accuracy to a certain degree.
Therefore, to obtain a higher convergence rate and computational efficiency, a novel update FDM based on Spline–Newtons interpolation is proposed in this paper. The most important difference of the proposed method compared with the existing methods is that the cubic spline interpolation method was utilized to handle the state item and the third-order Newton interpolation method was used to approximate the time-delay item. The remainder of the paper is organized as follows. In Section 2, a systematic mathematical model and algorithm are described in detail. In Section 3, the rate of convergence estimates of the proposed method is calculated compared with some existing method. In Section 4 and Section 5, two benchmark examples for a one degrees of freedom (DOF) milling model are given to illustrate the accuracy and efficiency of the proposed approach. Some verified experiments are conducted and analyzed in Section 6. Finally, some conclusions are presented.

2. Systematic Mathematical Model and Algorithm

For a conventional milling process, a schematic diagram of milling a thin-walled section while considering regenerative chatter is given in Figure 1. Without loss of generality, the dynamic model of the milling process system considering the regenerative effect can be expressed by a n-dimensional linear time periodic system with a single discrete time delay as follows:
X ˙ ( t ) = A 0 X ( t ) + A ( t ) [ X ( t ) X ( t T ) ]
where A0 is a constant matrix, A(t) is a time-periodic matrix, and A(t) = A(t + T), T is the time period which equals to the time delay, and X(t) is the relative displacement between the cutter and workpiece. In order to solve Equation (1), time period T can be equally divided into m small-time intervals, and T = mh, where m is an integer and h is the range of each interval. Then, the dynamic response of Equation (1) can be obtained by a direct integration method on each time interval kht ≤ (k + 1)h:
X ( t ) = e A 0 ( t k h ) X ( k h ) + k h t { e A 0 ( t ε ) [ A ( ε ) X ( ε ) A ( ε ) X ( ε T ) ] } d ε .
Let X(kh) = Xk with k = 0, 1,…, m, when t = (k + 1)h, Equation (2) can be equivalently converted to the following form:
X k + 1 = e A 0 h X k + 0 h e A 0 ε [ A ( k h + h ε ) X ( k h + h ε ) A ( k h + h ε ) X ( k h + h ε T ) ] d ε
Next, the integral term in Equation (3) should be handled. The state item X(kh + h − ε) can be approximately represented by cubic spline interpolation using Xk+1, XkXk−1, Xk−2. In addition, two other constraints X ˙ ( k h 2 h ) and X ˙ ( k h + h ) are used in cubic spline interpolation. Namely, let Ak stands for A(kh):
{ X ˙ ( k h 2 h + 0 ) = A 0 X k 2 + A k 2 ( X k 2 X k 2 m ) X ˙ ( k h + h 0 ) = A 0 X k + 1 + A k + 1 ( X k + 1 X k + 1 m ) .  
At the time interval [tk, tk+1], the state item X(kh + h − ε) can be approximated by cubic spline interpolation, resulting in:
X ( k h + h ε ) = μ 1 X k + 1 + μ 2 X k + μ 3 X k 1 + μ 4 X k 2
where
μ 1 = 11 h A 0 + 18 I 15 h 3 ε 3 + 26 h A 0 33 I 15 h 2 ε 2 A 0 ε + I
μ 2 = 9 I 5 h 3 ε 3 + 14 I 5 h 2 ε 2
μ 3 = 4 I 5 h 3 ε 3 4 I 5 h 2 ε 2
μ 4 = h A 0 3 I 15 h 3 ε 3 + h A 0 + 3 I 15 h 2 ε 2 .  
The time delay item X(kh + h − ε T) in Equation (3) can be approximately expressed at four points Xk−m+3, Xk−m+2, Xk−m+1, Xkm by the third-order Newton interpolation method as follows:
X ( k h + h ε T ) = λ 1 X k m + λ 2 X k m + 1 + λ 3 X k m + 2 + λ 4 X k m + 3
where
λ 1 = ε 3 6 h 3 + ε 2 2 h 2 + ε 3 h
λ 3 = ε 3 2 h 3 + ε 2 2 h 2 ε h
λ 3 = ε 3 2 h 3 + ε 2 2 h 2 ε h
λ 4 = ε 3 6 h 3 + ε 6 h .  
Subsequently, the time-periodic item A(kh + hε) in Equation (3) can be expressed by linear interpolation which using points A(kh) and A(kh + h), and:
A ( k h + h ε ) = A u + A v ε
where
A u = A k + 1
A v = ( A k A k + 1 ) / h .  
Then, substituting Equation (5), Equation (10), and Equation (15) into Equation (3) leads to the interpolated item X(kh + h ε), X(kh + hε − T) and A(kh + h ε) are taken into Equation (3), and the DDE is approximated by an ordinary differential equation (ODE), which can be simplified as follows:
P k + 1 X k + 1 = P k X k + P k 1 X k 1 + P k 2 X k 2 + P k m + 3 X k m + 3 + P k m + 2 X k m + 2 + P k m + 1 X k m + 1 + P k m X k m
where
P k + 1 = I ( a 1 A u + a 2 A v )
P k = Φ 0 + ( b 1 A u + b 2 A v )
P k 1 = c 1 A u + c 2 A v
P k 2 = d 1 A u + d 2 A v
P k m + 3 = e 1 A u + e 2 A v  
P k m + 2 = f 1 A u + f 2 A v
P k m + 1 = g 1 A u + g 2 A v
P k m = h 1 A u + h 2 A v .  
Define:
Φ 0 = e A 0 h
Φ 1 = 0 h e A 0 ε d ε = A 0 1 ( Φ 0 I )
Φ 2 = 0 h ε e A 0 ε d ε = A 0 1 ( h Φ 0 Φ 1 )
Φ 3 = 0 h ε 2 e A 0 ε d ε = A 0 1 ( h 2 Φ 0 2 Φ 2 )
Φ 4 = 0 h ε 3 e A 0 ε d ε = A 0 1 ( h 3 Φ 0 3 Φ 3 )
Φ 5 = 0 h ε 4 e A 0 ε d ε = A 0 1 ( h 4 Φ 0 4 Φ 4 ) .  
In addition, the coefficients in Equations (19)–(26) can be expressed as:
a 1 = 11 h A 0 + 18 I 15 h 3 Φ 4 + 26 h A 0 33 I 15 h 2 Φ 3 A 0 Φ 2 + Φ 1
a 2 = 11 h A 0 + 18 I 15 h 3 Φ 5 + 26 h A 0 33 I 15 h 2 Φ 4 A 0 Φ 3 + Φ 2
b 1 = 9 5 h 3 Φ 4 + 14 5 h 2 Φ 3
b 2 = 9 5 h 3 Φ 5 + 14 5 h 2 Φ 4
c 1 = 4 5 h 3 Φ 4 4 5 h 2 Φ 3
c 2 = 4 5 h 3 Φ 5 4 5 h 2 Φ 4
d 1 = h A 0 3 I 15 h 3 Φ 4 + h A 0 + 3 I 15 h 2 Φ 3
d 2 = h A 0 3 I 15 h 3 Φ 5 + h A 0 + 3 I 15 h 2 Φ 4
e 1 = 1 6 h 3 Φ 4 1 6 h Φ 2
e 2 = 1 6 h 3 Φ 5 1 6 h Φ 3
f 1 = 1 2 h 3 Φ 4 1 2 h 2 Φ 3 + 1 h Φ 2
f 2 = 1 2 h 3 Φ 5 1 2 h 2 Φ 4 + 1 h Φ 3
g 1 = 1 2 h 3 Φ 4 + 1 h 2 Φ 3 1 2 h Φ 2 Φ 1
g 2 = 1 2 h 3 Φ 5 + 1 h 2 Φ 4 1 2 h Φ 3 Φ 2
h 1 = 1 6 h 3 Φ 4 1 2 h 2 Φ 3 1 3 h Φ 2
h 2 = 1 6 h 3 Φ 5 1 2 h 2 Φ 4 1 3 h Φ 3 .  
Then, if the matrix Pk+1 in Equation (18) is nonsingular, Equation (18) can be given by:
X k + 1 = P k + 1 1 P k X k + P k + 1 1 P k 1 X k 1 + P k + 1 1 P k 2 X k 2 + P k + 1 1 P k m + 3 X k m + 3 + P k + 1 1 P k m + 2 X k m + 2 + P k + 1 1 P k m + 1 X k m + 1 + P k + 1 1 P k m X k m
According to Equation (49), a discrete map could be defined as:
Q N k + m = R N k
where
N k = [ X k m , X k m + 1 , , X k 1 , X k ] T .  
In addition, Q and R can be expressed as:
R = [ 0 0 0 0 0 0 0 0 0 I P k , m P k , m 1 P k , m 2 P k , m 3 0 0 0 P k , 2 P k , 1 0 0 P k + 1 , m P k + 1 , m 1 P k + 1 , m 2 P k + 1 , m 3 0 0 0 P k + 1 , 2 0 0 0 P k + 2 , m   P k + 2 , m 1 P k + 2 , m 2 P k + 2 , m 3 0 0 0 0 0 0 0 0 0 0 P k + m 3 , 2 P k + m 3 , 1 P k + m 3 , 0 P k + m 3 , 1 0 0 0 0 0 0 0 P k + m 2 , 2 P k + m 2 , 1 P k + m 2 , 0 0 0 0 0 0 0 0 0 P k + m 1 , 2 P k + m 1 , 1 ]
Q = [ I 0 0 0 0 0 0 0 0 0 P k , 0 P k , 1 0 0 0 0 0 0 0 0 P k + 1 , 1 P k + 1 , 0 P k + 1 , 1 0 0 0 0 0 0 0 P k + 2 , 2 P k + 2 , 1 P k + 2 , 0 P k + 2 , 1 0 0 0 0 0 0 0 0 0 0 P k + m 3 , 2 P k + m 3 , 1 P k + m 3 , 0 P k + m 3 , 1 0 0 0 0 0 0 0 P k + m 2 , 2 P k + m 2 , 1 P k + m 2 , 0 P k + m 2 , 1 0 0 P k + m 1 , m 2 P k + m 1 , m 3 0 0 0 P k + m 1 , 2 P k + m 1 , 1 P k + m 1 , 0 P k + m 1 , 1 ]
It is clear that Q and R are both a (2m + 2) × (2m + 2) dimensional matrix. Therefore, the transition matrix V in a single tooth passing period can be defined as:
V = Q 1 R .  
Now, according to Floquet theory [26], the stability of the system can be determined by judging whether the modulus of the eigenvalues of the transition matrixes are less than 1 or not. If not, the system is unstable, otherwise, it will be stable.
Remark: If Pk+1 is singular, the processing method in Ref. [20] can be utilized. From Equations (52) and (53), it can be seen that 8m variables need to be calculated complicatedly compared with the first-order and second-order FDM, which leads to the increase of calculating time. To enhance the calculation efficiency, the traditional algorithms compressed the 2m + 2 dimensional matrix into a m + 2 dimensional matrix, and calculated the eigenvalues of the transition matrix in one period in the whole region, then, the stability boundary is drawn. However, a novel algorithm is proposed to obviously improve computational efficiency. It is well known that the machining process is stable below the boundary of the SLD and is unstable on the upper boundary of the SLD, while on the stable boundary the eigenvalue of the transition matrix is 1, represented by the spindle speed and depth of cut. According to the constraints of the spindle speed and depth of cut, the modulus of the transition matrix eigenvalues are calculated and judged with 1. If the value is more than 1, the algorithm stops, otherwise it continues, which is only the modulus of transition matrix eigenvalues that are less than 1 and calculated for the stability boundary, which can greatly improve the computational efficiency by reducing the calculation of the eigenvalues in the unstable region.

3. Rate of Convergence Estimates

To verify the fast convergence rate of the proposed method, a one DOF dynamic milling system is selected, and the proposed method is compared to the 0th SDM, 1st FDM, 2nd FDM, 3rd FDM, and Hermite FDM.
The rate of convergence can be clearly determined by the local discretization error (LDE). As stated in the literature [19], for the 0th SDM, the LDE is O(h2). The LDE of FDM with the 1st, 2nd, 3rd, and Hermite are O(h2), O(h3), O(h4), and O(h3) [24,27], respectively. For the proposed method, the LDE is O(h4).
All operations are from the same computer environment: Matlab 2018b, Inter(R) Core(TM) i5-4210H CPU @ 2.90 GHz 2.90 GHz. The milling system parameters are derived from the Ding [20]: The damping ratio ζ, model mass mt, and natural frequency wn are 0.011, 0.03993 kg, and 1844 Hz, respectively. The cutter has two flutes. The cutting force coefficients are Ktc = 6 × 108 and Krc = 2 × 108.
The spindle speed n is selected as 5000 rpm, and the axial depths of cut ap is selected as 0.1 mm, 0.2 mm, 0.5 mm, and 0.80 mm, respectively. The exact eigenvalues | μ 0 | corresponding to different axial depths of cut are 0.7368, 0.8192, 1.0726, and 1.2880, respectively. The LDE can be known as the absolute value of difference between the current eigenvalue | μ | and exact eigenvalue | μ 0 | .
As shown in Figure 2, the LDE of the 0th SDM, 1st FDM, 2nd FDM, 3rd FDM, Hermite FDM, and the proposed method are analyzed. From Figure 2, it can be clearly seen that the proposed method has a faster convergence rate. It should be mentioned that the proposed method is able to converge to a sufficient accuracy when the discrete number m is small.

4. Computational Accuracy Analysis

Under low and high immersion ratio conditions, the effectiveness of the proposed method is verified in terms of both computational efficiency and accuracy of milling stability prediction by comparing with other methods. The modal parameter selection is consistent with Section 3 of this paper. Equation (1) is the dynamic state-space model of the milling system, where constant matrix A0 and time-periodic matrix A(t) can be express as:
A 0 = [ ζ w n   1 / m t m t ( ζ w n ) 2 m t w n 2 ζ w n ]
A ( t ) = [ 0 0 a p h ( t ) 0 ]
h ( t ) = j = 1 N g ( φ j ( t ) ) sin ( φ j ( t ) ) ( K t c cos ( φ j ( t ) ) + K r c sin ( φ j ( t ) ) )
φ j ( t ) = 2 π n 60 t + ( j 1 ) 2 π N
g ( φ j ( t ) ) = { 1 φ s t < φ j ( t ) < φ e x 0 o t h e r w i s e
where φj(t) is the angular position of the j-th cutter tooth, and φst and φex are the starting and exiting edge positions of the tool in contact with the workpiece, respectively. For down-milling, φst = arccos(2a/D − 1) and φex = π; for up-milling, φst = 0 and φex = arccos(1 − 2a/D), where a/D is the radial immersion ratio.
For a/D = 1, all methods are calculated over a 200 × 100-sized grid of parameters under the condition of n = 5000–10,000 rpm and ap = 0–4 mm. For a/D = 0.1, all methods are calculated over a 200 × 100-sized grid of parameters under the condition of n = 5000–12,000 rpm and ap = 0–5 mm. The SLDs for a/D = 1 and a/D = 0.1 are shown in Figure 3 and Figure 4. The red curves in Figure 3 and Figure 4 are the reference curves, calculated by the Hermite FDM at the discrete number m = 100. It can be seen that the proposed method has a high accuracy both at a low and high radial immersion ratio. In particular, the 0th SDM, 1st FDM, 2nd FDM, 3rd FDM, and Hermite FDM all have large errors with the reference curve when m = 20 at the a/D = 1. However, the proposed method almost coincides with the reference curve. For a/D = 0.1, the 1st FDM agrees best with the reference curve when m = 20, and the proposed method agrees with the reference curve immediately as m increases.

5. Computational Efficiency Analysis

To verify the computational efficiency of the proposed method, the time required for the computation of the FDM in Section 4 for different discrete numbers m is discussed. The time required for the calculation is shown in Figure 5. From Figure 5, it can be seen that the proposed method has a faster computational efficiency compared to other methods. When a/D = 1, the proposed method saves an average of 69.2%, 73.3%,75.4%, and 66.7% of time compared to the 1st FDM, 2nd FDM, 3rd FDM, and Hermite FDM, respectively. When a/D = 0.1, the proposed method saves an average of 53.3%, 58.8%, 63.8%, and 47.5% of time compared to the 1st FDM, 2nd FDM, 3rd FDM, and Hermite FDM, respectively. It can be seen that the proposed method has a higher computational linear efficiency when a/D = 1. The main reason is that when a/D = 1, the contact time between the cutter and workpiece is at a maximum, while the transfer matrix Equation (54) needs to be calculated multiple times, and the proposed method saves the calculation time required for the stability region.

6. Verification

To verify the effectiveness of the proposed method in the milling of the thin-walled plate, some experiments are conducted in this section. The dimension of the plate used in modal test and machining experiments is 80 × 40 × 3 mm, and all experiments were carried out on the three-axis milling center (VMC-850E), which is shown in Figure 6. The material properties of the workpiece and the cutter parameters are given in Table 1.

6.1. Cutting Force Coefficients Calibration

For cutting force coefficients calibrated, as is known to all, when full-immersion milling (slot milling) is used, the average milling forces are expressed as:
{ F ¯ x = N a 4 K r c f N a π K r e F ¯ y = N a 4 K t c f + N a π K t e .
Then, for slot milling, five groups of full-immersion milling experiments were carried out. The machining parameters are the spindle speed 1000 rpm, axial depth of cut at 0.5 mm, and feed rate at 40 mm/min, 80 mm/min, 120 mm/min, 160 mm/min, and 200 mm/min.
Therefore, the average milling forces at each feed rate are measured by Kistler9257B, and the cutting-edge components are estimated by a linear regression of the accumulated data. Next, based on the literature [28], the cutting force coefficients are evaluated as Ktc = 1120.8 N/mm2, Krc = 2285.6 N/mm2, Kte = 9.16 N/mm, and Kre = 13.21 N/mm.

6.2. Modal Parameters Identification

An impact experiment is conducted for obtaining the modal parameters of the thin-walled workpiece. The modal parameters of the milling system are obtained by an acquisition instrument DH5981, acceleration sensors (Ref. sensitivity 10.25 mV/g), and modal hammer (500 N).
In tests, for a different measured position on the workpiece, the dynamic response is different. Therefore, considering the clamping constraints, the impact measured points 1, 2, and 3 distributed on the thin-walled plate are shown in Figure 7a. Point 1 and point 3 are symmetric with respect to point 2, and point 2 locates the middle of the thin-walled plate edge. Next, all the vibration responses on the different measured points are obtained, and according to the experimental results, we found that the vibration response at point 1 is the same as that at point 3. In addition, considering the unstable state in the cut-in and cut-out region, the representative point 2 are chosen for measuring responses, as shown in Figure 7b,c, and the experimental setup is shown in Figure 6. Therefore, the modal parameters in point 2 are identified and listed in Table 2.

6.3. Machining Tests

In this section, in order to validate the accuracy of the proposed approach for quickly and accurately predicting the stability of the milling system, the four degree of freedom in the X and Y direction for the milling cutter and workpiece in the X and Y direction for thin-walled section are considered. According to the proposed method, the stability lobe diagram with the discrete number m = 40 at the a/D = 0.1 is calculated, and the milling parameters are determined based on the stability lobe diagram as shown in Figure 8. The milling parameters in points A(n = 1500 rpm, ap = 0.2 mm) and C(n = 2500 rpm, ap = 0.4 mm) are stable parameters, while the points B(n = 1500 rpm, ap = 0.6 mm) and D(n = 3000 rpm, ap = 0.4 mm) are located in the unstable cutting region. All the dynamic responses in different points are measured and investigated, and only the dynamic response and its spectrum in points A, B, C, and D are shown in Figure 9. From Figure 9a,c, it can be seen that there is only the tool tooth passing frequency (i.e., 200 Hz, 400 Hz, 600 Hz, 650 Hz, 666 Hz, 833 Hz, 875 Hz, and 917 Hz.). From Figure 9b,d, it can be seen that the chatter frequency (i.e., 520 Hz, 580 Hz, 620 Hz, 660 Hz, 690 Hz, 790 Hz, 860 Hz, and 910 Hz.) occurs besides at the tool tooth passing frequency (i.e., 100 Hz, 200 Hz, 300 Hz, 400 Hz, 500 Hz, 600 Hz).
In addition, in order to more clearly investigate the milling chatter, the morphologies of the machined surface at different points are shown in Figure 10. From Figure 10, we can see that the machining chatter occurs at observation points B and D, which can be observed from the rough surface quality and obvious vibration.

7. Conclusions

(1) A novel updated FDM is proposed to predict the milling SLD. The cubic-spline interpolation and the Newton interpolation are introduced to approximate the state item and time-delay item, respectively. A discrete map is established between the current state matrix and the previous state matrix, and the SLD is obtained based on the eigenvalue modulus judgement criterion of the transition matrix.
(2) An iterative algorithm is proposed to obviously improve computational efficiency. The calculation of the transition matrix eigenvalues in the chattering region is eliminated. The simulation results of a benchmark example with two different radial immersion ratios show that the algorithm has a faster computational efficiency than other methods, especially when the radial immersion ratio is large.
(3) The proposed method has obvious advantages in terms of computational accuracy and convergence speed than other methods. In terms of calculation accuracy, it already coincides with the reference curve when the discrete number m is small whether the radial immersion ratio is large or small. In addition, it has a faster convergence speed both in the stable or unstable region, and this part will be further studied in the future.
(4) A series of milling experiments under different spindle speeds are designed to verify the accuracy of the proposed method. The experimental results show that the proposed method is in good agreement with the experimental value.

Author Contributions

Conceptualization, J.M.; methodology, J.M. and Y.L.; software, Y.L.; formal analysis, Y.L. and X.P.; writing—original draft preparation, J.M. and Y.L.; investigation, Y.L., X.P. and G.W.; writing—review and editing, D.Z. and B.Z.; resources, J.M. and X.P.; funding acquisition, J.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 (no. 52005166), the Postdoc-toral Research Foundation of China (no. 2019M652534), the Henan Postdoctoral Foundation (no. 19030071), the Foundation of Henan Educational Committee (no. 20A460016), Young Backbone Teachers Foundation Scheme of Henan Polytechnic University (no. 2019XQG-01), the Science and Technology Department of Henan Province (no. 202102210082), and the Doctor Foundation from Henan Polytechnic University (no.B2019-49).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Insperger, T.; Stepan, G. Updated semi-discretization method for periodic delay-differential equat with discrete delay. Int. J. Numer. Meth. Eng. 2004, 61, 117–141. [Google Scholar] [CrossRef]
  2. Insperger, T.; Gradišek, J.; Kalveram, M.; Stépán, G.; Winert, K.; Govekar, E. Machine Tool chatter and surface location error in milling processes. J. Manuf. Sci. Eng. 2006, 128, 913–920. [Google Scholar] [CrossRef]
  3. Insperger, T.; Stepan, G.; Turi, J. State-dependent delay in regenerative turning processes. Nonlinnear Dynam. 2007, 47, 275–283. [Google Scholar] [CrossRef]
  4. Chandiramani, N.K.; Pothala, T. Dynamics of 2-dof regenerative chatter during turning. J. Sound Vib. 2006, 290, 448–464. [Google Scholar] [CrossRef] [Green Version]
  5. Eksioglu, C.; Kilic, Z.M.; Altintas, Y. Discrete-Time prediction of chatter stability, cutting forces, and surface location errors in flexible milling systems. J. Manuf. Sci. Eng. 2012, 134, 061006. [Google Scholar] [CrossRef]
  6. Dun, Y.C.; Zhu, L.D.; Wang, S.H. Multi-modal method for chatter stability prediction and control in milling of thin-walled workpiece. Appl. Math. Model. 2020, 80, 602–624. [Google Scholar] [CrossRef]
  7. Li, H.; Li, S.; Shao, M.L.; Ling, J.H.; Zhao, C.Y.; Wen, B.C. Analysis of the stability of intermittent cut with energy lost. Int. J. Mater. Prod. Technol. 2014, 48, 270–285. [Google Scholar] [CrossRef]
  8. Guo, M.X.; Zhu, L.D.; Yan, B.L.; Guan, Z.H. Research on the milling stability of thin-walled parts based on the semi-discretization method of improved Runge-Kutta method. Int. J. Adv. Manuf. Technol. 2021, 115, 2325–2342. [Google Scholar] [CrossRef]
  9. Ozoegwu, P.; Eberhard, P. Stability analysis of multi-discrete delay milling with helix effects using a general order full-discretization method updated with a generalized integral quadrature. Mathematics 2020, 8, 1003. [Google Scholar] [CrossRef]
  10. Liu, X.L.; Li, R.Y.; Wu, S.; Yang, L.; Yue, C.X. A prediction method of milling chatter stability for complex surface mold. Int. J. Adv. Manuf. Technol. 2016, 89, 2637–2648. [Google Scholar] [CrossRef]
  11. Grossi, N.; Scippa, A.; Sallese, L.; Sato, R.; Campatelli, G. Spindle speed ramp-up test: A novel experimental approach for chatter stability detection. Int. J. Mach. Tools Manuf. 2015, 89, 221–230. [Google Scholar] [CrossRef]
  12. Altintas, Y.; Budak, E. Analytical prediction of stability lobes in milling. CIRP Ann.-Manuf. Technol. 1995, 44, 357–362. [Google Scholar] [CrossRef]
  13. Merdol, S.D.; Altintas, Y. Multi frequency solution of chatter stability for low immersion milling. J. Manuf. Sci. Eng. 2004, 126, 459–466. [Google Scholar] [CrossRef]
  14. Faassen, R.P.H.; Wouw, N.; Nijmeijer, H.; Oosterling, J.A.J. An improved tool path model including periodic delay for chatter prediction in milling. J. Comput. Nonliner Dyn. 2007, 2, 167–179. [Google Scholar] [CrossRef]
  15. Ding, Y.; Zhu, L.M.; Zhang, X.J.; Ding, H. Numerical integration method for prediction of milling stability. J. Manuf. Sci. Eng. 2011, 133, 031005. [Google Scholar] [CrossRef]
  16. Liang, X.G.; Yao, Z.Q.; Luo, L.; Hu, J. An improved numerical integration method for predicting milling stability with varying time delay. Int. J. Adv. Manuf. Technol. 2013, 68, 1967–1976. [Google Scholar] [CrossRef]
  17. Li, Z.Q.; Yang, Z.K.; Peng, Y.R.; Zhu, F.; Ming, X.Z. Prediction of chatter stability for milling process using Runge-Kutta-based complete discretization method. Int. J. Adv. Manuf. Technol. 2016, 86, 943–952. [Google Scholar] [CrossRef]
  18. Dai, Y.N.; Li, H.K.; Xing, X.Y.; Hao, B.T. Prediction of chatter stability for milling process using precise integration method. Precis. Eng. 2018, 52, 152–157. [Google Scholar] [CrossRef]
  19. Insperger, T.; Stepan, G. Semi-discretization method for delayed systems. Int. J. Numer. Meth. Eng. 2002, 55, 503–518. [Google Scholar] [CrossRef]
  20. Ding, Y.; Zhu, L.M.; Zhang, X.J.; Ding, H. A full-discretization method for prediction of milling stability. Int. J. Mach. Tools Manuf. 2010, 50, 502–509. [Google Scholar] [CrossRef]
  21. Ding, Y.; Zhu, L.M.; Zhang, X.J.; Ding, H. Second-order full-discretization method for milling stability prediction. Int. J. Mach. Tools Manuf. 2010, 50, 926–932. [Google Scholar] [CrossRef]
  22. Quo, Q.; Sun, Y.W.; Jiang, Y. On the accurate calculation of milling stability limits using third-order full-discretization method. Int. J. Mach. Tools Manuf. 2012, 62, 61–66. [Google Scholar] [CrossRef]
  23. Ozoegwu, C.G.; Omenyi, S.N.; Ofochebe, S.M. Hyper-third order full-discretization methods in milling stability prediction. Int. J. Mach. Tools Manuf. 2015, 92, 1–9. [Google Scholar] [CrossRef]
  24. Liu, Y.L.; Zhang, D.H.; Wu, B.H. An efficient full-discretization method for prediction of milling stability. Int. J. Mach. Tools Manuf. 2012, 63, 44–48. [Google Scholar] [CrossRef]
  25. Tang, X.W.; Peng, F.Y.; Yan, R.; Gong, Y.H.; Li, Y.T.; Jiang, L.L. Accurate and efficient prediction of milling stability with updated full-discretization method. Int. J. Adv. Manuf. Technol. 2017, 88, 2357–2368. [Google Scholar] [CrossRef]
  26. Lakshmikantham, V.; Trigiante, D. Theory of Difference Equations: Numerical Methods and Applications; Academic Press: London, UK, 1998; pp. 90–99. [Google Scholar]
  27. Insperger, T. Full-discretization and semi-discretization for milling stability prediction: Some comments. Int. J. Mach. Tools Manuf. 2010, 50, 658–662. [Google Scholar] [CrossRef]
  28. Altintas, Y. Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibration, and CNC Design, 2nd ed.; Cambridge University Press: Cambridge, UK, 2000; pp. 44–47. [Google Scholar]
Figure 1. Dynamic model of milling a workpiece.
Figure 1. Dynamic model of milling a workpiece.
Micromachines 13 00160 g001
Figure 2. The convergence rate of eigenvalues for different discrete numbers m for the 0th SDM, 1st FDM, 2nd FDM, 3rd FDM, Hermite FDM, and the proposed method. (a) Axial depth of cut ap is 0.1 mm, exact eigenvalue |μ0| is 0.7368. (b) Axial depth of cut ap is 0.2 mm, exact eigenvalue |μ0| is 0.8192. (c) Axial depth of cut ap is 0.5 mm, exact eigenvalue |μ0| is 1.0726. (d) Axial depth of cut ap is 0.8 mm, exact eigenvalue |μ0| is 1.2880.
Figure 2. The convergence rate of eigenvalues for different discrete numbers m for the 0th SDM, 1st FDM, 2nd FDM, 3rd FDM, Hermite FDM, and the proposed method. (a) Axial depth of cut ap is 0.1 mm, exact eigenvalue |μ0| is 0.7368. (b) Axial depth of cut ap is 0.2 mm, exact eigenvalue |μ0| is 0.8192. (c) Axial depth of cut ap is 0.5 mm, exact eigenvalue |μ0| is 1.0726. (d) Axial depth of cut ap is 0.8 mm, exact eigenvalue |μ0| is 1.2880.
Micromachines 13 00160 g002
Figure 3. A comparison of computational accuracy in a/D = 1.
Figure 3. A comparison of computational accuracy in a/D = 1.
Micromachines 13 00160 g003
Figure 4. A comparison of computational accuracy in a/D = 0.1.
Figure 4. A comparison of computational accuracy in a/D = 0.1.
Micromachines 13 00160 g004
Figure 5. Comparison of calculated time among the 0th SDM, 1st FDM, 2nd FDM, 3rd FDM, Hermite FDM, and the proposed method with different radial immersion ratios. (a) Calculated time at a/D = 0.1 and (b) Calculated time at a/D = 1.
Figure 5. Comparison of calculated time among the 0th SDM, 1st FDM, 2nd FDM, 3rd FDM, Hermite FDM, and the proposed method with different radial immersion ratios. (a) Calculated time at a/D = 0.1 and (b) Calculated time at a/D = 1.
Micromachines 13 00160 g005
Figure 6. Configuration of the experiments.
Figure 6. Configuration of the experiments.
Micromachines 13 00160 g006
Figure 7. (a) Distributed points 1, 2, and 3 on the thin-walled plate for impact tests. (b) Position of the sensor and point 2 for measuring responses. (c) Position of the sensor and cutter region for machining.
Figure 7. (a) Distributed points 1, 2, and 3 on the thin-walled plate for impact tests. (b) Position of the sensor and point 2 for measuring responses. (c) Position of the sensor and cutter region for machining.
Micromachines 13 00160 g007
Figure 8. Stability lobe diagrams at the measured point 2.
Figure 8. Stability lobe diagrams at the measured point 2.
Micromachines 13 00160 g008
Figure 9. Acceleration signal and its spectrum at point A, B, C, and D in the milling of the thin-walled plate.
Figure 9. Acceleration signal and its spectrum at point A, B, C, and D in the milling of the thin-walled plate.
Micromachines 13 00160 g009
Figure 10. Surface morphology of the thin plate in the milling region.
Figure 10. Surface morphology of the thin plate in the milling region.
Micromachines 13 00160 g010
Table 1. The properties of the workpiece and the specifications of the cutter.
Table 1. The properties of the workpiece and the specifications of the cutter.
CutterDiameter (mm)Number of FlutesHelix Angle (°)Length (mm)
1223075
WorkpieceDensity (g/cm3)Possion’s RationYoung’s Modulus (GPa)Material
4.60.34108Ti-6Al-4V
Table 2. Modal parameters of the cutter and workpiece.
Table 2. Modal parameters of the cutter and workpiece.
SystemNatural Frequency ω (Hz)Damping Ratio ζStiffness k (N∙m−1)
Workpiece (Mode no.1)5750.0071.19 × 106
Workpiece (Mode no.2)18200.0121.31 × 107
Cutter in X direction21260.0361.45 × 108
Cutter in Y direction21340.0331.47 × 108
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ma, J.; Li, Y.; Zhang, D.; Zhao, B.; Wang, G.; Pang, X. A Novel Updated Full-Discretization Method for Prediction of Milling Stability. Micromachines 2022, 13, 160. https://doi.org/10.3390/mi13020160

AMA Style

Ma J, Li Y, Zhang D, Zhao B, Wang G, Pang X. A Novel Updated Full-Discretization Method for Prediction of Milling Stability. Micromachines. 2022; 13(2):160. https://doi.org/10.3390/mi13020160

Chicago/Turabian Style

Ma, Junjin, Yunfei Li, Dinghua Zhang, Bo Zhao, Geng Wang, and Xiaoyan Pang. 2022. "A Novel Updated Full-Discretization Method for Prediction of Milling Stability" Micromachines 13, no. 2: 160. https://doi.org/10.3390/mi13020160

APA Style

Ma, J., Li, Y., Zhang, D., Zhao, B., Wang, G., & Pang, X. (2022). A Novel Updated Full-Discretization Method for Prediction of Milling Stability. Micromachines, 13(2), 160. https://doi.org/10.3390/mi13020160

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