李忠群，劉 浪，段林升，肖檢冬，張偉峰

（湖南工業(yè)大學 機械工程學院，湖南 株洲 412007）

1 研究背景

作為一種高效的切削加工方式，銑削被廣泛應用于航空、航天及汽車制造等領域中。顫振是切削過程中的自激振動現(xiàn)象，它會制約銑削加工質量與加工效率，如在工件表面留下振紋，降低零件表面質量，嚴重時可能會損壞刀具和工件，甚至機床主軸。研究顫振的首要目的是根據(jù)穩(wěn)定性葉瓣圖選取切削參數(shù)，進而實現(xiàn)無顫振高效切削。

目前，顫振穩(wěn)定性分析方法可歸納為頻域分析和時域分析兩大類。頻域分析方法由于在分析過程中進行了近似處理，故處理速度較快，但是處理精度有限。時域分析方法由于綜合考慮了非線性切削力、刀齒跳離切削區(qū)以及過程阻尼等非線性因素的影響，因而分析精度較高[1]。時域分析方法主要有數(shù)值求解法、半離散和全離散求解法等。其中數(shù)值求解法應用較廣，它通過先求解銑削過程的動力學微分方程，獲得切削、刀尖位移等時間歷程參數(shù)，再施用某種穩(wěn)定性判據(jù)，進而確定當前切削條件的穩(wěn)定性。數(shù)值求解法直觀、簡單，但穩(wěn)定性判據(jù)的有效性則成為該方法是否可行的關鍵。

數(shù)值求解法常用的穩(wěn)定判據(jù)主要有FFT（fast Fourier transform）法、動態(tài)力比靜態(tài)力法[2]、刀尖位移統(tǒng)計法[3]以及峰值力法等，其中FFT法的運用較為廣泛。FFT法適合處理周期信號或平穩(wěn)信號，而不適合處理切削振動信號[4]。刀尖位移統(tǒng)計法將樣本方差作為顫振判定標準，但是在小徑向切深情況下預測精度欠佳。相較靜態(tài)力法，動態(tài)力只能對某一工況顫振穩(wěn)定性進行判斷，適應性較差。吳石等[5]以最大Lyapunov指數(shù)為判據(jù)，分別通過等高線法、輪廓法確定銑削顫振穩(wěn)定域。Sun Y.W.等[6]則將特征值是否大于1作為穩(wěn)定性判據(jù)。李忠群等[7]提出了一種用于預測銑削穩(wěn)定性的基于三階龍格-庫塔法的半解析方法，利用Floquet理論判斷特定切削狀態(tài)下的穩(wěn)定性。但是上述顫振穩(wěn)定性時域判據(jù)，都無法滿足在所有工況下都能準確識別顫振的要求，因此需進行改進。

鑒于顫振發(fā)生的典型特征是能量在頻率軸分布的改變，因此越來越來多的學者運用基于能量的方法在線判斷顫振現(xiàn)象。Rényi熵可準確反映切削過程中的能量變化，故本文先采用變步長龍格-庫塔法求解銑削動力學微分方程，然后采用Rényi熵分析其時間歷程參數(shù)，獲得銑削過程穩(wěn)定性葉瓣圖，以期為實現(xiàn)無顫振高效切削工藝提供參考。

2 變步長龍格-庫塔法求解銑削動力學微分方程

常用歐拉法或龍格-庫塔法求解銑削過程動力學微分方程。歐拉法將切線端點作為下一步起點，其計算效率較高但精度較低。而龍格-庫塔法因將求解區(qū)間細分，并對斜率進行加權平均，作為導數(shù)近似，故其在求解精度上有較大改善。本文采用變步長龍格-庫塔法求解銑削動力學微分方程，以確保計算精度，降低截斷誤差，減少計算量。

機床刀具系統(tǒng)的動力學微分方程可表示如下：

式中：X、Y分別為進給方向及其法方向；

mX、mY分別為X、Y方向系統(tǒng)模態(tài)質量；

cX、cY分別為X、Y方向系統(tǒng)阻尼；

kX、kY分別為X、Y方向系統(tǒng)模態(tài)剛度；

FX(t)、FY(t)分別為X、Y方向上刀具所受切削力。

以X方向為例，其微分方程可轉化為

其初始條件為

以變步長龍格-庫塔法求解銑削過程動力學微分方程描述[8]如下：

針對式（2），在[tn,tn+1]時間內，以步長h，運用龍格-庫塔法求振幅x的近似值，記為，此時截斷誤差為O(h5)，因此：

將步長h折半，[tn,tn+1]時間內，分兩次求得近似值，記為，每次截斷誤差為c(h/2)5，因此：

對比式（5）（6），可得兩者的誤差關系如下：

由式（7）可知，若求解區(qū)間步長縮短一半，則截斷誤差變?yōu)樵瓉淼?/16。故采用變步長龍格-庫塔法求得的振幅近似值為

而將步長折半后兩次計算的偏差為

對于給定精度ε，步長選擇也不同。當ε＞Δ時，將步長反復折半并計算近似值，直至誤差小于給定精度，此時的近似值即為最終結果；當ε＜Δ時，將步長反復加倍并計算近似值，直至誤差大于給定精度，再將此步長值進行折半處理并計算，所得近似值即為最終的結果。

同理，可以得到Y方向的最終計算結果。

通過上述方法，獲得給定切削條件下的銑削時域信號，通過施用某種時域穩(wěn)定性判據(jù)，就能判定當前的切削條件是否穩(wěn)定。按一定的增量改變主軸轉速和軸向切深，重復上述過程，就能得到穩(wěn)定性葉瓣圖，其仿真流程如圖1所示，其中動態(tài)切削厚度與瞬時切削力的詳細計算見參考文獻[9]。

圖1 銑削穩(wěn)定性數(shù)值求解流程圖Fig.1 Flow chart of numerical solution of chatter stability for milling process

為了驗證上文所述方法求解所得結果的準確性，對不同條件下的槽銑進行切削仿真分析。仿真所用刀具為φ12的硬質合金整體銑刀，刀齒數(shù)為2，螺旋角β=30°，鋁合金AL7075-T6工件材料，切削力辨識系數(shù)如下：切向切削力系數(shù)Ktc=796.0 N/mm2，徑向切削力系數(shù)Krc=168.0 N/mm2，軸向切削力系數(shù)Kae=222.0 N/mm2，切向刃口力系數(shù)Kte=27.7 N/mm，徑向刃口力系數(shù)Kre=30.8 N/mm，軸向刃口力系數(shù)Kae=1.5 N/mm。該條件下得到的仿真結果如圖2、3所示。

圖2 轉速為8 000 r/min、軸向切深為0.2 mm、每齒進給量為1.5 mm時的仿真結果Fig.2 Simulation results under specific cutting conditions（n=8 000 r/min，ap= 0.2 mm，fZ=1.5 mm）

銑銑削過程中，穩(wěn)定切削時，切削信號幅值波動收斂；而顫振時，切削信號幅值會大幅增長。由圖2所示穩(wěn)定切削時切削力和刀具振幅仿真結果可以得知，穩(wěn)定切削下X、Y、Z方向的切削力信號波峰分別為 100, 250, 70 N 左右，波谷約為-150, -10 , -5 N。X、Y方向刀具振幅波峰分別約為 0.008, 0.021 mm，波谷約為-0.002, -0.014 mm ，而且無論是切削力還是刀具振幅，在整個切削過程中，波峰波谷值變化不大。由圖 3 所示顫振時切削力和刀具振幅仿真結果，可以得知，0.018 s 前，切削力和刀具振幅信號的波峰和波谷均保持在一個穩(wěn)定值。之后，刀具振幅和切削力值大幅增長，并經過一段時間后穩(wěn)定，X、Y、Z方向切削力信號波峰分別約為 2 000, 4 000, 1 000 N，波谷約為 -2 800, -50, -10 N 。X、Y方向刀具振幅波峰分別約為 0.3, 0.5 mm，波谷約為-0.38, -0.40 mm。圖 2 和 3 所示仿真結果很好地反映了銑削過程的特點。因此，圖2 和 3 可證明本文所提求解銑削過程微分方程的方法是正確。

圖3 轉速為10 000 r/min，軸向切深為0.9 mm，每齒進給量為1.5 mm時的仿真結果Fig.3 Simulation results under specific cutting conditions（n=10 000 r/min，ap= 0.9 mm，fZ=1.5 mm）

3 基于Rényi熵的時域穩(wěn)定性判據(jù)

3.1 Rényi熵定義

熵在數(shù)學和信息學中分別用來表示問題的不確定性和系統(tǒng)的復雜性。問題的不確定性和復雜性越大，熵值越大，反之越小。作為熵推廣的Renyi熵，可用于判斷信號信息量和復雜度，被廣泛用于圖像配準[10-11]、無線電頻譜感知等信號處理領域。Rényi熵是一個無量綱指標，當概率集的所有值幾乎相等時，Rényi熵值較大，對應系統(tǒng)復雜度較大；若只有少數(shù)值較大，其他值較小，則Rényi熵值較小，對應系統(tǒng)復雜度較低。

時域信號x(t)的Rényi熵計算過程如下：

1）計算時域信號x(t)的幅度譜，為

式中：N為數(shù)據(jù)長度；

X(ω)為x(t)的頻譜。

2）通過歸一化處理，可估計出頻譜的概率密度函數(shù)，為

式中：s(fi)為頻率分量fi的頻譜能量；

pi為相應的概率密度。

3）計算相應的Rényi熵H，定義為

為了比較不同的工作條件，用系數(shù)log2N進行歸一化處理，可得歸一化處理后的Rényi熵E為

3.2 閾值確定

Rényi熵E經歸一化處理后是(0，1)范圍內的無量綱指標，其中，1表示不確定性最大，0表示不確定性最小。

為揭示Rényi熵是否可檢測顫振，并進一步確定顫振閾值，根據(jù)顫振穩(wěn)定性葉瓣圖分別選取穩(wěn)定狀態(tài)和非穩(wěn)定狀態(tài)所對應的切削參數(shù)進行銑削過程時域仿真，并計算其Rényi熵。仿真參數(shù)設置如下：銑刀直徑D=12 mm，刀齒數(shù)N=2，螺旋角β= 30°，Ktc=796.0 N/mm2，Krc=168.0 N/mm2，Kte= 27.7 N/mm，Kre= 30.8 N/mm，工藝系統(tǒng)模態(tài)參數(shù)見表1。穩(wěn)定與非穩(wěn)定狀態(tài)下刀尖振動位移及Rényi熵見圖4、5。

表1 工藝系統(tǒng)模態(tài)參數(shù)Table 1 Modal parameters of the process system

圖4 穩(wěn)定切削條件下的振幅及其Rényi熵Fig.4 Amplitude with its Rényi entropy under a stable cutting condition

圖4顯示，穩(wěn)定切削時，信號的Rényi熵除開始時存在少許波動外，都一直穩(wěn)定在0.9左右，能量主要消耗在切削頻率及其諧振頻率上，確定性較大。圖5顯示，在非穩(wěn)定即顫振切削時，Rényi熵值開始時約為0.9，但隨著切削進行，進入顫振狀態(tài)，能量不僅消耗在切削頻率及其諧振頻率上，并且有部分能量消耗在顫振頻率上，其Rényi熵值逐漸減少，直至在0.8左右上下波動，即銑削顫振發(fā)生時，Rényi熵的值為0.8左右。

圖5 顫振條件下的振幅及其Rényi熵Fig.5 Amplitude with its Rényi entropy under a cutting condition

為進一步揭示顫振發(fā)生時的閾值，改變銑削的軸向切深，計算不同軸向切深下對應的Rényi熵，并且計算不同切削狀態(tài)下從穩(wěn)定到顫振的Rényi熵，所得結果如表2所示。

表2 切削狀態(tài)與Rényi熵的關系表Table 2 Relationship between cutting state and Rényi entropy

表2顯示，穩(wěn)定切削時的Rényi熵在0.9左右；當顫振發(fā)生時，Rényi熵在0.8左右。因此，本文將發(fā)生顫振的Rényi熵閾值設為0.83，低于該值，則判定為有顫振現(xiàn)象發(fā)生。

4 仿真與試驗驗證

為驗證本文提出的Rényi熵穩(wěn)定性閾值及基于Rényi熵的銑削過程穩(wěn)定性預測方法，在銑削穩(wěn)定性時域仿真的基礎上進行穩(wěn)定性驗證試驗。驗證試驗在五軸加工中心Mikron HSM/600U上進行，所用刀具為直徑12 mm、螺旋角30°、2齒硬質合金整體銑刀，工件材料為AL7075-T6，由切削力系數(shù)辨識試驗得到切削力系數(shù)同前一節(jié)，通過錘擊試驗獲得的工藝系統(tǒng)模態(tài)參數(shù)如表1所示。

在進行顫振驗證試驗時，首先固定主軸轉速，將初始軸向切深設定為ap=0.07 mm，然后以0.03 mm的步長逐漸增加，直至監(jiān)測到切削顫振，然后改變主軸轉速重復上述過程。通過監(jiān)控切削再生及觀察已加工表面是否產生振紋來確定加工過程中是否產生顫振。基于Rényi熵的銑削穩(wěn)定性葉瓣圖預測結果與切削實驗結果對比如圖6所示。圖中，細實線為基于Rényi熵方法仿真得到的穩(wěn)定區(qū)與顫振區(qū)分界線，圖線下方為穩(wěn)定切削區(qū)，上方為顫振區(qū)；方塊點為試驗穩(wěn)定切削點，圓圈點為試驗顫振切削點。該圖顯示，銑削過程顫振穩(wěn)定性仿真結果與實驗結果較為一致，從而驗證了本文提出的基于Rényi熵的銑削過程顫振穩(wěn)定性預測方法及用于判定顫振的Rényi熵閾值是正確的。

圖6 基于Rényi熵的銑削穩(wěn)定性預測結果與實測結果對比Fig.6 Comparison of milling stability prediction results based on Rényi entropy with measured results

5 結論

1）考慮再生作用的銑削動力學微分方程，運用變步長龍格-庫塔法求解，可獲得諸如切削力、刀具振動等較為準確的銑削過程時域信號；

2）Rényi熵可以被作為一種判斷銑削過程穩(wěn)定性的有效判據(jù)；

3）作為銑削過程穩(wěn)定性判據(jù)的Rényi熵，其閾值為0.83。