Chatter identification in thin-wall milling using an adaptive variational mode decomposition method combined with the decision tree model

Abstract

Chatter is prone to occur in thin-wall part milling process due to the low stiffness and damping of the workpiece. It roughens the machining surface, shortens the tool life, and thus should be detected and prevented. However, the multimode and time-varying dynamics of thin-wall parts produces nonstationary and multicomponent cutting signals, which makes it challenging to accurately identify the chatter occurrence. In this article, an effective chatter identification method based on adaptive variational mode decomposition and decision tree is presented to tackle this problem. The adaptive variational mode decomposition is used to adaptively decompose the raw acoustic signal into several subsignals, and the decision tree is used to automatically determine the chatter threshold. First, the convergence of multicenter-frequency signal processing is analyzed and is proved to be closely related to the accuracy of variational mode decomposition. Afterward, a criterion is set up to initialize the center frequencies of variational mode decomposition, on the basis of which an adaptive energy ratio-based method with good computational efficiency is presented to extract the frequencies of the main components from the raw signal. The initial center frequencies and the number of modes of variational mode decomposition are simultaneously obtained. As a result, the raw signal is adaptively decomposed into several subsignals, which contain its principal components. Then, the normalized energy and sample entropy of the subsignals are selected to establish the decision tree model for automatic chatter identification. Several milling tests on a thin-wall plate are carried out to verify the proposed method. The results show that chatter can be identified accurately and efficiently using the proposed method.

Keywords

Chatter identification thin-wall part milling adaptive variational mode decomposition decision tree

Introduction

Milling chatter is an ever-existing problem accompanying the machining of thin-wall parts, which dramatically reduces the machining accuracy and productivity.^1,2 In thin-wall milling process, the time-varying dynamics of the workpiece makes it difficult to build the cutting system model accurately, such that the cutting parameters may be unsuitable by using the analytic methods.^3–7 Paralleled by the development of sensor technology, the signal produced in the milling process is widely used to identify chatter. Hence, a sensor-based method is essential to identify and suppress chatter, which can enhance the surface finish and guarantee the stable cutting process.

The sensor-based chatter identification method involves signal collecting and processing, feature extraction, and states recognition. For signal collecting and processing, the acoustic signal with large bandwidth and low delay is widely adopted compared to other signals.⁸ However, the noise may be unavoidably collected into the signal by a microphone. Therefore, it is necessary to develop an effective signal processing method for denoising and enhancing the robustness of the identification system. To date, extensive efforts have been made to improve the accuracy and efficiency of signal processing, and the corresponding methods include time domain,^9,10 frequency domain,^11–14 and time-frequency domain analyses.^15–19 In time domain analysis, Lamraoui et al.⁹ applied the neural network and the entropy-based feature extraction method to detect chatter, but cannot easily extract the useful features from the real noisy signal. By the Fast Fourier Transform (FFT), the signal can be transformed into the frequency domain so that the features unidentified in time domain may appear clearly. And the noise also can be removed from the signal by denoising methods. Li et al.¹¹ investigated the properties of the chatter occurrence in the micromilling process, and reported that the amplitudes and the dominant frequencies of the signal shift when chatter occurs. Polli et al.¹² utilized different sensors simultaneously to analyze the chatter property and found that the chatter frequency dominates the distribution of the signal spectrum. Nevertheless, the results in the frequency domain by FFT may be incorrect due to the signal properties of nonlinearity and nonstationarity. To overcome this deficiency, several advanced signal processing methods in the time-frequency domain are proposed. Wavelet transform with the characteristic of adaptive time-frequency resolution shows good performance in revealing the local states in the cutting process. Cao et al.¹⁵ applied wavelet packets and Hilbert–Huang transform to process the acceleration signal for chatter identification. Based on wavelet transform, Yoon and Chin¹⁶ proposed a general chatter detection method, which can be used in other cutting processes. Moreover, chaotic theory can probably be used in predicting chatter. Koohestani et al.²⁰ presented an effective chatter detection method by using the linear regression model and adopted it for monitoring the cutting states and detecting chatter occurrence online.

As an adaptive time-frequency domain method, empirical mode decomposition (EMD) is proposed to decompose the signal into intrinsic mode functions (IMFs).²¹ Using EMD and Hilbert-Huang transform (HHT), Fu et al.¹⁹ proposed an online chatter detection method to distinguish the energy aggregation. In the study of Cao et al.,²² EEMD was used in the milling process and nonlinear dimensionless indicators were adopted to identify the onset of chatter vibration. However, the application of the EMD method is restricted by the shortage of mathematical theory and mode mixing. To enhance the robustness to sampling and noise, variational mode decomposition (VMD) is therefore developed²³ and outperforms EMD in terms of band-limit and mode mixing. Wang et al.²⁴ investigated the equivalent filter bank properties of VMD and reported its effectiveness in signal decomposition. The initial conditions were also found to significantly influence the accuracy of VMD. However, the majority of existing research studies^25,26 focus on optimizing the penalty parameter and the mode number without fully considering the role of the initial center frequencies (ICFs). Jiang et al.²⁷ analyzed the ICFs using simulated signals and found a U-shaped convergent phenomenon of center frequencies and reported that the convergent speed of VMD depends on the distribution of ICFs. Apart from the convergence speed, the range of center frequencies that limits the convergence is also important and discussed first in this article. A limiting bandwidth of the convergence is presented to optimize the ICFs, which improves the accuracy and efficiency of VMD significantly.

In this article, an effective chatter identification method is proposed in the milling of thin-walled parts. The acoustic signal is decomposed into submodes by the adaptive VMD (AVMD) method. Then, the normalized energy and the sample entropy, as sensitive indicators, are calculated and input to the model of chatter identification. In the conventional threshold-based chatter identification method, the value of the threshold is set by empirical judgment rather than the specific algorithm. Here, decision tree is utilized to compute the value of the threshold by using the CART algorithm.²⁸ The threshold is objective without personal judgment. Henceforth, this article is organized as follows: in section “Methodology of chatter identification for thin-wall part,” an AVMD method to preprocess a raw signal is described in detail. The normalized energy and sample entropy are then computed to build the decision tree. In section “Experimental setup,” the experiment setup is introduced. In section “Results and discussion,” the results and discussion are shown to verify the proposed method. The conclusion is drawn in section “Conclusion.”

Methodology of chatter identification for thin-wall part

Adaptive variational mode decomposition method

To enhance rapidity and accuracy in chatter identification, an AVMD method is proposed to preprocess the acoustic signal. It contains two aspects. First, the convergence of VMD is analyzed and a certain convergent scope of the mode bandwidth is revealed for multicenter frequency signal processing. Second, an adaptive energy ratio-based method for acquiring the frequencies of principle modes is developed to optimize the ICFs.

Brief introduction of VMD

As an adaptive and effective method for nonlinear and nonstationary signal processing, VMD aims to decompose a signal into narrowband modes. The bandwidth of each mode is estimated by the constrained variational problem as follows²³

\begin{matrix} \begin{matrix} \min \\ {u_{k}}, {ω_{k}} \end{matrix} {\sum_{k} ‖ \partial_{t} [{(δ (t) + \frac{j}{π t})}^{*} u_{k} (t)] e^{- j ω_{k} t} ‖_{2}^{2}} \\ s . t . \sum_{k} u_{k} = f \end{matrix}

(1)

where $u_{k}$ is the kth mode of the signal, $ω_{k}$ indicates the center frequency, and $δ$ is the Dirac distribution.

The expanded Lagrangian $L$ of equation (1) is written as

\begin{matrix} L ({u_{k}}, {ω_{k}}, λ) α \sum_{k} ‖ \partial_{t} [{(δ (t) + \frac{j}{π t})}^{*} u_{k} (t)] e^{- j ω_{k} t} ‖_{2}^{2} \\ + ‖ f (t) - \sum_{k} u_{k} (t) ‖_{2}^{2} + λ (t), f (t) - \sum_{k} u_{k} (t) \end{matrix}

(2)

By means of an alternate direction method of multipliers, the above quadratic optimization problem is solved

{\hat{u}}_{k}^{n + 1} (ω) = \frac{\hat{f} (ω) - \sum_{i \neq k} {\hat{u}}_{i} (ω) + \frac{\hat{λ} (ω)}{2}}{1 + 2 α {(ω - ω_{k})}^{2}}

(3)

Then, the center frequency $ω_{k}$ in equation (3) is reestimated through the gravity of the corresponding mode’s power spectrum as

ω_{k}^{n + 1} = \frac{\int_{0}^{\infty} ω {| {\hat{u}}_{k} (ω) |}^{2} d ω}{\int_{0}^{\infty} {| {\hat{u}}_{k} (ω) |}^{2} d ω}

(4)

And the terminal condition of the iterative algorithm is shown as

\sum_{k} \frac{| | {\hat{u}}_{k}^{n + 1} - {\hat{u}}_{k}^{n} | |_{2}^{2}}{| | {\hat{u}}_{k}^{n} | |_{2}^{2}} < ϵ

(5)

Center frequency convergence of VMD

For multiband signals, the convergence of VMD is analyzed first and then a certain convergent scope of the mode bandwidth is found. As a result, a criterion is introduced to set up the ICFs. A raw signal with two modes is given to analyze the center frequency convergence of VMD

f (t) = A_{1} \cos 2 π ω_{1} t + A_{2} \cos 2 π ω_{2} t + WGN

(6)

where $A_{1}$ and $A_{2}$ are the amplitudes of subcomponents $u_{k}$ , $ω_{1}$ and $ω_{2}$ are the natural frequencies, and WGN is the white Gaussian noise. The initial conditions of VMD are given in Table 1.

Table 1.

Initial conditions of VMD.

Signal component $u_{k}^{1}$	Center frequency $f_{k}^{1}$	Penalty parameter $α$	Number K	Lagrangian multiplier $λ$
$[0 0]$	$[f_{1}^{1} f_{2}^{1}]$	3000	2	0

VMD: variational mode decomposition.

In the following parts, the derivation process of the convergent range is expressed in detail. In order to match the domain of center frequency estimation in VMD, the time series $f (t)$ is transformed into Fourier domain

\hat{f} (ω) = {\begin{matrix} A_{1}, & ω = ω_{1} \\ A_{2}, & ω = ω_{2} \\ 0, & others \end{matrix}

(7)

According to equations (3) and (4), the mode $u_{1}^{2}$ and the center frequency $f_{1}^{2}$ are calculated under the above initial conditions

u_{1}^{2} = {\begin{matrix} \frac{A_{1}}{1 + 2 α {(ω - f_{1}^{1})}^{2}}, & ω = ω_{1} \\ \frac{A_{2}}{1 + 2 α {(ω - f_{1}^{1})}^{2}}, & ω = ω_{2} \\ 0, & others \end{matrix}

(8)

f_{1}^{2} = \frac{ω_{1} {| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{1}^{1})}^{2}} |}^{2} + ω_{2} {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{1}^{1})}^{2}} |}^{2}}{{| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{1}^{1})}^{2}} |}^{2} + {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{1}^{1})}^{2}} |}^{2}}

(9)

According to equation (3), the new mode $u_{1}^{n + 1}$ is associated with the mode $u_{2}^{n}$ and the center frequency $f_{1}^{n}$ of the last step. Therefore, the mode $u_{2}^{2}$ should be calculated to get the mode $u_{1}^{3}$ , as follows

u_{2}^{2} = \frac{\hat{f} (ω) - u_{1}^{2} (ω)}{1 + 2 α {(ω - f_{2}^{1})}^{2}}

(10)

And the numerator of equation (10) is shown as

\hat{f} (ω) - u_{1}^{2} (ω) = {\begin{matrix} A_{1} - \frac{A_{1}}{1 + 2 α {(ω - f_{1}^{1})}^{2}} = A_{1} \frac{2 α {(ω_{1} - f_{1}^{1})}^{2}}{1 + 2 α {(ω_{1} - f_{1}^{1})}^{2}}, (ω = ω_{1}) \\ A_{2} - \frac{A_{2}}{1 + 2 α {(ω - f_{1}^{1})}^{2}} = A_{2} \frac{2 α {(ω_{2} - f_{1}^{1})}^{2}}{1 + 2 α {(ω_{2} - f_{1}^{1})}^{2}}, (ω = ω_{2}) \\ 0, others \end{matrix}

(11)

Under the premise of the proper penalty parameter $α$ , the data fidelity can be guaranteed.²³ Hence, a reasonable approximation is applied to simplify equation (11)

\frac{2 α {(ω_{1} - f_{1}^{1})}^{2}}{1 + 2 α {(ω_{1} - f_{1}^{1})}^{2}} \approx {\begin{matrix} 1, f_{1}^{1} \neq ω_{1} \\ 0, f_{1}^{1} = ω_{1} \end{matrix}

(12)

With this approximation, the mode $u_{2}^{2}$ and the center frequency $f_{2}^{2}$ are simplified as

u_{2}^{2} = {\begin{matrix} \frac{A_{1}}{1 + 2 α {(ω - f_{2}^{1})}^{2}}, (ω = ω_{1}, f_{1}^{1} \neq ω_{1}) \\ \frac{A_{2}}{1 + 2 α {(ω - f_{2}^{1})}^{2}}, (ω = ω_{2}, f_{1}^{1} \neq ω_{2}) \\ 0, others \end{matrix}

(13)

f_{2}^{2} = \frac{ω_{1} {| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{2}^{1})}^{2}} |}^{2} + ω_{2} {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{2}^{1})}^{2}} |}^{2}}{{| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{2}^{1})}^{2}} |}^{2} + {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{2}^{1})}^{2}} |}^{2}}

(14)

In accordance with the above formulas, the latest mode $u_{1}^{3}$ and its center frequency $f_{1}^{3}$ are expressed according to equation (3)

u_{1}^{3} = {\begin{matrix} \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{1}^{2})}^{2}}, (ω = ω_{1}, f_{2}^{1} \neq ω_{1}) \\ \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{1}^{2})}^{2}}, (ω = ω_{2}, f_{2}^{1} \neq ω_{2}) \\ 0, others \end{matrix}

(15)

f_{1}^{3} = \frac{ω_{1} {| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{1}^{2})}^{2}} |}^{2} + ω_{2} {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{1}^{2})}^{2}} |}^{2}}{{| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{1}^{2})}^{2}} |}^{2} + {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{1}^{2})}^{2}} |}^{2}}

(16)

The center frequency is gradually close to the first intrinsic frequency of the initial signal under the reasonable approximation in equation (12), and thus the mode $u_{1}$ is obtained at the end of iteration. Considering the principle of VMD, under the same iteration process, the next mode $u_{2}$ is calculated by subtracting the first mode $u_{1}$ from the initial signal. Hence, the general formula for center frequency convergence of the two-mode signal is deduced

f_{k}^{n + 1} = \frac{ω_{1} {| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{k}^{n})}^{2}} |}^{2} + ω_{2} {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{k}^{n})}^{2}} |}^{2}}{{| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{k}^{n})}^{2}} |}^{2} + {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{k}^{n})}^{2}} |}^{2}}

(17)

Based on the above analysis of a signal with two modes, it is easy to expand the above derivation to the multicomponent signal. Thus, the general formula of m components is shown as follows

{\bar{f}}_{k}^{n + 1} = \frac{ω_{1} {| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{k}^{n})}^{2}} |}^{2} + ω_{2} {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{k}^{n})}^{2}} |}^{2} + \dots + ω_{m} {| \frac{A_{m}}{1 + 2 α {(ω_{m} - f_{k}^{n})}^{2}} |}^{2}}{{| \frac{A_{1}}{1 + 2 α {(ω_{1} - f_{k}^{n})}^{2}} |}^{2} + {| \frac{A_{2}}{1 + 2 α {(ω_{2} - f_{k}^{n})}^{2}} |}^{2} + \dots + {| \frac{A_{m}}{1 + 2 α {(ω_{m} - f_{k}^{n})}^{2}} |}^{2}}

(18)

where ${\bar{f}}_{k}^{n}$ is the latest center frequency of the mode $u_{k}$ and $ω_{m}$ is the natural frequency of its corresponding mode.

According to equations (17) and (18), the convergent center frequencies are related to the amplitude of each intrinsic mode and ICFs obviously. To verify this result, we give the initial center frequencies in the range of 0–128 Hz, and then calculate the center frequencies according to equation (17). At the end of iteration, the center frequency may converge to the natural frequency of each mode. Figure 1 shows that the center frequency is gradually approximating to the natural frequency of the mode in a limiting bandwidth. And the amplitude of each mode influences the final center frequency, which corresponds to the gravity of power spectral density. The influence makes the curves in Figure 1(a) and (b) appear slightly different.

Figure 1.

The limit band of center frequency convergence of two-mode signal: (a) f = 64 Hz, A = 2, and f = 128 Hz, A = 4. (b) f = 64 Hz, A = 4, and f = 128 Hz, A = 2.

It is observed that the critical frequency $f_{a}$ determines the distribution of initial frequency to guarantee the result of frequency convergence according to Figure 1(a). Therefore, the initial frequencies should be set on both sides of the critical frequency as well as close to natural frequencies. In Figure 2, for a signal with three modes, the convergence of VMD is further illustrated and the influence of amplitude is also concerned. It investigates that the distribution of initial center frequencies should be set in the limiting bandwidths around the frequencies of the intrinsic modes, so as to guarantee the convergence.

Figure 2.

The limit band of center frequency convergence of a three-mode signal: (a) f = 16 Hz, A = 1, f = 32 Hz, A = 2, and f = 64 Hz, A = 4 and (b) f = 16 Hz, A = 1, f = 32 Hz, A = 4, and f = 64 Hz, A = 2.

Nevertheless, considering the spectral leakage of FFT, the gravity of the power density spectrum will shift a little. It causes a little error in the range of setting ICFs. In order to ensure the convergence, the limit band of ICFs is constrained within a little wide range properly

{f_{i}^{1} | f_{N, i} - r_{i} < f_{i}^{1} < f_{N, i} + r_{i}, i = N^{*}}

(19)

where $f_{i}^{1}$ is the initial center frequency of the i th mode, $f_{N, i}$ is the natural frequency of the i th mode, and $r_{i}$ is the small positive to constrain the range of $f_{i}^{1}$ , which ensures $f_{i}^{1}$ equal to or less than the critical frequency.

Frequency acquisition of signal main components

As for the VMD method, ICFs influence the result seriously although the other parameters are fine-tuning. In this section, an adaptive energy ratio-based sifting method is presented to find out the frequencies of the signal main components. As a consequence, the ICFs and the number of modes are simultaneously obtained.

First, an acoustic signal $x (t)$ collected from the thin-wall milling process is transferred into Fourier domain, and the PSD is calculated

S (f) = \frac{1}{2 π N} {| X (w) |}^{2}

(20)

where N represents the length of signal and $X (w)$ is the Fourier transform of input signal $x (t)$ .

In order to obtain the frequencies of the main components, two sifting procedures are applied. According to the curve of PSD in spectral domain, a series of local maximum is computed with the interval $d_{1}$ so that much useful information can be preserved.

The series of PSD is divided into frames as $X_{k} (f)$ , $k = 1, 2, 3, \dots, N$ . Thus, the local maximum $S_{\max, k} (f)$ and the corresponding frequency $f_{k}$ are gained

S_{\max, k} (f) = max {X_{k} (f)}

(21)

f_{k} = \arg {S_{\max, k} (f)}

(22)

Toward the curve of local maximum series, the values beyond the average are sifted out and its corresponding series are defined as $S_{s t_{1}, k_{1}} (f)$

S_{s t_{1}, k_{1}} (f) = \arg {S_{\max, k} (f) \geq mean (S_{\max, k} (f))}

(23)

f_{s, k 1} = \arg {S_{s t_{1}, k_{1}} (f)}

(24)

Due to the reduction of data points, the curve appears as several straight parts that become obstacles to search for peak values. To overcome this case, a cubic spline is used to fit the obtained series $S_{s t_{1}, k_{1}} (f)$ . As a result, the new curve, defined as $S_{s t_{2}, k_{2}} (f)$ , becomes smoother at cuspidal points and the uniqueness of the local maximum value is ensured.

In the second sifting procedure, the frequencies of the main components are obtained by using the energy ratio-based iterative sifting algorithm. The energy ratio (i.e. as the threshold) is set as 0.8, and the initial sifting value is set as the midvalue $V_{i}$ of the sifted curve $S_{s t_{2}, k_{2}} (f)$ . Peak points of the curve are found out, and then the targeted points, larger than mid-value, are selected. Subsequently, the energy $e_{p}$ around all sifted points is added up and the energy ratio is calculated. If the energy ratio exceeds the threshold, the algorithm goes to end; otherwise, the initial sifting value is reset by multiplying the adjusting factor. Finally, the principle frequency $f_{e}$ and the corresponding amplitude $S_{e} (f)$ are obtained. The process of the energy ratio-based iterative sifting algorithm is summarized as follows:

Step 1	Initialization of the second sifting procedure
	The energy ratio: $e_{r} = 0$
	The corresponding threshold: $threshold = 0.8$ .
	The initial sifting value:
	$V_{i} = 0.5 \cdot (max (S_{s t_{2}, k_{2}}) - min (S_{s t_{2}, k_{2}}))$ .
	The interval width h
	The total energy: $e_{t} = \sum^{(S_{s t_{2}, k_{2}})}$ .
Step 2	Implementation of the second procedure
	While $e_{r} \leq threshold$
	$j = 1$
	for i = 2 to length $(S_{s t_{2}, k_{2}})$ − 1
	if $S_{s t_{2}, i} (f) > max (S_{s t_{2}, i - 1} (f), S_{s t_{2}, i + 1} (f))$
	$and S_{s t_{2}, i} (f) \geq V_{i}$
	$S_{e, j} (f) = S_{s t_{2}, i} (f)$
	$f_{e, j} = \arg (S_{e, j} (f))$
	$j = j + 1$
	end
	end
	$e_{p} = 0$
	$f_{e, j}^{'} = {sort (S_{e, j} (f))}$
	$n = length (S_{e, j} (f))$
	$k = 0$
	for i = 1: n(local accumulation around the $f_{e, j}^{'}$ )
	$e_{p} = e_{p} + \sum_{f}^{f_{e, i \pm h}^{'}} S_{e, i} (f)$
	$e_{r} = \frac{e_{p}}{e_{t}}$
	if $e_{r} \geq threshold$
	$f_{e, k} = f_{e, i}^{'}$
	$k = k + 1$
	$break;$
	end
	end
	$V_{i} = {V_{i}}^{*} threshold$
	end
Step 3	The main frequencies $f_{e}$ and the relevant amplitude $S_{e}$ of the raw signal are acquired.

The frequencies of the signal main components can be obtained after the above sifting procedures with the appropriate threshold and interval width. It is automatic and time-saving that has great potential for online chatter detection. Through the preprocessor for VMD, the initial conditions that consist of ICFs and the number of modes are calculated. Hence, an adaptive center frequencies initialization-based VMD (AVMD) method is conducted. The subsignals with intrinsic narrowband and maximal energy are obtained through the AVMD method.

Decision tree construction with sensitive indicators

The energy is widely used in chatter identification, because it has an apparent change when chatter occurs. As a nondimensional value, the sample entropy is sensitive to signal complexity in time series and is widely adopted in the fault diagnosis of rotating machine recently.²⁵ Based on the two indicators, the overall trend and local complexity are revealed concurrently.

Through the AVMD method, several subsignals are successfully limited into narrowbands. The energy of the nth subsignal $x_{n} (t)$ with N frames is computed as

E_{n} (t) = X_{rms} (t) = \sqrt{\frac{\sum_{t = 1}^{N} x_{n} {(t)}^{2}}{N}}

(25)

While the energy signal contains much noise, an exponentially weighted moving average method is adopted to filter out the noise and enhance the changing tendency. This kind of filter is sensitive to the latest input data that can guarantee the real-time property. The filtered energy signal is expressed as $E_{n}^{f} (t)$ . In order to make the approach general, the filtered energy signal is normalized as shown in Figure 3

E_{n}^{norm} (t) = \frac{E_{n}^{f} (t) - μ_{n}}{s_{n}}

(26)

where $μ_{n} = \frac{1}{N} \sum_{i = 1}^{N} E_{n, i}^{f} (t)$ , $s_{n} = \sqrt{\frac{1}{N - 1} \sum_{i = 1}^{N} {(E_{n, i}^{f} (t) - μ_{n})}^{2}}$

Figure 3.

Normalized energy of subsignals.

The sample entropy is regarded as an assessment of the complexity of a time series by measuring the probability of generating new patterns.²⁹ It is independent of the length of data. Aiming at the subsignal series $x_{n} (t)$ with N frames that decomposed from the actual signal by the AVMD method, the sample entropy is calculated as follows:

(1) Construct a series of vectors T with $N - m + 1$ dimension

T = {X_{m} (1), \dots, X_{m} (i), \dots, X_{m} (N - m + 1)}

where $X_{m} (i) = [x (i), x (i + 1), \dots, x (i + m - 1)]$ , $1 \leq i \leq N - m + 1$ .

(2) Define the distance between $X_{m} (i)$ and $X_{m} (j)$

d_{ij} = max_{1 \leq k \leq m - 1} (| x (i + k) - x (j + k) |)

(27)

(3) Aiming at $X_{m} (i)$ , calculate the number of data points that $d_{ij} \leq r$ and is denoted as $B_{i}$ . Then, compute $B_{i}^{m} (r)$

B_{i}^{m} (r) = \frac{1}{N - m - 1} B_{i}

(28)

(4) Define $B^{(m)} (r)$ as

B^{(m)} (r) = \frac{1}{N - m} \sum_{i = 1}^{N - m} B_{i}^{m} (r)

(29)

(5) Add the dimension to m+ 1 and repeat the above steps to calculate $B^{(m + 1)} (r)$ , then rewritten as $A^{(m)} (r)$ .

(6) The sample entropy of a subsignal $x_{n} (t)$ is defined as

SampE n_{n} (m, r, N) = - \ln [\frac{A^{(m)} (r)}{B^{(m)} (r)}]

(30)

For each subsignal, the energy and sample entropy of each frame are computed. The global trend of the signal is shown by energy, while the local complexity is expressed by the sample entropy. Nevertheless, it is necessary to filter out the chatter-independent components. According to the trend of energy and sample entropy, the most relevant component to cutting state is retained to identify the location of chatter occurrence and shown in Figure 4.

Figure 4.

Energy and sample entropy of subsignal 2.

For chatter identification, the decision tree model is built by the CART algorithm,²⁸ rather than the empirical model constructed from the extensive experiments. To enhance accuracy of the CART model, the cutting states are divided into two categories: chatter and nonchatter. The training dataset is constructed by double sensitive indicators, and the corresponding feature set is labeled through supervising learning. Then, the CART algorithm is applied to build the decision tree for cutting states recognition in the milling process, and the details are shown in Figure 6.

Experimental setup

To demonstrate the proposed method, milling tests at multipaths with different axial cutting depths are carried out on a five-axis CNC machine shown in Figure 5. The axial depth of cut is set as 0.2, 0.5, 0.8, and 1 mm, while the radial cutting depth is 0.8 mm. A solid carbide end mill cutting tool with 16 mm diameter and 45° helix angle is used. In addition, the thin-walled work piece is made of Al 6061 with the dimension: 80 × 40 × 3 mm³. A microphone GRAS 46AE is applied to collect the acoustic signal. The sample frequency is set as 25,600 Hz. The spindle speed of the milling experiments is 5000 r/min.

Figure 5.

Experimental setup.

Results and discussion

To verify the proposed method, the cutting states in thin-wall milling were identified by following the procedure shown in Figure 6. The acoustic signal with axial depth of cut 0.8 mm is shown in Figure 7, which contains three categories of cutting states (i.e. air cut, stable, and chatter states). To demonstrate the property of time-varying and multimode in thin-wall milling process, the frequency spectrums at the regions (a), (b), and (c) are calculated as shown in Figure 8, which represents the different cutting states distributed in the corresponding location. Note that the spectral distribution and its amplitude change significantly when chatter happens. The frequencies of the signal main components are expressed as $f = [1500; 4482; 8501; 12, 407]$ close to the natural frequencies of the work piece. To verify the proposed AVMD method, the adaptive energy ratio-based signal sifting algorithm is used to find the frequencies of the main components automatically, and the result is shown in Figure 9. The main frequencies in Figure 9 depicted by red dots are obtained as $f = [1500; 4482; 8501; 12, 407]$ , which is consistent with the frequencies of the spectral analysis in Figure 8. The energy around the frequencies of the main components closely accounts for 80% of the total energy. Hence, the subsignals can represent the main components of the original signal. Simultaneously, the number of mode K and ICFs for VMD are gained.

Figure 6.

Flowchart of the proposed chatter identification method.

Figure 7.

The acoustic signal collected from thin-wall milling.

Figure 8.

The results of FFT at the corresponding regions (a), (b), and (c) in Figure 7.

Figure 9.

The results of the adaptive energy ratio-based sifting method.

Next, the VMD method with the optimal initial conditions is used to process the vibration signal. The number of mode ( $K = 4$ ) and the initial center frequencies ( $ICFs = [1500; 4482; 8501; 12, 407]$ ) are applied. In Figure 10, the subsignals and the results of FFT are presented. Each subsignal is limited into a narrowband successfully. By adding the subsignals, the reconstructed signal is built and the reconstructed error is computed. Figure 11 shows the results of reconstruction accuracy with different initialization methods: the proposed method and the uniform distribution method. Combined with the adaptive center frequencies initialization method, the accuracy of VMD has been improved greatly when chatter occurs. At chatter regions, the superiority of the proposed method is prominent and the absolute error is well limited in allowable error according to Figure 11(b). In comparison with the fundamental initialization method, the AVMD method is demonstrated to be feasible and effective for the signal that is composed of multiple modes.

Figure 10.

Subsignals obtained from the AVMD method and FFT.

Figure 11.

Results of the VMD method with different initial center frequencies.

To generate the decision tree, the normalized energy and the sample entropy of the subsignals are calculated. However, useless information are found along with the energy curves, which need to be eliminated. Figure 4 shows the selected curves of the energy and the sample entropy. It is obvious that the cutting states can be classified into chatter and nonchatter. In consideration of surface finish, the feature set is labeled under the supervised learning. The training dataset of the decision tree is set by the energy curve and sample entropy. Tree pruning, as a vital step, can enhance the generalization of decision tree to fit other conditions in a small scope.

After tree pruning, the decision tree is obtained to identify the cutting states into two aspects: chatter and nonchatter. At the root of tree model, the sample entropy is selected as classification criteria by the CART algorithm. Meanwhile, the normalized energy is calculated to classify the states at other nodes of the branches. The value of each threshold is computed and shown in Figure 12. Compared with the methods that determine threshold empirically, the thresholds obtained by the decision tree model are more reasonable, although it is suboptimal. To verify the accuracy of the decision tree model, another cutting test with a depth of cut of 0.8 mm is carried out that also contains the above states (i.e. chatter and nonchatter). With the proposed AVMD method, the new signal is decomposed into several modes. Then, the energy and its corresponding sample entropy after selecting are shown in Figure 13. Although there is a little difference between the two signals of two tests, the result of chatter identification is acceptable with a high percentage of accuracy. The rate of accuracy is 92.42% according to the statistics of cutting states. For the approaches of identifying chatter by threshold, the value of threshold is usually based on personal experience. Due to the subjective factors, it will not be optimal. Conversely, the threshold in the decision tree model is calculated by the CART algorithm objectively, where the cutting states along with tool path is fully taken into account. And the decision tree, as a suboptimal model, can make good prediction about unknown data.

Figure 12.

The decision tree model.

Figure 13.

The energy and sample entropy of the subsignal in another test.

Conclusion

An effective method by means of AVMD and decision tree is proposed to identify chatter in thin-wall milling. The measured acoustic signal that contains multiple components is difficult to handle for VMD. Therefore, an adaptive energy ratio-based sifting method is presented to improve initialization conditions of VMD. As a result, the above signal can be self-adaptively decomposed into submodes with the expected narrowband. Unlike the chatter identification method with an empirical threshold, decision tree is established to calculate the chatter threshold by the CART algorithm objectively. The results of milling tests validate that chatter can be identified at the allowable scale with the proposed model. In addition, the proposed method has great potential for application in online chatter detection. Furthermore, the identification of chatter should be directly mapped to the quality of the machined surface, rather than the derived threshold.

Footnotes

Appendix 1 Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research is supported by the National Key Research and Development Program of China (2018YFA0704603), the National Natural Science Foundation of China (91948203), the China Postdoctoral Science Foundation (BX20190054), and the Natural Science Foundation of Liaoning Province (2019-BS-045).

ORCID iD

Yuwen Sun

References

Altintas

Manufacturing automation: metal cutting mechanics, machine tool vibrations and CNC design. New York: Cambridge University Press, 2012.

Quintana

Ciurana

Chatter in machining processes: a review. Int J Mach Tools Manuf 2011; 51(5): 363–376.

Tang

Liu

ZQ.

Three-dimensional stability lobe and maximum material removal rate in end milling of thin-walled plate. Int J Adv Manuf Technol 2009; 43(1–2): 33–39.

Cao

Zhang

Huang

, et al. An improved semi-analytical approach for modeling of process damping in orthogonal cutting considering cutting edge radius. Proc IMechE Part B: J Engineering Manufacture 2020; 234(3): 641–653.

Sun

Jiang

SL.

Predictive modeling of chatter stability considering force-induced deformation effect in milling thin-walled parts. Int J Mach Tools Manuf 2018; 135: 38–52.

Yan

Luo

Tang

Multi-axis variable depth-of-cut machining of thin-walled workpieces based on the workpiece deflection constraint. Comput Des 2018; 100: 14–29.

Jiang

Sun

YW.

Chatter stability and surface location error predictions in milling with mode coupling and process damping. Proc IMechE Part B: J Engineering Manufacture 2019; 233(3): 686–698.

Delio

Tlusty

Smith

Use of audio signals for chatter detection and control. J Manuf Sci Eng Trans ASME 1992; 114(2): 146–157.

Lamraoui

Barakat

Thomas

, et al. Chatter detection in milling machines by neural network classification and feature selection. J Vib Cont 2015; 21(7): 1251–1266.

10.

Huang

Sun

, et al. Vibration analysis in milling titanium alloy based on signal processing of cutting force. Int J Adv Manuf Technol 2013; 64(5–8): 613–621.

11.

Jing

Wang

Detection and analysis of chatter occurrence in micro-milling process. Proc IMechE Part B: J Engineering Manufacture 2014; 228(11): 1359–1371.

12.

Polli

Weingaertner

Schroeter

, et al. Analysis of high-speed milling dynamic stability through sound pressure, machining force and tool displacement measurements. Proc IMechE Part B: J Engineering Manufacture 2012; 226(11): 1774–1783.

13.

Cao

Zhang

Chen

XF.

The concept and progress of intelligent spindles: a review. Int J Mach Tools Manuf 2017; 112: 21–52.

14.

Kang

Feng

, et al. Research on chatter prediction and monitor based on DHMM pattern recognition theory. Proc IEEE Int Conf Autom Logist ICAL 2007 2007(1): 1368–1372.

15.

Cao

Lei

ZJ.

Chatter identification in end milling process using wavelet packets and Hilbert-Huang transform. Int J Mach Tools Manuf 2013; 69: 11–19.

16.

Yoon

Chin

DH.

Cutting force monitoring in the end milling operation for chatter detection. Proc IMechE Part B: J Engineering Manufacture 2005; 219(6): 455–465.

17.

Yuan

Jing

, et al. Chatter detection based on wavelet coherence functions in micro-end-milling processes. Proc IMechE Part B: J Engineering Manufacture 2019; 233(9): 1934–1945.

18.

Lei

Lin

, et al. A review on empirical mode decomposition in fault diagnosis of rotating machinery. Mech Syst Signal Process 2013; 35(1–2): 108–126.

19.

Zhang

Zhou

, et al. Timely online chatter detection in end milling process. Mech Syst Signal Process 2016; 75: 668–688.

20.

Koohestani

JPT

Yang

Application of linear regression model on chatter threshold delineation. Proc IMechE Part B: J Engineering Manufacture 2015; 229(8): 1337–1347.

21.

Huang

Shen

Long

, et al. The empirical mode decomposition and the Hubert spectrum for nonlinear and non-stationary time series analysis. Proc R Soc A Math Phys Eng Sci 1998; 454(1971): 903–995.

22.

Cao

Zhou

Chen

XF.

Chatter identification in end milling process based on EEMD and nonlinear dimensionless indicators. Int J Mach Tools Manuf 2015; 92: 52–59.

23.

Dragomiretskiy

Zosso

Variational mode decomposition. IEEE Trans Signal Process 2014; 62(3): 531–544.

24.

Wang

Markert

Xiang

, et al. Research on variational mode decomposition and its application in detecting rub-impact fault of the rotor system. Mech Syst Signal Process 2015; 60: 243–251.

25.

Yang

Wang

Dong

, et al. Early chatter identification based on an optimized variational mode decomposition. Mech Syst Signal Process 2019; 115: 238–254.

26.

Shi

Yang

WX.

Precise feature extraction from wind turbine condition monitoring signals by using optimised variational mode decomposition. IET Renew Power Gener 2017; 11(3): 245–252.

27.

Jiang

Shen

Shi

, et al. Initial center frequency-guided VMD for fault diagnosis of rotating machines. J Sound Vib 2018; 435: 36–55.

28.

Breiman

Classification and regression trees. New York: Routledge, 2017.

29.

Yin

Shang

Feng

GC.

Modified multiscale cross-sample entropy for complex time series. Appl Math Comput 2016; 289: 98–110.