OCSVM與LOF演算法用於智慧主軸之即時監測與顫振分析抑制

(1)

國立臺灣大學工學院機械工程學研究所碩士論文

Graduate Institute of Mechanical Engineering College of Engineering

National Taiwan University Master Thesis

OCSVM 與 LOF 演算法用於智慧主軸之即時監測與顫振分析抑制

Real-time Chatter Detection, Analysis and Suppression for Intelligent Spindles Based on One-class Support Vector Machine

and Local Outlier Factor

姚泳辰 Yung-Chen Yao

指導教授：施文彬博士、劉建豪博士 Advisor: Wen-Pin Shih, Ph.D.

Chien-Hao Liu, Ph.D.

中華民國 107 年 6 月

June, 2018

(2)

(3)

誌謝

能完成這篇論文，首先，我要特別感謝我的指導教授施文彬老師，施老師總能在討論中一次次點出我的不足之處，並給予適當建議，同時也時常鼓勵我學習不同領域的知識，相信這對我不論在研究上或是未來發展上皆有很深的影響，真的非常感謝施老師。

其次，要感謝劉建豪老師，在施老師去美國的這一年給予我的指導與幫助，協助我尋找研究所需的資源，並給予論文建議。也謝謝施博仁老師，總能點出一些理論與實驗上的細節，使我的研究更完善。另外，感謝應力所的張培仁教授能引薦我至工研院做相關研究，感謝工研院的林錦德博士及工程師陳韻巧小姐，提供機台讓我做實驗，並總能協助我突破在實驗上遇到的困難及給予建議。

感謝實驗室學長瑞鴻、彥廷、明華、胤禎，學長們是實驗室的好榜樣，讓實驗室的學習氣氛非常良好，也總是能在我遇到瓶頸時提供我不同角度的看法。感謝同屆的同學善謙、星宇、家倫、建璋、尚軒、歆儒，在兩年期間一同學習，一起度過各種難關與喜怒哀樂。感謝實驗室學弟宇軒、炫奇、健宏，和你們交流總有出乎意料外的收穫，有你們在也為實驗室注入不少活力。

另外，感謝應力所的家銘與品蓉，做為不同子計畫的負責人時常能給我與不同方面的建議與鼓勵。感謝廖運炫教授、宗昇學長與譽升學弟，能夠借用機台讓我做實驗，也給予很多操作機台上的協助。

最後，必須感謝我的父母、家人、女友馨云，時常給予我鼓勵與陪伴，使我更有動力走完這段學業歷程，謝謝。

謝謝所有人的幫助，有你們我才能完成這篇論文，謝謝。

(4)

中文摘要

隨著工業4.0 的時代到來，對智慧工具機的需求與日俱增。智慧工具機強調機

台本身具有自我監測、判斷與修復的能力，而工具機主軸切削顫振即為其中一項重要議題。顫振為主軸切削過程中所引發的自激式振動，顫振發生時會破壞工件表面精度且加速刀具磨耗，嚴重甚至會造成斷刀等現象。過去已有多位學者提出顫振偵測方法，然而大多數研究僅能分辨出成長完全的顫振訊號，無法分辨轉換狀態資料，

且顫振偵測的閾值設定也多採用經驗法則，較難達到即時偵測與抑制的效果。因此，

本論文提出一套智慧主軸顫振抑制模組，透過建立切削振動力學模型、小波包能量熵、異常偵測演算法及顫振抑制策略達到即時偵測及抑制顫振的效果。首先，建立切削力學模型以推導主軸顫振頻率，並繪製了三維系統穩定性圖。此模型考慮了進刀速率所造成的時間延遲影響，使系統穩定性圖增加了進刀速率的維度。其次，以小波包轉換計算相對小波包能量熵做為偵測顫振的特徵。小波包轉換對高頻有良好的靈敏度，可使顫振更早被預測。另外，以一類支持向量機與局部密度因子分辨顫振資料並比較兩種演算法的結果，實驗結果顯示局部密度因子代價函數為較低

的 0.0684，約略可視為有 93.16%的分辨率，有較好的分類效果。最後，基於系統

穩定性圖，制定顫振抑制策略，優先確保加工效率的提升，同時達到抑制顫振的效果。本論文從預處理至即時監控提出一套完整對策，成功開發出偵測顫振轉換狀態資料，並抑制顫振發生之系統。

關鍵字：智慧主軸、顫振抑制、即時監控、系統穩定性圖、小波包能量熵、一類支持向量機、局部密度因子

(5)

ABSTRACT

With the advent of industry 4.0 era, the demand for intelligent machine tools is increasing. Intelligent machine tools emphasize that machine has the capability of self- monitoring, self-abnormal diagnosis and self-repairing. One of the significant issues of intelligent machine tools is chattering of spindle. Chatter is a self-excited vibratio n generated during cutting process. When it occurs, the accuracy of the workpiece surface is destroyed and the tool wear out rapidly. Even worse, it may cause the cutting tool broken. Different chatter detection methods have been proposed in the past. However, most researches could only identify fully developed chatter data rather than transitio n data. Moreover, setting threshold values of chatter often relied on experiences and experiments. Therefore, this thesis presented an intelligent-spindle chatter-suppressio n module to achieve real-time detections and suppressions of chatter through the establishment of spindle vibration model, relative wavelet packet energy entropy, anomaly detection algorithm and chatter suppression strategy. First, the spindle vibratio n model was established to derive the chatter frequency of spindle and construct the three- dimensional stability lobe diagram. This model took the variable time delay caused by feed rate into account, which increase the dimension of feed rate in stability lobe diagram.

Second, the relative wavelet packet energy entropy calculated by wavelet packet transform was used as the feature of chatter detection. Wavelet packet transform had good sensitivity in high frequency which lead to early prediction of chatter. In addition, one- class support vector machine (OCSVM) and local outlier factor (LOF) were used to identify transition data. Based on the experimental results obtained from the two algorithms, local outlier factor demonstrated 0.0684 cost, which could be considered as

(6)

93.16% accuracy based on the weight of each class and was better than the one-class support vector machine. Last, based on the stability lobe diagram, a chatter suppression strategy was developed to improve the efficiency of machining and achieve effective chatter suppression. In conclusion, this thesis presented an intelligent chatter-suppressed system composed of preprocessing, real-time detection, and immediate suppressions by detecting chatter in transition state and suppressing chattering automatically.

Keywords: intelligent spindle, chatter suppression, real-time detection, stability lobe diagram, relative wavelet packet energy entropy, one-class support vector machine, local outlier factor

(7)

SYMBOL TABLE

M modal mass matrix C damping matrix K stiffness matrix q displacement of tool F cutting force

m x mass in x direction my mass in y direction c x damping in x direction cy damping in y direction k x stiffness in x direction ky stiffness in y direction t time

x displacement in x direction y displacement in y direction

,

F i x cutting force of i^th tooth in x direction

,y

F i cutting force of i^th tooth in y direction

,

F i t cutting force of i^th tooth in tangential direction

,

F i n cutting force of i^th tooth in normal direction

i angular position of i^th tooth Kt tangential cutting force coefficient Kn normal cutting force coefficient

(8)

XF exponent terms

gi contact function of i^th tooth ap axial depth of cut

hi chip thickness

st enter angle

ex exit angle f chip load



time delay T time period

R radius of tool N number of tooth

 spindle speed

aD radius depth of cut ratio Q0 periodic orbit

X perturbation

t time step

,

lrj i residual ratio of the time delay Φst transition matrix

 eigenvalue of transition matrix

fc basic chatter frequency f tpe tooth path frequency f H Hopf bifurcation frequency

(9)

fpd period-doubling frequency

fPD period-doubling bifurcation frequency f T cyclic-fold bifurcation frequency

 mother wavelet a scale factor b shift factor

J maximum level of decomposition K length of filter

signal

f signal function

 scaling function h low-pass filter g high-pass filter

W wavelet packet function

, , j n k

c wavelet packet coefficient

,

Ej n wavelet packet energy

Etot total energy of all frequency band

,

pj n relative energy SE j Shannon entropy

,

RWPEE j n relative wavelet packet energy entropy w weight of support vector machine

b bias of support vector machine xn input data position

(10)

yn the correctness of data HP optimal margin of hyperplane

 feature transform function

z n data in feature space d feature space

d original space

n first Lagrange Multipliers

n second Lagrange Multipliers L Lagrange function

K kernel function

n slack variable

C the constant controls the trade-off between large margin and margin violation

 the distance between hyperplane and origin in feature space



control the effect of slack variable -distance( )

k p the k-th nearest neighbor from the point p

k( )

N p the number of data set in the k-distance( )p area reach-dist_k reachability-distance

lrd_k local reachability density LOF_k local outlier factor

T p true positive T n true negative

F p false positive

(11)

Fn false negative PR precision RE recall

MA missing alarm rate A

F false alarm rate

C st cost function

CM weight of missing alarm rate CF weight of false alarm rate F x mean cutting force in x direction Fy mean cutting force in y direction r the ratio between K_t and K_n

f n natural frequency

 damping ratio

(12)

LIST OF FIGURES

Fig 1.1 An integrated concept of intelligent spindles [1]...2

Fig 1.2 Measured vibration signal and the corresponding workpiece surface [3]. ...3

Fig 1.3 Stability lobe diagram compared with experimentally performed cuts [7]. ...4

Fig 1.4 Stability lobe diagrams obtained from different methods [11]. ...5

Fig 1.5 FFT analysis on force signatures in y direction [14]. ...6

Fig 1.6 Behavior of the detail coefficients wavelet transformed in levels 1 to 9 for cutting force [17]. ...7

Fig 1.7 Illustration of chatter detection based on k-means clustering. Cx is the threshold determined by k-means [23]. ...8

Fig 1.8 The operation of multi-ART2 neural network [24]. ...8

Fig 2.1 Spindle cutting system model ... 11

Fig 2.2 regenerative effect [27] ...12

Fig 2.3 The geometric approximation of tooth path [30]...13

Fig 2.4 Discretization of time period. Each time delay is separated into single time step then computed by numerical integration. ...16

Fig 2.5 Floquet bifurcations: (a) Secondary Hopf bifurcation. (b) Period-doubling bifurcation. (c) Cyclic-fold bifurcation. [33] ...19

Fig 2.6 Stability lobe diagram (a) three-dimensional perspective. (b) font view...21

Fig 3.1 Daubechies 10 wavelet is one of mother wavelets [36]. ...23

Fig 3.2 Relationship between scale factor and time frequency resolution. ...24

Fig 3.3 Wavelet decomposition of time domain signal [22]. ...25 Fig 3.4 Wavelet packet tree with corresponding high-pass and low-pass filters. The

(16)

shaded nodes indicate the node not to be produced by DWT [38]. ...26

Fig 3.5 Acceleration during stability transition of cutting process. ...30

Fig 3.6 Schematic of SVM [41]. To distinguish blue points and red points, data dimension is converted from 2D to 3D space, and the black hyperplane separates these two classes. ...36

Fig 3.7 Illustration of OCSVM [43]. The data points on the hyperplane are support vectors. ...38

Fig 3.8 k-distance( )p d( , )p o for k 5 ,and N_k( )p 6 . d( , )p o is the distance between p and o...40

Fig 3.9 Reachability-distance for k 4...41

Fig 3.10 The conception of LOF. For k 4, local reachability density of point p is much lower than its neighbors a b c, , and d , and local outlier factor is much higher than 1. ...42

Fig 3.11 Stratified k-fold cross validation for k=3. Each subset has about 75% A class data and 25% B class data the same as the total dataset. ...45

Fig 3.12 Chatter suppression algorithm ...46

Fig 3.13 Control flow of system ...48

Fig 4.1 Deta type C tapping center ...49

Fig 4.2 Experimental setup of cutting force measurement ...52

Fig 4.3 Mean cutting force in x direction...52

Fig 4.4 Mean cutting force in y direction...53

Fig 4.5 Cutting force fitted result with a_p 0.2mm in x direction ...53

Fig 4.6 Cutting force fitted result with a_p 0.2mm in y direction ...54

(17)

Fig 4.8 Schematic diagram of impact test ...55

Fig 4.9 Frequency response function and coherence in x direction. ...56

Fig 4.10 Frequency response function and coherence in y direction. ...56

Fig 4.11 Frequency response function in x direction ...57

Fig 4.12 Frequency response function in y direction ...57

Fig 4.13 First mode and second mode stability lobe diagram ...59

Fig 4.14 SLD’s comparison of semi-discretization and ZOA...60

Fig 4.15 3D stability lobe diagram...61

Fig 4.16 Experimental setup of acceleration measurement. ...62

Fig 4.17 Cost of OCSVM. v was the parameter in equation (3.48) and  is the parameter in equation (3.37). ...64

Fig 4.18 Cost of LOF. k was the numbers of neighbors in equation (3.54) ...65

Fig 4.19 The actual distribution of entropy in the dataset. Entropy increases when chatter occurs. ...66

Fig 4.20 True answer of the dataset. ...67

Fig 4.21 Predicted answer of the dataset...67

Fig 4.22 User interface of intelligent chatter suppression module. ...68

Fig 4.23 Acceleration and RWPEE with 4700rpm initial spindle speed. The blue dotted line was the time when chatter was detected. (a) signal in x direction. (b) signal in y direction. (c) signal in z direction. ...69

Fig 4.24 Acceleration and RWPEE with 5800rpm initial spindle speed. The blue dotted line was the time when chatter was detected. (a) signal in x direction. (b) signal in y direction. (c) signal in z direction. ...70 Fig 4.25 Acceleration and RWPEE with 6000rpm initial spindle speed. The blue dotted line was the time when chatter was detected. (a) signal in x direction. (b)

(18)

signal in y direction. (c) signal in z direction. ...71 Fig 4.26 Acceleration and RWPEE with 6700rpm initial spindle speed. The blue dotted line was the time when chatter was detected. (a) signal in x direction. (b) signal in y direction. (c) signal in z direction. ...72 Fig 4.27 Acceleration and RWPEE with 6800rpm initial spindle speed. The blue dotted line was the time when chatter was detected. (a) signal in x direction. (b) signal in y direction. (c) signal in z direction. ...73

(19)

LIST OF TABLES

Table 3.1 Parameter’s representation in chatter detection ...43 Table 4.1 The parameters of FRF in x direction ...58 Table 4.2 The parameters of FRF in y direction ...58

(20)

Chapter 1 Introduction

1.1 Background and motivation

The Industry 4.0 Era has arrived. This era emphasizes that machine can integrate all aspects of intelligent technology components, including machine has customized service capabilities, machine, plants, and even companies can link through the cloud mutua lly, and machine has self-awareness, which contains self-monitoring and self-abnorma l diagnosis and self-repairing capabilities. Smart manufacturing is one of the branches of the Industry 4.0. It aims to improve processing yield, accuracy and reliability by collecting comprehensive data and exchanging information with each machine. In response to changes in the times, the demand for developing smart manufacturing has gradually emerged, especially in intelligent spindles.

According to Cao et al. [1], in Fig 1.1, intelligent spindles, or named intellige nt machine tools, should have three major capabilities: sensing, decision making and control.

First, whenever necessary, spindles should sense or monitor themselves in several aspects, including temperature, spindle balance, chatter, and etc. The more information of spindle, the more understanding of the spindle status. Second, if any abnormal condition is detected, spindles should make decision by themselves to troubleshoot. It could be imaged that the spindles have its own mind to make decision. Third, after making decision, they should control by themselves to achieve the goal. These three capabilities form a cycle that gives the spindle intelligence, so they are called intelligent spindles. This thesis will focus on applying these capabilities on chatter suppressions.

(21)

Fig 1.1 An integrated concept of intelligent spindles [1].

Chatter has been one of the significant issues in cutting process of spindles for the past decades. As a classic but tricky problem, Taylor [2] gave it a statement that “chatter is the most obscure and delicate of all problems facing the machinist” in 1907. This phenomenon is mainly caused by the self-excited vibration that the vibration on tool tip has a phase delay from the wave of workpiece surface left before. When chatter occurs during the cutting process, it will not only leave an uneven workpiece surface shown in Fig 1.2, but also make serious noises. Even worse, the tool will wear out and the structure of workpiece will fatigue or even be damaged. Chatter severely limits the productivities and precision pf manufacturing; as a result, this thesis proposes a systematic solution including chatter detection and chatter suppression.

(22)

Fig 1.2 Measured vibration signal and the corresponding workpiece surface [3].

1.2 Literature review

1.2.1 Stability lobe diagram

In the 1980s, Tlusty et al. first defined cutting process as a non-linear system [4] and visualized a stability lobe diagram (SLD) to classify stable cut and unstable cut at differe nt spindle speed and depth of cut [5]. It was not only the prediction of stability but also the basis of chatter suppression strategy. In 1995, Altintas et al. presented a new method that used Fourier series expression of time variant cutting force model to calculate analytica l solutions and redraw the SLD [6]. Although this method simplified the cutting force function and saved computing time, it was still not reliable to predict chatter based on the experimental results of Faassen et al. [7] shown in Fig 1.3.

(23)

Fig 1.3 Stability lobe diagram compared with experimentally performed cuts [7].

After that, scholars dedicated themselves to making SLD more accurate one after another. Solis et al. redefined the real and imaginary components of the transfer functio n to improve the vibration model [8]. Merdol and Altintas detailed a multi frequency solution [9]. However, all above methods cannot get rid of analytical method. Until 2002, Insperger et al. presented a new numerical method called semi-discretization to solve the spindle vibration model [10]. According to Muñoa et al. [11] shown in Fig 1.4, semi- discretization and multi frequency method demonstrated the same solution but were more accurate than zero order approximation reported in [6]. Therefore, this research will calculate spindle vibration model and SLD based on semi-discretization.

(24)

Fig 1.4 Stability lobe diagrams obtained from different methods [11].

1.2.2 Chatter detection

Since the concept of chatter was proposed, several techniques of chatter detection have been developed. In time domain, R. Du et al. detected chatter based on the probability distribution of cutting force signal [12]. Soliman and Ismail presented a strategy by monitoring R-value computed from spindle drive current [13]. In frequency domain, fast Fourier transform (FFT) has been extensively applied on chatter detection.

Toh used FFT to analyze optimum cutter path orientation to avoid chatter [14]. For example, Fig 1.5 showed the chatter frequency was dominant based in FFT.

Tangjitsitcharoen defined three feature ratios by power spectrum density of dynamic cutting force [15]. Han et al. chosen 16 sampling points to compute FFT as fast as possible [16]. Although time domain and frequency domain detection were usable, their sensitivity

(25)

in high frequency was insufficient to detect chatter earlier.

Fig 1.5 FFT analysis on force signatures in y direction [14].

To improve the sensitivity of detection, some scholars studied signal processing methods in time frequency domain, including wavelet transform [17], wavelet packet transform [18], [19], [20] and Hilbert–Huang transform [21]. Yen et al. first introduced the concept of wavelet into vibration monitoring in 2000 [22]. After that, more and more scholars used wavelet transform or wavelet packet transform as the feature to detect chatter. Yoon and Chin verified wavelet transform is more reliable than FFT on chatter detection [17]. For example, Fig 1.6 showed the chatter property was evident in level 3 and level 4 based on wavelet transform. Yao et al. simultaneously used wavelet transform and wavelet packet transform as the feature to monitor chatter [18]. Sun et al. defined a weighted wavelet packet entropy to detect chatter and verified its effect on earlier detection [19], [20]. On the other hand, Cao et al. presented Hilbert–Huang transform, which worked through Intrinsic Mode Functions and Empirical Mode Decomposition to acquire instantaneous frequency for chatter detection [21].

(26)

Fig 1.6 Behavior of the detail coefficients wavelet transformed in levels 1 to 9 for cutting force [17].

After having good features from signal processing, it is essential to defined the threshold to identify stable data and chatter data. In the past, setting threshold values often relied on experience and experiments, which was inefficient and uncertain. Therefore, classification and clustering algorithm were presented in threshold criteria [18], [23], [24], [25], [26]. Yao et al. used support vector machine, which classified data in hyperspace, to separate vibration signal into stable, transition and chatter data [18]. Tangjitsitcharoen et al. classified features based on k-means clustering, which would find representative data as the cluster centers to classify data [23]. Li, Lamraoui, HINO et al. respectively presented different types of neural networks to identify vibration signal [24], [25], [26].

However, most above algorithms could only detect fully developed chatter. It was too late that the workpiece surface or the tool might had been damaged. As a result, this thesis presents two anomaly detection algorithms to monitor cutting process by detecting abnormal signals before it has been fully developed for early detection.

(27)

Fig 1.7 Illustration of chatter detection based on k-means clustering. Cx is the threshold determined by k-means [23].

Fig 1.8 The operation of multi-ART2 neural network [24].

(28)

1.3 Thesis organization

This thesis establishes a complete chatter-suppressed module includ ing preprocessing, real-time detection, and automatic suppression. First of all, in Chapter 2, the spindle vibration model is established and is analyzed via semi-discretization. The influence of feed rate is derived in detail in order to close to the reality. This model is used to calculate chatter frequency and construct stability lobe diagram. In Chapter 3, the relative wavelet packet energy entropy is derived as the feature of chatter detection.

Wavelet packet feature has been proven to have better sensitivity in the high frequency band. After that, two anomaly detection algorithms－One-class support vector machine and local outlier factor are used to set the threshold between stable data and chatter data based on the feature. After the SLD and the chatter detection system are constructed, the strategy of chatter suppression will be decided. All above will be integrated into the intelligent chatter suppression module.

Chapter 4 starts with the experiments of spindle parameters, which are used as the known parameters of SLD. Furthermore, the model used to classify data is trained and validated here, too. After preprocessing finish, it can proceed the cutting test. This part verifies the reliability of the module in real-time chatter detection and suppression. Last, conclusion and future work are placed in Chapter 5.

(29)

Chapter 2 Spindle chatter theory

2.1 Spindle vibration model with variable time delay

Assume 2-DOF model of end-milling shown in Fig 2.1. System boundary is fixed on the bearing of spindle and only workpiece can be moved for cutting. All the system, including cutting tool, bearing, vise and etc., are assumed as rigid bodies except workpiece and spindle. With these assumptions, the equation of motion can be written as

( )t  ( )t  ( )t  ( )t

Mq Cq Kq F (2.1)

Where M , C , K are modal mass matrix, damping matrix, and stiffness matrix, respectively. q(t) is the displacement of the tool, and F( )t is the cutting force in x direction and y direction, respectively defined as

0 0

x y

m m

 

  

 

M (2.2)

0 0

x y

c c

 

  

 

C (2.3)

0 0

x y

k k

 

  

 

K (2.4)

( ) ( )

( ) t x t

y t

 

  

 

q (2.5)

( ) ( )

( )

x y

t F t F t

 

  

 

F (2.6)

(30)

Fig 2.1 Spindle cutting system model

By separating into normal and tangential components, the i^th tooth of cutting force can be represent as

, , ,

,y , ,

( ) (t)cos ( ) ( ) sin ( ) ( ) (t)sin ( ) ( ) cos ( )

i x i t i i n i

i i t i i n i

F t F t F t t

 

 

   

 (2.7)

In most cases [27], [28], [29], the cutting force model of the milling process are in exponential forms expressed as

,

( ) ( ) ( )

F

X

i t i t p i

X

i n i n p i

F t g t K a h t F t g t K a h t

 

 

 (2.8)

Where K_t and K_n are tangential and normal cutting force coefficients, respective ly,

ap is axial depth of cut, X_F is exponent terms, and g t_i( ) is the contact functio n shown below

1 , if ( ) ( ) ( ) ( ) 0

st i ex

i

t t t

g t    

  (2.9)

( )t F

(31)

Where _i( )t is the angular position of i^th tooth. _st( ) and t _ex( )t are enter angle and exit angle, which are time-variant and will be discussed later.

Next, h t_i( ) is chip thickness of i^th tooth related to chip load f and tooth path with perturbation [27]. Assuming model of chip thickness is the circular tooth path, it can be approximated as

( ) [ ( ) ( )]sin ( ) [ ( ) ( )]cos ( )

i i i

h t  f x t  x t  t  y t  y t  t (2.10) Where



is time delay of each tooth. In the cutting process, due to inconsistency of chip thickness, the new cutting surface will be out of phase with the old surface, which is called regenerative effect and is shown in Fig 2.2. The self-excited vibration generated by this effect indirectly leads to the occurrence of chatter.

Fig 2.2 regenerative effect [27]

In most of previous papers, time delay



is generally assumed to be equal to time period T . However, it’s effected by chip load because of tool vibration causing the variation of time delay [30]. According to the geometric approximation of tooth path shown in Fig 2.3, it is given as

( ) 120

(2 cos (t))

i

t R

N R fN

  

 

 

  (2.11)

Where R is radius of tool, N is the number of tooth, and  is spindle speed. In Fig

(32)

2.3, during a period of time delay, the center of cutting tool move from O to ₁ O₂, and the contact point of tooth move form B to A . Besides, the time delay related to chip load also change the contact angle between tool end and cutting surface of workpiece.

( ) ( )

120

i i

t f N

t R

  ^

  (2.12)

Therefore, it can determine the cutting angles for both up-milling and down-milling.

For up-milling

( ) ( ), ( ) cos (1 21 )

st t i t ex t aD

     ^  (2.13)

For down-milling

( ) cos (21 1), ( ) ( )

st t aD ex t i t

  ^       (2.14)

Where a_D is radius depth of cut ratio. _st and _ex are the angular position where the cutting teeth enter and exit the workpiece, respectively.

Fig 2.3 The geometric approximation of tooth path [30].

(33)

Substituting equation (2.8), (2.10) and (2.11) into (2.7), the cutting force can be written as

, 1

,y 1

( )[( ( ( )) ( )) sin ( ) ( ( ( )) ( )) cos ( )]

( )

( )[( ( ( )) ( )) sin ( ) ( ( ( )) ( )) cos ( )]

F

N

X

i x i i i i

i N

X

i i i i i

i

S t f x t t x t t y t t y t t

t

S t f x t t x t t y t t y t t

   



      

 

      





F (2.15)

Where

,

,y

( ) ( )(K cos ( ) sin ( ))

( ) ( ) ( )( K sin ( ) cos ( ))

i x p i t i n i

i

i p i t i n i

S t a g t t K t

t S t a g t t K t

 

 

    

S (2.16)

Using Taylor series expansion and neglecting high ordered terms, the cutting forces are simplified into the linear form

1

( ) ( )( sin ( )) {1 [( ( ( )) ( )) ( ( ( )) ( )) cot ( )]}

[ ( )( sin ( )) ( )( ( ( )) ( ))]

F

N

X

i i i F i i i

i N

X

i i i i

i

t t f t X f x t t x t y t t y t t

t f t t t t t

   

 





      

   



F S

S G q q

(2.17)

Where

1 1 1

, ,

1 1 1

,y ,y

( ) sin ( ) ( ) sin ( ) cos ( )

( ) ( ) sin ( ) ( ) sin ( ) cos ( )

F F F F

X X X X

i x F i i x F i i

i X X X X

i F i i F i i

S t X f t S t X f t t

t S t X f t S t X f t t

  

  

 

  

 

 

G (2.18)

It can be seen from the regenerative effect that the state of the cutting surface at cutting process is inconsistent. Therefore, it is known that the spindle cutting system model is a delay-differential equation, which is abbreviated as DDE. The solution of DDE is periodic and the initial conditions depend on solution of previous period. Then, the equation (2.1) can be transformed into the state-space form

(34)

1 1

1

1 -1

1

0

( ) ( ) ( ) ( ) ( ( ))

( ) sin ( )

where

0 0

( ) ( ) , ( )

[ ( )]

( )

0 0

( ) ( ) 0

F F

N

i i X X

i i i

i

N i i

i

t t t t t t

t f t

t t t

t t

 ^ 









 

 

     

 

 

 

   

     

 

  

 



Q A Q B Q

M S

Q q A

M K G M C

q

B M G

(2.19)

Next section will introduce the method to solve DDE.

2.2 Semi-discretization

In practice, because it’s a nonlinear and discontinuous system, it is difficult to find an analytical solution for DDE. Instead, the approximate solutions computed by semi- discretization method is widely used [30], [31], [32].

Assume the system is composed of periodic orbit Q₀( )t and perturbation X( )t . The system can be written as

( )t  0( )t  ( )t

Q Q X (2.20)

There is only the periodic orbit in ideal case of stable milling when no self-excited vibrations occur. Substituting equation (2.20) into (2.19), the perturbation is written as

1

( ) ( ) ( ) ( ) ( ( ))

N

i i

i

t t t t t  t



 





X A X B X (2.21)

In order to discretize the time period T , let T be distributed into k intervals.

t T

  k (2.22)

Therefore, in single time step, the coefficient matrices from (2.19) are approximated as 1

(35)

1

,

( ) 1 ^j ( )

j

t

j i t t i t dt

t

 





B B (2.24)

and also, time delay becomes

,

120

[2 cos ( 0.5 )]

j i

i j

R

N R fN t t

 

 

     (2.25)

Fig 2.4 Discretization of time period. Each time delay is separated into single time step then computed by numerical integration.

The discrete relationship between time delay and time step is expressed in Fig 2.4 and is given as

, ( , , 0.5)

j i lj i lrj i t

     (2.26)

, ,

mod( ^{j i} 0.5 )

j i

lr t

t

  

  (2.27)

(36)

,

, ^{j i} , 0.5

j i j i

l lr

t

  

 (2.28)

Where lr_{j i}_, is residual ratio of the time delay and mod is the function that returns the remainder after a number is divided by a divisor. Depending on the initial conditions, the time delay



may greater than or less than time period T .

According to the lr_{j i}_, , the delay terms are approximated by distributing the weights of time step

, ,

, ,i

, , 1

( ) ( 0.5 ) ( )

(1 ) (t ) ( )

j i j i

j i j j j l lr

j i j j j i j l

t t t t

lr lr t

  _ _

  

     

  

X X X

X X (2.29)

After substituting equation (2.29) and the initial value X( )t_j X_j into equation (2.21) and solving the equation with the assumption that A_j( )t is invertible for all j , the analytical solution over the interval of single time step is written as

, , , ,

1 , , 1 1

1

( (t ) (t ))

j i j i j i j i

N

j j j j l j l j l j l

i

    



 





X P X R X R X (2.30)

Where the coefficient matrices are

j t j e^A^

P (2.31)

,

1

, ( ^j ) ,(1 ,)

j i

t

j l  e^A^  ^j j i lrj i

R I A B (2.32)

,

1

, 1 ( ^j ) , ,

j i

t

j l _  e^A^  ^j j ilrj i

R I A B (2.33)

Then, equation (2.30) can be reconstructed into the discrete map

1

j  j j

z D z (2.34)

where

1 max

max max

[ , , , ]

mod( ) 1

T T T T

j j j j l

l  

 



  

z X X X

(2.35)

(37)

and

, ,i

, , 1

1

0 0

0 0 0 0

0 0 0 0 0 0 0 0

j i j

j l j l

j

N j

i





 

 

   

 

   



R R

P I D I

I I

(2.36)

From the discrete map, it can be mapped as

0 k 

z Φz (2.37)

Where Φ_st is the transition matrix.

1 1 0

st  k

Φ D D D (2.38)

2.3 Chatter frequency and stability lobe diagram

According to Floquet theory [33], the eigenvalues of transition matrix , which is called Floquet multipliers, can resolve the stability of DDE (2.19). If 1, the periodic solution of system is asymptotically stable. On contrast, if 1, the system is unstable.

By distinguishing the value of Floquet multipliers from the unit circle on complex coordinate, three states, (a)   1 Im( ) 0, (b)  1 and (c) 1 are made, which is shown in Fig 2.5. The state that crosses unit circle in different properties is called bifurcation. Discuss the above three states respectively:

(38)

Fig 2.5 Floquet bifurcations: (a) Secondary Hopf bifurcation. (b) Period-doubling bifurcation. (c) Cyclic-fold bifurcation. [33]

(1) Secondary Hopf bifurcation

When   1 Im( ) 0 , Floquet multipliers crosses the unit circle in the form of complex conjugate, which is shown in Fig 2.5(a) and is called secondary Hopf bifurcation. It’s the most common bifurcation in cutting process. The periodic solution is asymptotically stable before the occurrence of bifurcation. After bifurcation happens, a new frequency f_c, which is called basic chatter frequency, will be superposed on the origin frequency. f_c is not related to tooth path frequency and is only affected by .

Im(ln )

c 2

f T



  (2.39)

Then the secondary Hopf bifurcation frequency consists of tooth path frequency f_tpe and basic chatter frequency f_c.

H c tpe,

f   f kf kZ (2.40) (2) Period-doubling bifurcation

When  1, period-doubling bifurcation, shown in Fig 2.5(b) is formed. A new bifurcation frequency f_pd is generated and is half of the tooth pass frequency.

(39)

2

tpe pd

f  f (2.41)

So the period-doubling bifurcation frequency is

PD pd tpe,

f  f kf kZ (2.42)

(3) Cyclic-fold bifurcation

When  1, cyclic-fold bifurcation, shown in Fig 2.5(c) is formed. Cyclic- fold bifurcation frequency is on tooth pass frequency and its harmonics.

T tpe,

f kf kZ (2.43)

These bifurcation frequencies superpose on natural frequency of spindle system which result in abnormal vibration in cutting process [34]. This abnormal frequency is dominant chatter frequency. Generally, only first mode of natural frequency will occur.

Based on this conception, which frequency should be monitored is determined.

Furthermore, choosing spindle speed, depth of cut, and chip load as control parameters, the SLD can be illustrated by computing the Floquet multipliers of transitio n matrix. Taking parameters in reference [31] as an example, with a_D 0.05 , N 2 ,

1

107 10 N/m6 ^X^F

Kt   ^ , K_n 40 10 N/m ⁶ ¹^^X^F , X_F 0.75 and f 0.1 mm/tooth , the SLD is shown in Fig 2.6. the area below the critical surface is stable. In contrast, the area above the surface is unstable.

After calculating the above-mentioned bifurcations with the measured spindle parameters, the stability of the spindle system can be determined. The next chapter will establish the spindle chatter threshold by the chatter frequency, which serves as a basis for monitoring the cutting state of the spindle.

(40)

(41)

Chapter 3 Real-time detection and suppression

3.1 Monitoring strategy

3.1.1 Wavelet transform

According to spindle vibration model derived in chapter 2, it can be seen that chatter frequency of spindle is related to natural frequency. FFT is one of the popular frequency analysis for chatter detection [14], [15], [16]. By measuring acceleration signal during cutting process and using FFT, occurrence of chatter will be detected. However, when chatter is detected by FFT, generally the chatter frequency has become to dominant frequency, and also the workpiece surface has been damaged. Therefore, the sensitivity of chatter frequency must be increased in order to truly achieve real-time chatter detection.

Wavelet transform has been proposed by Grossmann and Morlet et al. since 1984 [35]. In recent years, many scholars have applied it to chatter detection for two reasons.

According to Heisenberg's uncertainty principle, it is impossible to have both high resolution in time domain and frequency domain with any frequency analysis. For instance, for high resolution in frequency domain, it’s necessary to extend time to acquire data for more samples. Since chatter frequency and tooth path frequency were determined, it only had to monitor the frequency band which chatter frequency was located to determine whether chatter occurs. In other words, it is feasible to improve time resolutio n by sacrificing frequency resolution. Second, unlike FFT, wavelet transform uses a different base called mother wavelet, which has higher sensitivity in high frequencies than FFT with a sinusoidal-function base. One of the mother wavelets is shown in Fig 3.1. As a consequence, it is possible to detect chatter earlier by wavelet transform.

(42)

Fig 3.1 Daubechies 10 wavelet is one of mother wavelets [36].

The definition of mother wavelet is shown below.

,

( ) 1 ( )

a b

t b

t a a

   ^ (3.1)

Where a b, R are scale and shift factors, respectively. When a decreases, mother wavelet is compressed. It has a small support in time axis, and frequency spectrum move toward high frequencies. Conversely, when a increases, mother wavelet becomes wider. Time axis is extended, and frequency spectrum move toward low frequencies. The relationship, as shown in Fig 3.2, can be used as the basis for real-time detection by utilizing the characteristics of high time resolution at high frequencies.

(43)

Fig 3.2 Relationship between scale factor and time frequency resolution.

In practical application, for finite-length discrete-time signal, the scale and shift factors can be discretized in binary a_j 2^^j, b_{j k}_, 2^^jk, where j k, Z. In Hilbert space ( )t L²( ) , the mother wavelet can be represented as

/ 2 ,

1, 2, 3,..., ( ) 2 (2 ),

1, 2, 3,...,

j j

j k

j J

t t k

k K

    ^ ^{   } ^ (3.2)

Where J is maximum level of decomposition and K is length of filter. Thus, signal

signal( )

f t can be expressed as

, ,

,

(t) ( ), ( ) ( )

signal signal j k j k

j k

f 



f t  t  t dt ^(3.3)

Where , is inner product. This expression, which is called discrete wavelet transform (DWT), shows any signal can represent as mother wavelet form by wavelet coefficie nts

( ), , ( )

signal j k

f t  t in wavelet series, and also implies wavelet coefficients contain all information from original signal f_signal(t). However, since DWT was an iteration process with heavy computing work, Mallat proposed the fast algorithm to compute much

(44)

efficiently [37]. This algorithm ignores the analysis of scaling functio n and mother wavelet function but constructs in hierarchical structure, which is called multiresolutio n analysis (MRA). According to MRA, the scaling function ( )t is defined with expressions shown below.

( ) 2 ( ) ( )

2 _k

t h k t k

 



  ^(3.4)

( ) 2 ( ) ( )

2 _k

t g k t k

 



  ^(3.5)

Where h k( ) and g k( ) are low-pass and high-pass filters, respectively. It can be seen that signal would be decomposed into high and low frequency bands, which limits the number of iteration with these expressions, instead transforming signal by iterating wavelet with all possible scale factors. As a result, time domain signal would be represented as J levels decomposed signal, which is shown in Fig 3.3.

Fig 3.3 Wavelet decomposition of time domain signal [22].

(45)

3.1.2 Wavelet packet transform

Despite wavelet transform splits data into different time-frequency resolution, it still retains a larger frequency band in the high frequency than in the low frequency, which could possibly influence chatter resolution. Therefore, wavelet packet transform has been developed. Wavelet packet transform, shown in Fig 3.4, decomposes data in both low- pass coefficient and high-pass coefficient respectively by generalizing MRA. Compared to wavelet transform, it provides more decomposition in the high frequency to obtain more detailed information. Therefore, wavelet packet transform (WPT) is more suitable to monitor chatter frequency in high frequency band.

Fig 3.4 Wavelet packet tree with corresponding high-pass and low-pass filters. The shaded nodes indicate the node not to be produced by DWT [38].

(46)

Wavelet packet function is written as

/2

, ( ) 2 (2 )

n j n j

Wj k t  W t k (3.6)

Where n is the parameter.

2 ^j 1

n ^  (3.7)

And the first two wavelet packet functions are:

0

0,0( ) ( )

W t  t (3.8)

1

0,0( ) ( )

W t  t (3.9)

For n2,3, 4,..., wavelet packet would be written with a recursive relationship.

2

1( ) 2 ( ) (2 )

n n

j j

k

W_ t 



h k W tk ^(3.10)

2 1

1 ( ) 2 ( ) (2 )

n n

j j

k

W_^ t 



g k W tk ^(3.11)

Similar to DWT, wavelet packet coefficient is represented as

, , ( ), ⁿ, ( )

j n k j k

c  f t W t (3.12)

Wavelet packet coefficient contains all information in the specific frequency band.

Because of this property, it will be applied in next section as a feature to detect chatter.

3.1.3 Relative wavelet packet energy entropy

In 1948, C. E. Shannon proposed a new concept of entropy in communication field [39]. In this new concept, according to R. Y. Liu et al. [40], entropy is defined as a measurement of uncertainty in a probability distribution, which would be used to estimate the information of random signals, and it is widely used in fault diagnosis. Hence, this research will utilize Shannon entropy as a feature to detect chatter.

Similar to FFT, any signal can be expressed in terms of mother wavelet via the

OCSVM與LOF演算法用於智慧主軸之即時監測與顫振分析抑制

國立臺灣大學工學院機械工程學研究所 碩士論文

Graduate Institute of Mechanical Engineering College of Engineering

National Taiwan University Master Thesis

OCSVM 與 LOF 演算法用於智慧主軸之即時監測與顫振分析 抑制

Real-time Chatter Detection, Analysis and Suppression for Intelligent Spindles Based on One-class Support Vector Machine

and Local Outlier Factor

姚泳辰 Yung-Chen Yao

指導教授：施文彬 博士、劉建豪 博士 Advisor: Wen-Pin Shih, Ph.D.

Chien-Hao Liu, Ph.D.

中華民國 107 年 6 月

June, 2018

誌謝

中文摘要

ABSTRACT

SYMBOL TABLE





CONTENTS

LIST OF FIGURES

LIST OF TABLES

Chapter 1 Introduction

1.1 Background and motivation

1.2 Literature review

1.2.1 Stability lobe diagram

1.2.2 Chatter detection

1.3 Thesis organization

Chapter 2 Spindle chatter theory

2.1 Spindle vibration model with variable time delay













 



2.2 Semi-discretization











2.3 Chatter frequency and stability lobe diagram

Chapter 3 Real-time detection and suppression

3.1 Monitoring strategy

3.1.1 Wavelet transform







3.1.2 Wavelet packet transform





3.1.3 Relative wavelet packet energy entropy

國立臺灣大學工學院機械工程學研究所碩士論文

OCSVM 與 LOF 演算法用於智慧主軸之即時監測與顫振分析抑制

指導教授：施文彬博士、劉建豪博士 Advisor: Wen-Pin Shih, Ph.D.