Next Article in Journal
Wavelength Modulation Characteristics of Metal Gratings on Si-Based Blocked-Impurity-Band (BIB) Terahertz Detectors
Next Article in Special Issue
Study on the Material Removal Mechanism of Ultrasonic Elliptical Vibration Cutting of Medical β Titanium Alloy
Previous Article in Journal
A 66–76 GHz Wide Dynamic Range GaAs Transceiver for Channel Emulator Application
Previous Article in Special Issue
Investigation on the Effect of Annealing Temperature on the Side Ohmic Contact Characteristics for Double Channel GaN/AlGaN Epitaxial Layer
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Predicting Milling Stability Based on Composite Cotes-Based and Simpson’s 3/8-Based Methods

1
Mechatronic Institute, Zhejiang Sci-Tech University, Hangzhou 310018, China
2
School of Mechanical Engineering, Hangzhou Dianzi University, Hangzhou 310018, China
*
Author to whom correspondence should be addressed.
Submission received: 28 April 2022 / Revised: 20 May 2022 / Accepted: 21 May 2022 / Published: 23 May 2022
(This article belongs to the Special Issue Advanced Manufacturing Technology and Systems)

Abstract

:
Avoiding chatter in milling processes is critical for obtaining machined parts with high surface quality. In this paper, we propose two methods for predicting the milling stability based on the composite Cotes and Simpson’s 3/8 formulas. First, a time-delay differential equation is established, wherein the regenerative effects are considered. Subsequently, it is discretized into a series of integral equations. Based on these integral equations, a transition matrix is determined using the composite Cotes formula. Finally, the system stability is analyzed according to the Floquet theory to obtain the milling stability lobe diagrams. The simulation results demonstrate that for the single degree of freedom (single-DOF) model, the convergence speed of the composite Cotes-based method is higher than that of the semi-discrete method and the Simpson’s equation method. In addition, the composite Cotes-based method demonstrates high computational efficiency. Moreover, to further improve the convergence speed, a second method based on the Simpson’s 3/8 formula is proposed. The simulation results show that the Simpson’s 3/8-based method has the fastest convergence speed when the radial immersion ratio is large; for the two degrees of freedom (two-DOF) model, it performs better in terms of calculation accuracy and efficiency.

1. Introduction

The three main types of vibration that occur during high-speed milling are free vibration, excited vibration, and self-excited vibration, and the regenerative chatter in self-excited vibration is the main cause of instability in the machining process [1]. Because this chatter typically leads to poor surface quality of the machined parts, aggravated tool wear, and even reduced service life of machine tools, it must be avoided to solve or overcome these problems [2]. In the analysis of chatter, dynamic milling processes that consider the delay effect are generally described using delay differential equations with respect to the time-periodic coefficients [3,4]. The chatter stability based on these differential equations is an important index for achieving high-performance machining [5].
In recent decades, experimental, analytical, and numerical methods have been studied to obtain stability lobe diagrams, which can be used to determine exact values avoiding chatter. Wu et al. [6] utilized the Lyapunov index to measure whether chatter occurred; however, they could only predict the system stability under specific parameters. Davies et al. [7] used the time-domain calculation method to predict unstable regions in the stability lobe diagrams, which was only applicable to small radial depths of cut. Altintas and Budak et al. [8] first presented a zero-order approximation method to quickly obtain the stability lobe diagrams, where the real and imaginary parts of the feature are used to determine the cutting parameters. It is worth mentioning that this method was only applicable to high radial depths of cut. Bayly et al. [9] employed the time-finite element method to calculate tool motion and utilized the weighted residual method to determine the Floquet transition matrix. Based on the delay differential equations, Insperger et al. [10] proposed the semi-discretization method (SDM) and the first-order semi-discretization method (1stSDM), which are widely used as the basis for evaluating other stability analysis methods in the time domain. The SDM and 1stSDM only discretized delay states and periodic coefficients. Based on direct integration, Ding et al. [11] presented a full-discretization method (FDM) to predict milling stability. Unlike the SDM and 1stSDM, the FDM carried out the linear interpolation on the state terms to improve computational efficiency. Further, Insperger [12] compared the FDM with the SDM and 1stSDM in terms of the convergence speed and computational efficiency. It was found that the convergence speed of the 1stSDM was faster than those of the SDM and FDM, while the computational efficiency of the FDM was higher than those of the SDM and 1stSDM. After the comparisons, some FDM-based methods were used to predict the milling stability. In addition to the above methods, the numerical integration methods are widely utilized to obtain the stability lobe diagrams. Ding et al. [13] approximated the integral terms with the Newton–Cotes and Gauss integral formulas to predict the milling stability. Lu et al. [14] approximated the solution of the delay differential equation with the direct integration technique. Qin et al. [15] approximated the state terms of the delay differential equation with the Chebyshev wavelets of the second kind. Moreover, Liu et al. [16] combined the numerical solution of the delay differential equation with the Simpson formula to predict milling stability (SEM), which was only suitable for large radial immersion. However, the aforementioned numerical methods cannot simultaneously guarantee fast convergence speed and high computational efficiency.
To this end, we present two numerical analysis methods, namely the composite Cotes-based method (CCM) and Simpson’s 3/8-based method (S38M), in this paper to improve the convergence speed and computational efficiency at the same time. The remainder of this paper is organized as follows. Section 2 presents the mathematical model of system motion with two degrees of freedom (DOF) and the state-space. The composite Cortes-based method and its simulation and analysis for the single-DOF model are presented in Section 3. Subsequently, the Simpson’s 3/8-based method and its simulation and analysis for the single-DOF and two-DOF models are presented in Section 4. Finally, the conclusions are stated in Section 5.

2. Mathematical Model

Taking down-milling as an example, the dynamic system of end-milling with two degrees of freedom is shown in Figure 1, where the workpiece is assumed to be rigid, and the milling cutter is assumed to be evenly distributed. Note that the helix angle of the milling cutter is not taken into account. Considering the regenerative effect, the mathematical model of system motion can be expressed as a second-order differential equation [17] as follows:
M q ¨ ( t ) + C q ˙ ( t ) + K q ( t ) = a p K c ( t ) [ q ( t ) q ( t T ) ]
where q ( t ) , q ˙ ( t ) , and q ¨ ( t ) represent the displacement vector, the first derivative of q(t) w.r.t t, and the second derivative of q(t) w.r.t t, respectively. M, C, K, and ap represent the mass matrix, damping matrix, stiffness matrix, and depth of cut, respectively. T = 60/(Nv) represents the delay period of the delay differential equations, where N represents the cutter tooth number, and v represents the spindle speed. The radial cutting force coefficient matrix Kc(t) can be expressed as follows:
K c ( t ) = [ h x x ( t ) h x y ( t ) h y x ( t ) h y y ( t ) ]
where
h x x ( t ) = j = 1 N g ( ϕ j ( t ) ) sin ( ϕ j ( t ) ) [ K t cos ( ϕ j ( t ) ) + K n sin ( ϕ j ( t ) ) ]
h x y ( t ) = j = 1 N g ( ϕ j ( t ) ) cos ( ϕ j ( t ) ) [ K t cos ( ϕ j ( t ) ) + K n sin ( ϕ j ( t ) ) ]
h y x ( t ) = j = 1 N g ( ϕ j ( t ) ) sin ( ϕ j ( t ) ) [ K t sin ( ϕ j ( t ) ) + K n cos ( ϕ j ( t ) ) ]
and
h y y ( t ) = j = 1 N g ( ϕ j ( t ) ) cos ( ϕ j ( t ) ) [ K t sin ( ϕ j ( t ) ) + K n cos ( ϕ j ( t ) ) ]
where Kt and Kn represent the tangential and normal cutting force coefficients, respectively. ϕ j ( t ) represents the angular position of tooth j, as shown in Figure 1, and is defined as follows:
ϕ j ( t ) = 2 π v 60 t + 2 π ( j 1 ) N
g ( ϕ j ( t ) ) is a piecewise function used to judge whether the tooth j is cutting or not and is defined as follows:
g ( ϕ j ( t ) ) = { 1 ϕ s t < ϕ j ( t ) < ϕ e x 0 otherwise
where ϕ s t and ϕ e x represent the start and exit cutting angles of the tooth j, respectively. In the up-milling process, ϕ s t = 0 and ϕ e x = arccos ( 1 2 a / D ) ; in the down-milling process, ϕ s t = arccos ( 2 a / D 1 ) and ϕ e x = π , where a/D represents the radial immersion ratio.

3. Composite Cotes-Based Method

3.1. State-Space Expression

x(t) can be defined as x ( t ) = [ q ( t ) M q ˙ ( t ) + C q ( t ) / 2 ] . Then, the state expression of (1) can be expressed as follows:
x ˙ ( t ) = A x ( t ) + a p B ( t ) [ x ( t ) x ( t T ) ]
where
A = [ M 1 C / 2 M 1 C M 1 C / 4 K C M 1 / 2 ]
B ( t ) = [ 0 0 K c ( t ) 0 ]
If a p B ( t ) [ x ( t ) x ( t T ) ] is regarded as the inhomogeneous term of x ˙ ( t ) = A x ( t ) , then the following equation can be derived:
x ( t ) = e A ( t t 0 ) x ( t 0 ) + a p t 0 t { e A ( t ξ ) B ( ξ ) [ x ( ξ ) x ( ξ T ) ] } d ξ
where t0 represents the start time. It is worth noting that x(ξ) is unknown such that (12) is a Volterra integral equation of the second kind.

3.2. Numerical Algorithm

To solve (9) using the composite Cotes formula, the integral equation f(x) is assumed to be continuous in the interval [c, d]. Then, [c, d] is divided into m isometric intervals; that is, the span of each interval h is equal to (d − c)/m. Based on the isometric nodes uk = c + (k − 1)/h, where k = 1, 2, …, m + 1, we obtain the following:
I m = ( d c ) k = 0 m C k ( m ) f ( u k )
where C k ( m ) represents the Cotes coefficients [18], expressed as follows:
C k ( m ) = ( 1 ) m k n k ! ( n k ) 0 m j = 0 , j k m ( t j ) d t , k = 0 , 1 , , m
It is worth noting that C k ( m ) only depends on m. When m = 4, (14) can be expressed as follows:
C = 1 90 [ 7 f ( u 0 ) + 32 f ( u 1 ) + 12 f ( u 2 ) + 32 f ( u 3 ) + 7 f ( u 4 ) ]
where T consists of the durations of the free and excited vibrations [15]. The durations of the free vibration and excited vibration can be defined as tf and tc, respectively. For the free vibration, B(ξ) = 0, and (12) can be expressed as follows:
x ( t ) = e A ( t t 0 ) x ( t 0 )
Therefore, the state equation can be solved when t = t0 + tf, and it can be expressed as follows:
x ( t 0 + t f ) = e A t f x ( t 0 )
For the excited vibration, its time lies in [t0 + tf, t0 + T]. The interval is divided into n isometric intervals, and the span of each interval l is equal to (T − tf)/n. The discrete time can be expressed as ti = t0 + tf + (i − 1)h, where i = 1, 2, …, n + 1. According to the numerical integral solution of the Volterra integral equation of the second kind [19], we obtain the following:
x ( t i + 1 ) = e A ( t i + 1 t i ) x ( t i ) + a p t i t i + 1 { e A ( t i + 1 ξ ) B ( ξ ) [ x ( ξ ) x ( ξ t ) ] } d ξ
If [ti, ti + 1] is divided into four isometric intervals, then ti, ti+1/4, ti+1/2, ti+3/4, and ti+1 can be derived. Combining (13) and (15), the integral equation can be expressed as follows:
t i t i + 1 f ( ξ ) d ξ = t i + 1 t i 90 [ 7 f ( t i ) + 32 f ( t i + 1 / 4 ) + 12 f ( t i + 1 / 2 ) + 32 f ( t i + 3 / 4 ) + 7 f ( t i + 1 ) ]
Upon substituting (19) into (18), we can see that x(ti+1/4), x(ti+1/2), and x(ti+3/4) cannot be solved directly. The barycentric Lagrange interpolation method [20] is introduced here to approximate x(ti+1/4), x(ti+1/2), and x(ti+3/4) and obtain the following:
x ( t i + 1 / 4 ) 21 x ( t i ) + 14 x ( t i + 1 ) 3 x ( t i + 2 ) 32
x ( t i + 1 / 2 ) 3 x ( t i ) + 6 x ( t i + 1 ) x ( t i + 2 ) 8
x ( t i + 3 / 4 ) 5 x ( t i ) + 30 x ( t i + 1 ) 3 x ( t i + 2 ) 32
Substituting (19)–(22) into (18), x(ti) can be expressed as follows:
x ( t i + 1 ) = e A h x ( t i ) + a p h 6 { 5 2 e A h B ( t i ) ( x ( t i ) x ( t i T ) ) + 4 B ( t i + 1 ) ( x ( t i + 1 ) x ( t i + 1 T ) ) 1 2 e A h B ( t i + 2 ) ( x ( t i + 2 ) x ( t i + 2 T ) ) }
It is worth noting that B(tn+2) cannot be calculated directly using (9) when i = n. A Newton Cotes formula [13] is introduced here to calculate x(tn+1) and is expressed as follows:
x ( t n + 1 ) = e A h x ( t n ) + a p h 2 [ e A h B ( t n ) ( x ( t n ) x ( t n T ) ) + B ( t n + 1 ) ( x ( t n + 1 ) x ( t n + 1 T ) ) ]
Combining (17), (23), and (24), the transition matrix can be constructed as follows:
C 1 = [ 0 0 0 0 0 0 e A h 0 0 0 0 0 0 e A h 0 0 0 0 0 0 0 0 0 0 0 0 0 e A h 0 0 0 0 0 0 e A h 0 ] ( 2 n + 2 ) × ( 2 n + 2 )
D 1 = [ 0 0 0 0 0 0 75 2 e A h B 1 60 B 2 15 2 e A h B 3 0 0 0 0 75 2 e A h B 2 60 B 3 0 0 0 0 0 0 60 B i 1 15 2 e A h B i 0 0 0 0 75 2 e A h B i 1 60 B i 15 2 e A h B i + 1 0 0 0 0 45 e A h B i 45 B i + 1 ] ( 2 n + 2 ) × ( 2 n + 2 )
E 1 = [ 0 0 0 e A t f 0 0 0 0 0 0 0 0 0 0 0 0 ] ( 2 n + 2 ) × ( 2 n + 2 )
where Bi represents B(ti), i = 1, 2, …, n + 1.
According to (25)–(27), the dynamic mapping of the discrete system can be determined as follows:
( I 1 + C 1 + a p h 90 D 1 ) [ x ( t 1 ) x ( t 2 ) x ( t n + 1 ) ] = ( a p h 90 D 1 + E 1 ) [ x ( t 1 T ) x ( t 2 T ) x ( t n + 1 T ) ]
where I1 is a (2n + 1) × (2n + 1) identity matrix.
Therefore, the transition matrix of the dynamic system in a tool pass cycle is obtained as follows:
Φ 1 = ( I 1 + C 1 + a p h 90 D 1 ) 1 ( a p h 90 D 1 + E 1 )
It is worth noting that the chatter stability must be determined using the Floquet theory. If the modulus of any eigenvalue in Φ1 exceeds 1, the system is unstable; if the moduli of all eigenvalues in Φ1 are less than 1, the system is stable [21,22]. Therefore, the boundary curve between the unstable and stable regions in the lobe diagram can be used as the criterion for judging whether chatter occurs.

3.3. Simulation and Analysis

In this section, SDM [10] and SEM [16] are compared. For objective comparison, the CCM, SDM, and SEM share the same parameters and machining conditions. A benchmark example of the single-DOF milling model is utilized to validate and analyze the CCM. The single-DOF milling mathematical model is expressed as follows [23]:
x ¨ ( t ) + 2 ζ ω n x ˙ ( t ) + ω n 2 x ( t ) = a p h ( t ) m t ( x ( t ) x ( t T ) )
where ζ denotes the relative damping, ωn denotes the angular natural frequency, and h(t) denotes the specific cutting force coefficient.
The state-space for single-DOF milling mathematical model is expressed as follows:
x ˙ ( t ) = A x ( t ) + a p B ( t ) [ x ( t ) x ( t T ) ]
where
A = [ ζ ω n 1 m t m t ω n 2 ( ζ 2 1 ) ζ ω n ]
B ( t ) = [ 0 0 h ( t ) 0 ]
where h(t) is equal to hxx(t), defined in (3).
The parameters in (32) and (33) are defined as follows: fn = 922 Hz, ξ = 0.011, mt = 0.3993 kg, Kt = 6 × 108 N/m2 and Kt = 2 × 108 N/m2. The adopted machining condition is down-milling. The radial immersion ratio a/D is set as 1 to avoid intermittent milling. The spindle speed is set as 5000 rpm (v = 5000 rpm). The depths of cut are set as 0.001 m, 0.0005 mm, and 0.0002 mm, respectively. All programs in this study are executed in MATLAB R2019a and run on a personal computer (AMD Ryzen 5 5600H; CPU 4.0 GHz, 16 GB). The maximum modulus of the eigenvalues of Φ is labelled as | λ | , and the maximum modulus of the eigenvalues of Φ with SDM is labelled as | λ 0 | . In this study, | λ 0 | is treated as the exact value.
The convergence rate comparisons of the CCM, SDM, and SEM are shown in Figure 2. It is seen from the figure that when n is small, the convergence rate of the CCM is significantly higher than those of the SDM and SEM regardless of the depth of cut. As n increases, the convergence rates of the three methods approach gradually, and the convergence rates of the CCM and SEM are higher than that of SDM. The figure (a) is enlarged to observe the differences better. It is worth noting that the convergence rate of the CCM is significantly higher than those of the SDM and SEM when the discrete number changes from 75 to 100. The results demonstrate that the convergence rate of the CCM outweighs those of the SDM and SEM.
The stability lobe diagrams can be used to compare the calculation accuracy among the CCM, SDM, and SEM. As a reference, the discrete number is set as 500. The other parameters are set as follows: the discrete numbers are set as 25, 40, and 55; the radial immersion ratios are set as 0.05, 0.5, and 1; the spindle speed varies from 0.5 × 104 to 2.5 × 104 rpm, and 200 equally distributed sampling points are selected within this range; the depth of cut varies from 0 to 0.01 m, and 100 equally distributed sampling points are selected within this range.
In the single-DOF system, the stability lobe diagrams obtained with the CCM, SDM, and SEM when a/D = 1 are given in Table 1. It can be seen that the stability lobe diagrams of the CCM are closer to the reference ones than the SDM and SEM. Note that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25. Overall, the CCM is more accurate than the SDM and SEM. To further compare the calculation accuracy, two indicators are introduced, i.e., the arithmetic mean of relative error (AMRE) and mean squared error (MSE) [2]. The AMRE and MSE are defined as follows [2]:
A M R E = 1 n i = 1 n | a i a i 0 | a i 0
M S E = 1 n i = 1 n ( a i a i 0 ) 2
where ai and ai0 denote the predicted axial depth and reference axial depth, respectively. n denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the AMRE and MSE are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the AMRE and MSE of the SEM cannot be attained. Take n = 40 as an example. The AMRE obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The MSE obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07 × 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM.
Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when a/D = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06 s to 30.13 s, and 81.21 s to 35.18 s when the discrete numbers are 40, 45, 50, and 55, respectively. Furthermore, the stability lobe diagrams, AMRE and MSE, as well as the computational time when a/D = 0.05 and a/D = 0.5 are given in Table 2 and Table 3 and Figure 5 and Figure 6. It could be noted that the stability lobe diagrams obtained with the CCM are closer to the reference ones than those obtained with the SDM and SEM, and the CCM has the highest calculation accuracy. The results also demonstrate that the CCM has significant advantages in the computational efficiency.

4. Simpson’s 3/8-Based Method

4.1. State-Space Expression

To match with S38M, a new state-space expression is defined as follows [17]:
x ˙ ( t ) = A x ( t ) + a p B ( t ) [ x ( t ) x ( t T ) ]
where
A = [ 0 I M 1 K M 1 C ]
and
B ( t ) = [ 0 0 M 1 K c ( t ) 0 ]
Based on the numerical integration method of the Volterra integral equation of the second kind [19], the state-space equation can be deduced as follows:
x ( t ) = e A ( t t s t ) x ( t s t ) + a p t s t t { e A ( t ξ ) B ( ξ ) [ x ( ξ ) x ( ξ T ) ] } d ξ
where tst represents the start time and defaults to t1.
When t = t3, we obtain the following:
x ( t 3 ) = e A ( t 3 t 1 ) x ( t 1 ) + a p t s t t 1 { e A ( t 1 ξ ) B ( ξ ) [ x ( ξ ) x ( ξ T ) ] } d ξ

4.2. Numerical Algorithm

The Lagrange polynomial denoted by PL(x) [16,24] can be used to approximate f(x) as follows:
P L ( x ) = k = 0 L f ( x k ) η L , k ( x )
where
η L , k ( x ) = 0 m k m L x x m x L x m
It is worth noting that xi, xi+1, and xi+2 only exist when L = 2, and then, the Simpson’s 3/8 formula (41) can be derived using P2(x).
x i x i + 2 f ( x ) d x = h 3 ( f ( x i ) + 4 f ( x i + 1 ) + f ( x i + 2 ) )
Based on (41), (38) can be rewritten as follows:
x ( t 3 ) = e 2 A h x ( t 1 ) + a p h 3 { e 2 A h B ( t 1 ) [ x ( t 1 ) x ( t 1 T ) ] + 4 e A h B ( t 2 ) [ x ( t 2 ) x ( t 2 T ) ] + B ( t 3 ) [ x ( t 3 ) x ( t 3 T ) ] }
which can be expressed as follows:
x ( t i + 2 ) = e 2 A h x ( t i ) + a p h 3 { e 2 A h B ( t i ) [ x ( t i ) x ( t i T ) ] + 4 e A h B ( t i + 1 ) [ x ( t i + 1 ) x ( t i + 1 T ) ] + B ( t i + 2 ) [ x ( t i + 2 ) x ( t i + 2 T ) ] }
As B(ti+2) cannot be solved directly when i = n, the classical numerical integration method is introduced to determine x(tn+1) as follows:
x ( t n + 1 ) = e A h x ( t n ) + a p h 2 { e A h B ( t n ) [ x ( t n ) x ( t n T ) ] + B ( t n + 1 ) [ x ( t n + 1 ) x ( t n + 1 T ) ] }
Combining (17), (43), and (44), the transition matrix can be constructed as follows:
C 2 = [ 0 0 0 0 0 0 e 2 A h 0 0 0 0 0 0 e 2 A h 0 0 0 0 0 0 0 0 0 0 0 0 0 e 2 A h 0 0 0 0 0 0 e A h 0 ] ( 2 n + 2 ) × ( 2 n + 2 )
D 2 = [ 0 0 0 0 0 0 e 2 A h B 1 4 e A h B 2 B 3 0 0 0 0 e 2 A h B 2 4 e A h B 3 0 0 0 0 0 0 4 e A h B i 1 B i 0 0 0 0 e 2 A h B i 1 4 e A h B i B i + 1 0 0 0 0 3 2 e A h B i 3 2 B i + 1 ] ( 2 n + 2 ) × ( 2 n + 2 )
E 2 = [ 0 0 0 e A t f 0 0 0 0 0 0 0 0 0 0 0 0 ] ( 2 n + 2 ) × ( 2 n + 2 )
I 2 = [ 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 ] ( 2 n + 2 ) × ( 2 n + 2 )
According to (45)–(48), the dynamic mapping of the discrete system can be determined as follows:
( I 2 + C 2 + a p h 3 D 2 ) [ x ( t 1 ) x ( t 2 ) x ( t n + 1 ) ] = ( a p h 3 D 2 + E 2 ) [ x ( t 1 T ) x ( t 2 T ) x ( t n + 1 T ) ]
Therefore, the transition matrix of the dynamic system in a tool pass cycle is obtained as follows:
Φ 2 = ( I 2 + C 2 + a p h 3 D 2 ) 1 ( a p h 3 D 2 + E 2 )

4.3. Simulation and Analysis

4.3.1. Single-DOF Milling Model

A benchmark example of the single-DOF milling model is utilized to validate and analyze S38M. Its state-space expression is introduced as (31), where
A = [ 0 1 ω n 2 2 ζ ω n ]
B ( t ) = [ 0 0 h ( t ) m t 0 ]
To analyze the convergence rate of S38M, the radial immersion ratio a/D is set as 1 to avoid intermittent milling. The depths of cut are set as 0.0003 m, 0.0008 m, and 0.0015 m. Figure 7 shows the convergence rate comparisons of the S38M, CCM, SDM, and SEM when the spindle speed is 8000 rpm. It is worth noting that the S38M exhibits the highest convergence rate compared with other chatter analysis methods. According to the enlarged figures on the right, the S38M has a good convergence rate even when the discrete number is large. Additionally, the results show that the convergence rate of the SEM fluctuates evidently, and the SEM is not as stable as the other methods. To further analyze the convergence rate of the S38M, the spindle speed is raised to 10,000 rpm, and the other parameters remain unchanged. Figure 8 shows the resultant convergence rates. It can be seen that these results are consistent with the previous results.
The low-immersion condition can be utilized to verify the stability of the convergence rate [25]. The radial immersion ratio is set as 0.05; that is, a/D = 0.05. The depths of cut are set as 0.0001 m, 0.0002 m, and 0.0003 m. The spindle speed was set to 5000 rpm. Figure 9 shows the convergence rate comparisons of S38M, CCM, SDM, and SEM under low immersion. The convergence rates arranged from high to low are in the following order: S38M, CCM, SDM, and SEM. Moreover, it is seen that the convergence rate of SEM fluctuated remarkably.
The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.
As described above, the fluctuations appear in Figure 9. It is found that for the S38M, the convergence rate is often a local maximum when the discrete number is an odd one; the convergence rate is often a local minimum when the discrete number is an even one. To further study this problem, the discrete numbers are set as 30, 40, 50, 60, 70, and 80. The corresponding AMRE and MSE are plotted in Figure 12. The results show that the S38M attains a better computational accuracy when the discrete number is an even one.

4.3.2. Two-DOF Milling Model

The two-DOF milling mathematical model is expressed as follows [17]:
[ m t 0 0 m t ] [ x ¨ ( t ) y ¨ ( t ) ] + [ 2 m t ζ ω n 0 0 2 m t ζ ω n ] [ x ˙ ( t ) y ˙ ( t ) ] + [ m t ω n 2 0 0 m t ω n 2 ] [ x ( t ) y ( t ) ] = a p [ h x x ( t ) h x y ( t ) h y x ( t ) h y y ( t ) ] × { [ x ( t ) y ( t ) ] [ x ( t T ) y ( t T ) ] }
where the tool is symmetrical by default. mt, ζ, and ωn in the x-direction are identical to those in the y-direction. hxx(t), hxy(t), hyx(t), and hyy(t) are the same as those defined in (3)–(6).
A new state vector is defined as r ˙ ( t ) = [ x ( t ) , y ( t ) , x ˙ ( t ) , y ˙ ( t ) ] T . Then, (53) can be expressed as follows:
r ˙ ( t ) = A r ( t ) + a p B ( t ) [ r ( t ) r ( t T ) ]
where
A = [ 0 0 1 0 0 0 0 1 ω n 2 0 2 ζ ω n 0 0 ω n 2 0 2 ζ ω n ]
and
B = [ 0 0 0 0 0 0 0 0 h x x ( t ) / m t h x y ( t ) / m t 0 0 h y x ( t ) / m t h y y ( t ) / m t 0 0 ]
The following deduction is the same as that in Section 4.1.
The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.

5. Conclusions

In this study, we focused on the prediction of milling stability using composite Cotes-based and Simpson’s 3/8-based methods. First, the composite Cotes-based method is proposed for preventing chatter in milling processes. The three steps for obtaining the milling stability lobe diagrams are as follows: (1) establish the time-delay differential equation while considering the regenerative effects; (2) determine the transition matrix using the integral equations; (3) analyze the system stability according to the Floquet theory. To measure the calculation accuracy, the AMRE and MSE are adopted in our work. The results demonstrate that for the single-DOF model, when the discrete number is 40, the AMRE and MSE can be respectively reduced from 0.076 to 0.041 and from 4.24 × 10−8 to 2.62 × 10−8, and the calculation time can be reduced from 55.63s to 20.21s by using the proposed composite Cotes-based method. In addition, the Simpson’s 3/8-based method is proposed to further improve the calculation efficiency. The results demonstrate that for the single-DOF model, the proposed Simpson’s 3/8-based method significantly improves the convergence speed while sharing the same computation time with the composite Cotes-based method when the radial immersion ratio is large; for the two-DOF model, the accuracy of the stability lobe diagrams of the Simpson’s 3/8-based method are better than those of the SDM, and the computation time is remarkably saved using the Simpson’s 3/8-based method.

Author Contributions

Conceptualization, X.D. and J.Z.; methodology, X.D.; software, P.R.; validation, X.D., P.R. and J.Z.; formal analysis, P.R.; investigation, X.D.; resources, X.D.; data curation, P.R.; writing—original draft preparation, X.D.; writing—review and editing, J.Z.; visualization, X.D.; supervision, J.Z.; project administration, J.Z.; funding acquisition, J.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by [National Natural Science Foundation of China] grant number [52005142] and [Open Fund of State Key Laboratory of Robotics and System] grant number [SKLRS-2021-KF-08]. The authors would like to thank the Mechatronic Institute, Zhejiang Sci-Tech University.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. 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. [Google Scholar] [CrossRef] [PubMed]
  2. Yan, Z.; Zhang, C.; Jia, J.; Ma, B.; Jiang, X.; Wang, D.; Zhu, T. High-order semi-discretization methods for stability analysis in milling based on precise integration. Precis. Eng. 2021, 73, 71–92. [Google Scholar] [CrossRef]
  3. Catania, G.; Mancinelli, N. Theoretical–experimental modeling of milling machines for the prediction of chatter vibration. Int. J. Mach. Tools Manuf. 2011, 51, 339–348. [Google Scholar] [CrossRef]
  4. Fu, Z.; Zhang, X.; Wang, X.; Yang, W. Analytical modeling of chatter vibration in orthogonal cutting using a predictive force model. Int. J. Mech. Sci. 2014, 88, 145–153. [Google Scholar] [CrossRef]
  5. Gradišek, J.; Kalveram, M.; Insperger, T.; Weinert, K.; Stepan, G.; Govekar, E.; Grabec, I. On stability prediction for milling. Int. J. Mach. Tools Manuf. 2005, 45, 769–781. [Google Scholar] [CrossRef]
  6. Wu, S.; Li, R.; Liu, X.; Yang, L.; Zhu, M. Experimental Study of Thin Wall Milling Chatter Stability Nonlinear Criterion. Procedia CIRP 2016, 56, 422–427. [Google Scholar] [CrossRef] [Green Version]
  7. Davies, M.; Pratt, J.; Dutterer, B.; Burns, T. The Stability of Low Radial Immersion Milling. CIRP Ann. 2000, 49, 37–40. [Google Scholar] [CrossRef]
  8. Altintaş, Y.; Budak, E. Analytical Prediction of Stability Lobes in Milling. CIRP Ann. 1995, 44, 357–362. [Google Scholar] [CrossRef]
  9. Bayly, P.V.; Mann, B.P.; Schmitz, T.L.; Peters, D.A.; Stepan, G.; Insperger, T. Effects of Radial Immersion and Cutting Direction on Chatter Instability in End-Milling. Manufacturing 2002, 3641, 351–363. [Google Scholar] [CrossRef]
  10. Insperger, T.; Stépán, G. Updated semi-discretization method for periodic delay-differential equations with discrete delay. Int. J. Numer. Methods Eng. 2010, 61, 117–141. [Google Scholar] [CrossRef]
  11. 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]
  12. 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]
  13. 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]
  14. Lu, Y.; Ding, Y.; Peng, Z.; Chen, Z.C.; Zhu, L. A spline-based method for stability analysis of milling processes. Int. J. Adv. Manuf. Technol. 2017, 89, 2571–2586. [Google Scholar] [CrossRef]
  15. Qin, C.; Tao, J.; Shi, H.; Xiao, D.; Li, B.; Liu, C. A novel Chebyshev-wavelet-based approach for accurate and fast prediction of milling stability. Precis. Eng. 2020, 62, 244–255. [Google Scholar] [CrossRef]
  16. Liu, C.; Tang, D.; Li, S.; Ding, G. Simpson’s 3/8–based method stability analysis for milling processes. Int. J. Adv. Manuf. Technol. 2021, 114, 671–682. [Google Scholar] [CrossRef]
  17. Niu, J.; Ding, Y.; Zhu, L.; Ding, H. Runge–Kutta methods for a semi-analytical prediction of milling stability. Nonlinear Dyn. 2013, 76, 289–304. [Google Scholar] [CrossRef]
  18. Considine, G.D.; Kulik, P.H. (Eds.) Newton-Cotes formula. In Van Nostrand’s Scientific Encyclopedia; John Wiley & Sons: Chichester, UK, 2005. [Google Scholar] [CrossRef]
  19. Odibat, Z.M. Differential transform method for solving Volterra integral equation with separable kernels. Math. Comput. Model. 2008, 48, 1144–1149. [Google Scholar] [CrossRef]
  20. Berrut, J.P.; Trefethen, L.N. Barycentric lagrange lnterpolation. Siam Rev. 2004, 46, 501. [Google Scholar] [CrossRef] [Green Version]
  21. Hajdu, D.; Borgioli, F.; Michiels, W.; Insperger, T.; Stepan, G. Robust stability of milling operations based on pseudospectral approach. Int. J. Mach. Tools Manuf. 2019, 149, 103516. [Google Scholar] [CrossRef]
  22. Qin, C.; Tao, J.; Liu, C. A predictor-corrector-based holistic-discretization method for accurate and efficient milling stability analysis. Int. J. Adv. Manuf. Technol. 2018, 96, 2043–2054. [Google Scholar] [CrossRef]
  23. Niu, J.; Ding, Y.; Zhu, L.; Ding, H. Stability Analysis of Milling Processes With Periodic Spindle Speed Variation Via the Variable-Step Numerical Integration Method. J. Manuf. Sci. Eng. 2016, 138, 114501. [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. 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]
  26. Tang, X.; Peng, F.; Yan, R.; Gong, Y.; Li, Y.; Jiang, L. Accurate and efficient prediction of milling stability with updated full-discretization method. Int. J. Adv. Manuf. Technol. 2016, 88, 2357–2368. [Google Scholar] [CrossRef]
Figure 1. Dynamic system of end-milling with two degrees of freedom.
Figure 1. Dynamic system of end-milling with two degrees of freedom.
Micromachines 13 00810 g001
Figure 2. Convergence rate comparisons of the CCM, SDM and SEM when a/D = 1. (a) ap = 0.0002 m and n [30, 100]; (b) ap = 0.0002 m and n [75, 100]; (c) ap = 0.0005 m and n [30, 100]; (d) ap = 0.001 m and n [30, 100].
Figure 2. Convergence rate comparisons of the CCM, SDM and SEM when a/D = 1. (a) ap = 0.0002 m and n [30, 100]; (b) ap = 0.0002 m and n [75, 100]; (c) ap = 0.0005 m and n [30, 100]; (d) ap = 0.001 m and n [30, 100].
Micromachines 13 00810 g002
Figure 3. AMRE and MSE of the CCM, SDM, and SEM when a/D = 1. (a) AMRE; (b) MSE.
Figure 3. AMRE and MSE of the CCM, SDM, and SEM when a/D = 1. (a) AMRE; (b) MSE.
Micromachines 13 00810 g003
Figure 4. Computational time of the CCM, SDM, and SEM when a/D = 1.
Figure 4. Computational time of the CCM, SDM, and SEM when a/D = 1.
Micromachines 13 00810 g004
Figure 5. AMRE and MSE of the CCM, SDM, and SEM. (a) AMRE when a/D = 0.05; (b) AMRE when a/D = 0.5; (c) MSE when a/D = 0.05; (d) MSE when a/D = 0.5.
Figure 5. AMRE and MSE of the CCM, SDM, and SEM. (a) AMRE when a/D = 0.05; (b) AMRE when a/D = 0.5; (c) MSE when a/D = 0.05; (d) MSE when a/D = 0.5.
Micromachines 13 00810 g005
Figure 6. Computational time of the CCM, SDM, and SEM. (a) a/D = 0.05; (b) a/D = 0.5.
Figure 6. Computational time of the CCM, SDM, and SEM. (a) a/D = 0.05; (b) a/D = 0.5.
Micromachines 13 00810 g006
Figure 7. Convergence rate comparisons of the S38M, CCM, SDM, and SEM when v = 8000 rpm and a/D = 1. (a) ap = 0.0003 m and n [30, 100]; (b) ap = 0.0003 m and n [75, 100]; (c) ap = 0.0008 m and n [30, 100]; (d) ap = 0.0008 m and n [75, 100]; (e) ap = 0.0015 m and n [30, 100]; (f) ap = 0.0015 m and n [75, 100].
Figure 7. Convergence rate comparisons of the S38M, CCM, SDM, and SEM when v = 8000 rpm and a/D = 1. (a) ap = 0.0003 m and n [30, 100]; (b) ap = 0.0003 m and n [75, 100]; (c) ap = 0.0008 m and n [30, 100]; (d) ap = 0.0008 m and n [75, 100]; (e) ap = 0.0015 m and n [30, 100]; (f) ap = 0.0015 m and n [75, 100].
Micromachines 13 00810 g007
Figure 8. Convergence rate comparisons of the S38M, CCM, SDM, and SEM when v = 10,000 rpm and a/D = 1. (a) ap = 0.0003 m and n [30, 100]; (b) ap = 0.0003 m and n [75, 100]; (c) ap = 0.0008 m and n [30, 100]; (d) ap = 0.0008 m and n [75, 100]; (e) ap = 0.0015 m and n [30, 100]; (f) ap = 0.0015 m and n [75, 100].
Figure 8. Convergence rate comparisons of the S38M, CCM, SDM, and SEM when v = 10,000 rpm and a/D = 1. (a) ap = 0.0003 m and n [30, 100]; (b) ap = 0.0003 m and n [75, 100]; (c) ap = 0.0008 m and n [30, 100]; (d) ap = 0.0008 m and n [75, 100]; (e) ap = 0.0015 m and n [30, 100]; (f) ap = 0.0015 m and n [75, 100].
Micromachines 13 00810 g008
Figure 9. Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and a/D = 0.05. (a) ap = 0.0001 m; (b) ap = 0.0002 m; (c) ap = 0.0003 m.
Figure 9. Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and a/D = 0.05. (a) ap = 0.0001 m; (b) ap = 0.0002 m; (c) ap = 0.0003 m.
Micromachines 13 00810 g009
Figure 10. AMRE and MSE of the S38M, CCM, SDM, and SEM. (a) AMRE when a/D = 0.05; (b) AMRE when a/D = 1; (c) MSE when a/D = 0.05; (d) MSE when a/D = 1.
Figure 10. AMRE and MSE of the S38M, CCM, SDM, and SEM. (a) AMRE when a/D = 0.05; (b) AMRE when a/D = 1; (c) MSE when a/D = 0.05; (d) MSE when a/D = 1.
Micromachines 13 00810 g010
Figure 11. Computational time of the S38M, CCM, SDM, and SEM. (a) a/D = 0.05; (b) a/D = 1.
Figure 11. Computational time of the S38M, CCM, SDM, and SEM. (a) a/D = 0.05; (b) a/D = 1.
Micromachines 13 00810 g011
Figure 12. AMRE and MSE of the S38M for the fluctuation analysis when a/D = 0.05. (a) AMRE; (b) MSE.
Figure 12. AMRE and MSE of the S38M for the fluctuation analysis when a/D = 0.05. (a) AMRE; (b) MSE.
Micromachines 13 00810 g012
Table 1. Stability lobe diagrams of the CCM, SDM, and SEM when a/D = 1.
Table 1. Stability lobe diagrams of the CCM, SDM, and SEM when a/D = 1.
n254055
CCM Micromachines 13 00810 i001 Micromachines 13 00810 i002 Micromachines 13 00810 i003
SDM Micromachines 13 00810 i004 Micromachines 13 00810 i005 Micromachines 13 00810 i006
SEM Micromachines 13 00810 i007 Micromachines 13 00810 i008 Micromachines 13 00810 i009
Table 2. Stability lobe diagrams of the CCM, SDM, and SEM when a/D = 0.05.
Table 2. Stability lobe diagrams of the CCM, SDM, and SEM when a/D = 0.05.
n254055
CCM Micromachines 13 00810 i010 Micromachines 13 00810 i011 Micromachines 13 00810 i012
SDM Micromachines 13 00810 i013 Micromachines 13 00810 i014 Micromachines 13 00810 i015
SEM Micromachines 13 00810 i016 Micromachines 13 00810 i017 Micromachines 13 00810 i018
Table 3. Stability lobe diagrams of the CCM, SDM, and SEM when a/D = 0.5.
Table 3. Stability lobe diagrams of the CCM, SDM, and SEM when a/D = 0.5.
n254055
CCM Micromachines 13 00810 i019 Micromachines 13 00810 i020 Micromachines 13 00810 i021
SDM Micromachines 13 00810 i022 Micromachines 13 00810 i023 Micromachines 13 00810 i024
SEM Micromachines 13 00810 i025 Micromachines 13 00810 i026 Micromachines 13 00810 i027
Table 4. Stability lobe diagrams of the S38M.
Table 4. Stability lobe diagrams of the S38M.
n254055
a/D = 0.05 Micromachines 13 00810 i028 Micromachines 13 00810 i029 Micromachines 13 00810 i030
a/D = 0.5 Micromachines 13 00810 i031 Micromachines 13 00810 i032 Micromachines 13 00810 i033
a/D = 1 Micromachines 13 00810 i034 Micromachines 13 00810 i035 Micromachines 13 00810 i036
Table 5. Comparisons of S38M and SDM for the two-DOF milling model.
Table 5. Comparisons of S38M and SDM for the two-DOF milling model.
S38MSDM
a/D = 0.05 Micromachines 13 00810 i037 Micromachines 13 00810 i038
a/D = 0.5 Micromachines 13 00810 i039 Micromachines 13 00810 i040
a/D = 1 Micromachines 13 00810 i041 Micromachines 13 00810 i042
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Du, X.; Ren, P.; Zheng, J. Predicting Milling Stability Based on Composite Cotes-Based and Simpson’s 3/8-Based Methods. Micromachines 2022, 13, 810. https://0-doi-org.brum.beds.ac.uk/10.3390/mi13050810

AMA Style

Du X, Ren P, Zheng J. Predicting Milling Stability Based on Composite Cotes-Based and Simpson’s 3/8-Based Methods. Micromachines. 2022; 13(5):810. https://0-doi-org.brum.beds.ac.uk/10.3390/mi13050810

Chicago/Turabian Style

Du, Xu, Pengfei Ren, and Junqiang Zheng. 2022. "Predicting Milling Stability Based on Composite Cotes-Based and Simpson’s 3/8-Based Methods" Micromachines 13, no. 5: 810. https://0-doi-org.brum.beds.ac.uk/10.3390/mi13050810

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