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

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: 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 kh ≤ t ≤ (k + 1)h:

Figure 1

Dynamic model of milling a workpiece.

Let X(kh) = X with k = 0, 1,…, m, when t = (k + 1)h, Equation (2) can be equivalently converted to the following form:

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 X1, X, X1, X2. In addition, two other constraints and are used in cubic spline interpolation. Namely, let A stands for A(kh):

At the time interval [t, t+1], the state item X(kh + h − ε) can be approximated by cubic spline interpolation, resulting in: where

The time delay item X(kh + h − ε − T) in Equation (3) can be approximately expressed at four points X, X, X+1, X− by the third-order Newton interpolation method as follows: where

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

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

Define:

In addition, the coefficients in Equations (19)–(26) can be expressed as:

Then, if the matrix P+1 in Equation (18) is nonsingular, Equation (18) can be given by:

According to Equation (49), a discrete map could be defined as: where

In addition, Q and R can be expressed as:

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:

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 P+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 (h2). The LDE of FDM with the 1st, 2nd, 3rd, and Hermite are (h2), (h3), (h4), and (h3) [24,27], respectively. For the proposed method, the LDE is (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 m, and natural frequency w are 0.011, 0.03993 kg, and 1844 Hz, respectively. The cutter has two flutes. The cutting force coefficients are K = 6 × 108 and K = 2 × 108.

The spindle speed n is selected as 5000 rpm, and the axial depths of cut a is selected as 0.1 mm, 0.2 mm, 0.5 mm, and 0.80 mm, respectively. The exact eigenvalues 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 .

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.

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 a is 0.1 mm, exact eigenvalue |μ0| is 0.7368. (b) Axial depth of cut a is 0.2 mm, exact eigenvalue |μ0| is 0.8192. (c) Axial depth of cut a is 0.5 mm, exact eigenvalue |μ0| is 1.0726. (d) Axial depth of cut a is 0.8 mm, exact eigenvalue |μ0| is 1.2880.

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: where φ(t) is the angular position of the j-th cutter tooth, and φ and φ are the starting and exiting edge positions of the tool in contact with the workpiece, respectively. For down-milling, φ = arccos(2a/D − 1) and φ = π; for up-milling, φ = 0 and φ = 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 a = 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 a = 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.

Figure 3

A comparison of computational accuracy in a/D = 1.

Figure 4

A comparison of computational accuracy in a/D = 0.1.

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.

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.

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.

Figure 6

Configuration of the experiments.

Table 1

The properties of the workpiece and the specifications of the cutter.

Cutter	Diameter (mm)	Number of Flutes	Helix Angle (°)	Length (mm)
	12	2	30	75
Workpiece	Density (g/cm3)	Possion’s Ration	Young’s Modulus (GPa)	Material
	4.6	0.34	108	Ti-6Al-4V

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:

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 K = 1120.8 N/mm2, K = 2285.6 N/mm2, K = 9.16 N/mm, and K = 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.

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.

Table 2

Modal parameters of the cutter and workpiece.

System	Natural Frequency ω (Hz)	Damping Ratio ζ	Stiffness k (N∙m−1)
Workpiece (Mode no.1)	575	0.007	1.19 × 106
Workpiece (Mode no.2)	1820	0.012	1.31 × 107
Cutter in X direction	2126	0.036	1.45 × 108
Cutter in Y direction	2134	0.033	1.47 × 108

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, a = 0.2 mm) and C(n = 2500 rpm, a = 0.4 mm) are stable parameters, while the points B(n = 1500 rpm, a = 0.6 mm) and D(n = 3000 rpm, a = 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).

Figure 8

Stability lobe diagrams at the measured point 2.

Figure 9

Acceleration signal and its spectrum at point A, B, C, and D in the milling of the thin-walled plate.

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.

Figure 10

Surface morphology of the thin plate in the milling region.

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.