一種基于多步回溯算法的銑削穩(wěn)定性預測方法

于福航， 李茂月， 嚴復鋼， 劉獻禮， 梁越昇

(1. 哈爾濱理工大學 機械動力工程學院， 哈爾濱 150080；2. 佐治亞理工學院 喬治W.伍德拉夫機械工程學院， 亞特蘭大 30318)

薄壁類零件由于其輕便、高效的結構特性被廣泛地應用于航空航天、汽車制造和國防領域[1]。然而，薄壁類零件由于其低剛度，弱阻尼等特點往往不利于銑削過程中的穩(wěn)定切削，甚至引發(fā)顫振，最后影響加工工件的表面質量和加工效率等[2]。在薄壁類零件銑削過程中，顫振產生的最主要原因之一就是再生效應[3]。因此，考慮再生效應的薄壁工件銑削顫振預測仍是目前學術界和工業(yè)研究的熱點問題[4]。

銑削穩(wěn)定性預測在選擇合適的加工參數以實現無顫振加工中起到了重要的作用。最常用的穩(wěn)定性預測方法就是通過預測穩(wěn)定性葉瓣圖以顯示切削深度與主軸轉速之間的臨界關系，從而在穩(wěn)定切削區(qū)域選擇切削參數[5]。早在20世紀，由Tobias等[6]最早提出再生顫振，隨后又由Tlusty等[7-8]闡明了自激振動現象。此后，各國學者據此針對銑削穩(wěn)定性開展了大量的研究工作，并提出了一系列地穩(wěn)定性預測方法。這些方法主要聚焦于如何更加快速和準確地構建穩(wěn)定性葉瓣圖，可分為頻域解析法和時域法。

頻域解析方法以Altintas等[9]提出的零階頻域法(單頻率法，ZOA)為代表，該方法通過將周期銑削力系數與近似系統(tǒng)進行變換，再根據穩(wěn)定性判據計算軸向切深臨界值，隨后Ozturk等[10]將該方法推廣到五軸球頭刀銑削加工穩(wěn)定性分析，Altintas團隊也通過這種方法分析了不等齒距[11]、變螺旋角[12]等因素對銑削穩(wěn)定性的影響，但由于忽略了傅里葉級數中的高次項，這種方法的計算精度無法滿足小徑向切深時的要求[13]。針對這一問題，Merdol等[14]基于已完成的工作，提出了一種多頻率法，該方法被開發(fā)用于低浸沒式銑削，可以同時滿足小徑向切深和大徑向切深的要求，但由于使用了高階項進行近似，因此計算時間被大大地加長了。

時域法最初由Insperger等[15-16]提出，即半離散法(Semi-Discretization Method,SDM)，該算法是通過對銑削時滯微分方程中延遲項進行離散化處理后進行計算，隨后又提出了一種高階的半離散法[17]來近似周期性時滯系統(tǒng)。Ding等[18]通過提出一種全離散法(Full-Discretization Method，FDM)來預測銑削穩(wěn)定性，包括一階全離散法(1stFDM)和二階全離散法(2ndFDM)，這種方法是通過線性插值對狀態(tài)項和延遲項進行離散化處理，進而求出近似的數值解的方法。Guo等[19]通過將全離散法中的Euler插值變成Newton插值，提出了一種三階全離散法，并具有比一階、二階全離散法更高的計算精度和更快的計算效率。Li等[20]通過對微分方程的狀態(tài)項、延遲項、微分項全部離散化，提出了一種完全離散法(CDSEM)。Li等[21]基于完全離散法，提出了一種完全離散化的龍格庫塔方法(RKCDM)來預測銑削過程中的穩(wěn)定性問題。Jiang等[22]基于改進的精確時間積分算法，開發(fā)了一種高效的二階半離散法(2ndSDM)，大大提高了原半離散算法的計算效率。

在隨后的研究中，為了在滿足一定預測精度的同時提高預測效率，各國學者對銑削穩(wěn)定性的預測方法逐漸由線性單步法轉變?yōu)榫€性多步法。Eksioglu等[23]最早提出將Simpson公式應用于銑削穩(wěn)定性的求解過程中，隨后Zhang等[24]詳細地比較了基于Simpson公式的銑削穩(wěn)定性預測方法(SBM)與其他算法的計算精度和計算效率。Zhang等[25]基于泰勒公式，采用差分法和外推法的思想，提出一種數值微分的銑削穩(wěn)定性預測方法。智紅英等[26]根據隱式Adams公式提出了一種銑削穩(wěn)定性預測方法，該方法較1stFDM和2ndFDM方法具有收斂速度快的優(yōu)點，隨后智紅英等[27]又提出了一種基于Hamming公式的銑削穩(wěn)定性預測算法，并得到了高于1stFDM和2ndFDM的計算精度。Dai等[28]提出一種顯式精細積分法(PIM)來預測銑削過程的穩(wěn)定性，該方法可以顯著降低由差分格式引起離散化誤差，從而提高計算精度。

在考慮過程阻尼銑削穩(wěn)定性預測算法的研究中，Tlusty等[29]最先指出了隨著切削速度降低，穩(wěn)定性邊界有所提高，并將其歸結為過程阻尼效應的影響。隨后，Wu[30]指出過程阻尼主要來源于刀具后刀面壓入工件的材料體積，并提出了一種基于壓入體積和過程阻尼系數的過程阻尼力模型。Budak等[31]將過程阻尼當作附加阻尼納入動力學模型中，隨后采用簡單實驗法獲取穩(wěn)定性極限。近兩年，Wan等[32]提出了一種通過分離實測銑削力中的靜態(tài)過程阻尼力直接辨識過程阻尼系數的通用方法。Li等[33]提出了一種考慮過程阻尼效應的二階半離散法用于預測銑削過程的穩(wěn)定性。

然而，目前考慮過程阻尼的銑削穩(wěn)定性預測算法仍然較少，且對于附加過程阻尼力的求解比較耗時。針對這一問題，本文提出考慮過程阻尼效應的一種多步回溯算法，能夠有效地提高求解的效率和精度。

1 銑削動力學模型

切削過程顫振發(fā)生的原因可分為再生型振動、摩擦型振動和模態(tài)耦合型振動等，通常認為再生型振動為顫振發(fā)生的主要原因，如圖1所示。

圖1 二自由度銑削模型中再生效應示意圖

而在低速銑削區(qū)域，銑削穩(wěn)定性極大地受到源于刀具和工件間接觸處產生的過程阻尼效應的影響[34]。過程阻尼主要來源于刀具后刀面與工件波紋表面之間的干涉作用，通常被認為是一種能量轉移或消散機制，可以有效提高低速銑削的穩(wěn)定性，如圖2所示。

圖2 過程阻尼產生機理示意圖

Ahmadi等[34]指出，過程阻尼可以用等效黏性阻尼表示為

(1)

式中：Ceq為等效阻尼系數；Ksp為壓痕系數；W為刀具后刀面的磨損帶寬度；ap為軸向切削深度；切削速度v=πDΩ/60，D表示刀具直徑。

通過銑削動力學模型對再生型振動和過程阻尼進行建模，表達式如式(2)所示

(2)

式中：M、C和K分別是刀具的模態(tài)質量、模態(tài)阻尼和模態(tài)剛度矩陣；ap是切削過程中的軸向切深；t為切削的瞬時時刻；q(t)為刀具模態(tài)坐標；H(t)為周期切削力系數矩陣；G(t)為過程阻尼力系數矩陣；T為時滯量且等于刀齒切削周期，T=60/(NΩ)，且N為刀具齒數，Ω為刀具主軸轉速，單位為r/min。

通過考慮過程阻尼的銑削模型來預測銑削穩(wěn)定性，采用剛性刀具-柔性工件系統(tǒng)，以X和Y方向上具有相等的模態(tài)參數的雙自由度銑削動力學模型為例，可簡化為一個兩自由度銑削系統(tǒng)，如圖1所示。其銑削動力學模型可由式(2)表示為

其中，切削力系數表達式為

Knsin(φj(t))],

Knsin(φj(t))],

過程阻尼系數表達式為

gxx(t)=

gxy(t)=

gyx(t)=

gyy(t)=

式中：φst是刀具的切入角；φex是刀具的切出角。順銑時，φst=arccos(2a/D-1)，φex=π；逆銑時，φst=0，φex=arccos(1-2a/D)，其中a/D為徑向切深與刀具直徑之比。

(3)

式中，A(t)和B(t)為由再生效應和過程阻尼效應決定的周期系數矩陣。其中，

(4)

(5)

綜上，對于考慮再生顫振和過程阻尼的銑削穩(wěn)定性預測可以轉換為對式(3)的求解。

2 一種多步回溯算法

在一些難加工(如鈦合金)薄壁件的銑削穩(wěn)定性加工中，由于材料加工特性的影響，導致其加工的切削速度不易過高，進而在一定程度上限制了銑削加工的主軸轉速，而由于低速區(qū)域穩(wěn)定性邊界起伏變化比中速和高速區(qū)域快的多，即導致同一種算法的一個離散步內(假設均以單位1為離散步長)，高速區(qū)域達到精度要求時，在低速區(qū)域達不到精度要求，如圖3所示。

(a) 穩(wěn)定性邊界低速區(qū)示意圖

(b) 穩(wěn)定性邊界高速區(qū)示意圖

本文從這個角度出發(fā)，通過采用多步回溯插值算法，在確保算法在全局穩(wěn)定性預測精度的基礎上，提高其在主軸轉速低速區(qū)間的收斂速度，并以收斂率及計算精度作為算法的評價指標。在全離散法的基礎上，對狀態(tài)方程的時滯項進行多步插值近似，從而提出一種多步回溯算法來預測銑削過程中的穩(wěn)定性。

設銑削過程的初始時刻為t0，通過狀態(tài)空間變換方法，式(3)可以解為如下形式

(6)

式中：T為刀齒切削時間周期且等于時滯量；τ為將T等距間隔離散后的每個離散區(qū)間。

由于刀齒切削的時間周期T可分為兩個階段，即自由振動階段t∈[t0,t0+tf]和強迫振動階段t∈[t0+tf,t0+T]。當刀具處于自由振動階段，狀態(tài)值如下所示

x(t)=eA(t-t0)x(t0)

(7)

當刀具處于強迫振動階段時，即t∈[t0+tf,t0+T]，通過離散化的思想將強迫振動的時間等距分成m個時間區(qū)間，且每個時間區(qū)間的長度可以表示為τ=(T-tf)/m，則每個區(qū)間的時間點為

ti=t0+tf+(i-1)h(i=1，…，m+1)

(8)

當t∈[ti,ti+1]時，代入式(6)，原方程可化為

x(t)=eA(t-ti)x(ti)+

(9)

采用線性多步法的思想對狀態(tài)量x(t)進行插值計算，由于本文提出的多步回溯算法步數k(k≤4,且k∈N*)取值的不同會影響計算的精度和計算效率，故本節(jié)給出不同k取值時算法的計算式，并重點推導回溯步數k=4時算法的計算過程。

當t=t1時，即處于初始時刻時，狀態(tài)量

x(t1)=x(t0+tf)=eAtfx(t0)=eAtfx(tm+1-T)

(10)

為簡化表達式，用Ai+1表示A(ti+1)，則不同回溯步數k的第i個離散點算法計算式可表示為

(11)

xi=

(k=2)

(12)

xi=

(13)

xi=

(14)

本文以k=4為例，則式(14)可以被簡化為

xi=(I+Fi)-1[(F0+Fi-1)xi-1+Fi-2xi-2+

(15)

其中，F0=eAτ，

式中，

Bi1=Bi+(Bi-1-Bi)/τ，

Bi2=Bi-1+(Bi-2-Bi-1)/τ，

Bi3=Bi-2+(Bi-3-Bi-2)/τ，

Bi4=Bi-3+(Bi-4-Bi-3)/τ。

若式中[I-Fi]-1非奇異，令F=(I+Fi)-1，將式(15)代入到式(3)中可得離散映射為：yi+1=Diyi。

其中，向量yi可以定義為

yi=col(xi,…,xi-1,…,xi-m,…,xi-4-m)

式中，col(xi,…,xi-1,…,xi-m,…,xi-4-m)表示將括號里的列向量xi，…，xi-1，…，xi-4-m等集合成一個矩陣yi。

則Di可以定義為

其中，當k=4時，Di可表示為

綜上，在一個切削時間周期內的狀態(tài)傳遞矩陣D可以通過Floquet理論構建為

ym=Dm-1，…，D1D0y0=Dy0

(16)

即：D=Dm-1，…，D1D0

(17)

根據李雅普諾夫理論可知，系統(tǒng)的穩(wěn)定性與狀態(tài)傳遞矩陣的特征值有關，據此可獲得銑削系統(tǒng)的穩(wěn)定性葉瓣圖。

首先建立考慮過程阻尼效應的銑削動力學模型，通過將過程阻尼力看成是阻尼矩陣的附加項進行求解，并提出一種多步回溯算法計算模型的狀態(tài)傳遞矩陣，算法計算流程如圖4所示。

其中，eig()表示取特征值，abs()表示取絕對值，max()表示取最大值。

3 算法數值仿真分析

3.1 不同算法收斂率比較

圖4 基于多步回溯算法預測銑削穩(wěn)定性流程圖

同時，為揭示本文提出算法的收斂率，本節(jié)采用雙自由度銑削動力學模型進行計算，由于四階龍格庫塔法和完全離散法具有較高的收斂速度，下面選擇這兩種算法作為參照組與本文提出的算法進行比較，并將加工的系統(tǒng)參數設置為與參考文獻[18]中所采用的參數相同。本文所用的具體系統(tǒng)參數及加工參數設置如表1所示。得到不同算法不同離散數m的收斂率曲線如圖5所示。

表1 系統(tǒng)參數及加工參數

圖5揭示了不同算法之間收斂率的快慢。由圖5不難看出，在相同的切削條件下，所提出的多步回溯算法比四階龍格庫塔法和完全離散法具有更高的收斂效率。具體來說，當m=20時，多步回溯算法的收斂率已基本趨于平穩(wěn)，而圖中其它算法的收斂率在m=50時才基本趨于平穩(wěn)，且誤差值較m=20時的多步回溯算法大。

3.2 算法穩(wěn)定性預測誤差分析

在實際預測中，在低速銑削區(qū)域，由于過程阻尼的影響常常使得穩(wěn)定性區(qū)域進一步地擴大。常用的具有較高收斂速度的四階龍格庫塔法和完全離散法無法直接進行考慮過程阻尼的銑削動力學方程的求解。故本文選擇近兩年提出的二階半離散法作為對照組進行仿真計算，并將參數與參考文獻[33]設置相同，實驗參數如表2所示，得到穩(wěn)定性預測葉瓣圖如圖6所示。

(a) ap=0.6 mm, |λ0|=0.973 4

(b) ap=1.2 mm,|λ0|=0.093 7

表2 銑削穩(wěn)定性驗證實驗參數[33]

圖6(a)為引用參考文獻[33]二階半離散法的預測圖，圖6(b)為本文提出的多步回溯算法的預測圖，圖中的顫振點、邊界點、穩(wěn)定點均為文獻[33]中的實驗數據點。計算時間由MATLAB軟件在同一臺計算機[Intel?Core(TM) i5-3210M 2.50 GHz]上運算得出。

(a) 二階半離散法[33]m=80,時間t=2 138 s

(b) 多步回溯算法m=50,時間t=763 s

為進一步分析所提出算法的預測精度，本文引入了邊界預測誤差(Boundary Prediction Error, BPE)的概念。邊界預測誤差為預測邊界上出現實測顫振點和穩(wěn)定點占實測總數的百分比。

邊界預測誤差的計算如式(18)所示

(18)

式中：Ne為預測邊界出現實測非邊界點的個數；Nt為實測總點數。

從圖6(a)和6(b)中不難看出，文獻[33]中的二階半離散算法的預測邊界出現實測非邊界點的個數Ne=9，文中提出的多步回溯算法預測邊界出現實測非邊界點的個數Ne=4。根據式(18)可求出不同算法的BPE，如表3所示。

表3 考慮過程阻尼效應不同算法的計算時間及BPE

由表3可知，當多步回溯算法的離散數m=50時，其邊界預測誤差較二階半離散法減少了約12.5%，其計算時間也縮短了約64.3%。

4 結 論

本文重點研究了使用基于多步回溯算法對銑削穩(wěn)定性的預測，并從該研究中得出以下結論：

(1) 將動態(tài)銑削系統(tǒng)抽象為具有單個時間延遲的微分方程，通過微分方程理論可導出系統(tǒng)的一般解，并采用根據多步插值的思想對一般解進行離散化，最后采用Floquet理論確定銑削條件的穩(wěn)定性，最后獲得穩(wěn)定性葉瓣圖。

(2) 通過與四階龍格庫塔法和完全離散法在收斂速度與計算精度方面進行比較，結果表明：多步回溯算法可以在確保全局計算精度的前提下，對于低主軸轉速銑削區(qū)域具有更高的收斂速度，這種特性尤其適合一些銑削速度不高的難加工材料(例如鈦合金)的銑削預測。

(3) 通過多步回溯算法與考慮過程阻尼效應的二階全離散法的比較，結果表明：如果設置小于二階全離散法(m=80)的離散數，多步回溯算法(m=50)的邊界預測誤差較二階半離散法減少了約12.5%，其計算時間也縮短了約64.3%。

綜上所述，本文提出的多步回溯算法具有一定地快速收斂及高計算精度等特點，尤其在低速銑削的穩(wěn)定性預測中具有良好的應用前景。