基于三階牛頓–埃爾米特插值全離散法的銑削穩(wěn)定性預(yù)測(cè)*

黃 超，楊文安，黃久超，解菲菲，楊有成

（1.上海航天精密機(jī)械研究所，上海 201600；2.南京航空航天大學(xué)機(jī)電學(xué)院，南京 210016）

銑削作為一種高效的切削加工方法，具有加工精度高、材料去除率大和加工成本低等優(yōu)點(diǎn)。在銑削過程中，主要存在的振動(dòng)形式有自由振動(dòng)、受迫振動(dòng)及自激振動(dòng)，其中自激振動(dòng)中的再生顫振是引起加工過程不穩(wěn)定的主要因素[1–2]。銑削顫振通常會(huì)導(dǎo)致加工零件表面質(zhì)量差、刀具磨損加劇以及數(shù)控機(jī)床壽命降低等問題[3]。此外，考慮再生顫振的銑削動(dòng)力學(xué)模型可以近似地描述為時(shí)滯微分方程，通過求解時(shí)滯微分方程可以獲得穩(wěn)定性葉瓣圖，達(dá)到銑削穩(wěn)定性預(yù)測(cè)的目的[4]。因此，如何快速精確地獲取銑削穩(wěn)定性葉瓣圖顯得尤為重要。

迄今為止，預(yù)測(cè)銑削穩(wěn)定性的方法主要分為數(shù)值法[5–6]和解析法[7–16]。Davies等[5]提出了一種改進(jìn)的時(shí)域計(jì)算方法預(yù)測(cè)出了小徑向切深切削下的額外不穩(wěn)定區(qū)域。Li等[6]提出了基于Runge–Kutta法的銑削穩(wěn)定性時(shí)域分析方法，但計(jì)算效率太低。Altintas等[7]提出了一種經(jīng)典的零階解析方法，但是該方法不能預(yù)測(cè)小徑向切深切削下的銑削穩(wěn)定性。Merdol等[8]提出了考慮了周期力系數(shù)矩陣的高階展開項(xiàng)多頻率法，實(shí)現(xiàn)了小徑向切深切削下的穩(wěn)定性預(yù)測(cè)。Insperger等[9]通過對(duì)時(shí)滯項(xiàng)進(jìn)行離散處理，提出了具有高精度的時(shí)域數(shù)值穩(wěn)定性預(yù)測(cè)方法。Ding等[10–11]基于直接積分思想提出了包括一階全離散法（1stFDM）和二階全離散法的穩(wěn)定性預(yù)測(cè)方法。全離散法利用線性插值同時(shí)逼近時(shí)滯微分方程中的狀態(tài)項(xiàng)和時(shí)滯項(xiàng)，極大提高了方法計(jì)算效率。Ding等[12]隨后提出了一種基于牛頓–柯特斯公式的數(shù)值積分法，該方法通過直接離散積分方程的積分項(xiàng)構(gòu)造出相鄰周期的轉(zhuǎn)移矩陣來預(yù)測(cè)銑削穩(wěn)定性。隨后，Zhang等[13]在此基礎(chǔ)上通過辛普森公式提出了一種全新的銑削穩(wěn)定性預(yù)測(cè)方法。最近，Yan等[14]基于全離散法，提出了一種改進(jìn)的三階全離散法（3rdUFDM），得到了比一階、二階全離散法更高的預(yù)測(cè)精度和計(jì)算效率。Dai等[15]提出了一種基于黃金比例搜索的改進(jìn)全離散方法來預(yù)測(cè)銑削穩(wěn)定性，并且引入精細(xì)積分法來計(jì)算矩陣指數(shù)，極大地提高了該方法的計(jì)算效率。Li等[16]利用泰勒公式對(duì)非齊次項(xiàng)二次展開，提出了顯式精細(xì)積分法預(yù)測(cè)銑削穩(wěn)定性。

本文基于再生顫振理論構(gòu)建的銑削動(dòng)力學(xué)模型，提出一種高精度、高效率的預(yù)測(cè)銑削穩(wěn)定性的三階牛頓–埃爾米特插值全離散法。

1 動(dòng)力學(xué)模型

1.1 單自由度動(dòng)力學(xué)模型

一般來說，考慮再生顫振的單自由度銑削動(dòng)力學(xué)模型表述為時(shí)滯微分方程[9]：

式中，mt，ζ和ωn分別表示刀具的模態(tài)質(zhì)量、刀具的阻尼比和自然頻率；w為軸向切削深度；x(t)為刀具模態(tài)坐標(biāo)；時(shí)滯量T等于刀齒切削周期，即T=60/(NΩ)，N為刀具齒數(shù)，Ω為主軸轉(zhuǎn)速（r/min）；切削力系數(shù)h(t)表示為：

式中，Kt和Kn分別為切向和徑向切削力系數(shù)；?j(t)為第j個(gè)刀齒轉(zhuǎn)角：

函數(shù)g(?j(t))定義為：

其中，?st和?ex分別是刀齒的切入與切出角。對(duì)于順銑，?st=arcos(2a/D–1)，?ex=π；對(duì)于逆銑，?st=0，?ex=arcos(1–2a/D)，其中a/D為徑向切深與刀具直徑之比。

令x(t)=[x(t)mx.(t)+mζωnx(t)]T，單自由度銑削模型的狀態(tài)空間形式為：

1.2 雙自由度銑削模型

不失一般性，考慮再生顫振的雙自由度銑削動(dòng)力學(xué)模型表述為時(shí)滯微分方程[9]：

式中，mt，ζ和ωn分別為刀具的模態(tài)質(zhì)量、刀具的阻尼比和自然頻率；w為軸向切削深度；并假設(shè)它們?cè)趚和y兩個(gè)方向是相等的，hxx(t)，hxy(t)，hyx(t)和hyy(t)為如下函數(shù)：

同樣地，雙自由度銑削模型的狀態(tài)空間形式為：

其中，

2 三階牛頓–埃爾米特插值全離散法

由前文可知，考慮再生顫振的銑削動(dòng)力學(xué)方程可由以下n維狀態(tài)空間方表述：

式中，A0是表示系統(tǒng)時(shí)不變性質(zhì)的常數(shù)矩陣，B(t)為周期系數(shù)矩陣，并滿足B(t+T)=B(t)。

三階牛頓–埃爾米特插值全離散法的第1步是將刀齒切削周期T離散為m個(gè)相等的時(shí)間區(qū)段，即T=mτ（其中m為正整數(shù)）。在每一個(gè)時(shí)間區(qū)段kτ≤t≤(k+1)τ，（k=0，…，m）上，以xk=x(kτ)為初始條件，式（9）的響應(yīng)可以表示為：

令f(ξ)=B(ξ)x(ξ)，v(ξ)=B(ξ)x(ξ–T)，式（10）可以等效表示為

其中，δ=ξ–kτ，δ∈[0，τ]。求解式（11）中的積分項(xiàng)時(shí)首先構(gòu)造牛頓插值多項(xiàng)式：

其中，系數(shù)ai（i=0，1，…，n）由插值條件得

由插值條件Qn(x0)=f(x0)得：

由插值條件Qn(x1)=f(x1)得：

一般可得：

這樣得到滿足插值條件Qn(x1)=f(x1)，(i=0，1，…，n)的n次牛頓插值多項(xiàng)式

其中，

其中，Bk表示B(t)在時(shí)間點(diǎn)tk=kτ的取值。然后構(gòu)造埃爾米特插值多項(xiàng)式，給定f(x)在n+1個(gè)節(jié)點(diǎn)a=x0＜x1＜…＜xn=b上的函數(shù)值f(x0)，f(x1)，…，f(xn)，以及一階導(dǎo)數(shù)值f'(x0)，f'(x1)，…，f'(xn)。則可在每個(gè)小區(qū)間[xi，xi+1]上構(gòu)造3次埃爾米特插值多項(xiàng)式

其中，

所以，在時(shí)間區(qū)間[(k–m)τ，(k+1–m)τ]中，時(shí)滯項(xiàng)v(δ–T)可通過埃爾米特插值獲得：

其中，

將式（18）～（25）代入式（11），簡化的銑削系統(tǒng)表達(dá)式為：

其中，

基于矩陣指數(shù)的加法定理，矩陣指數(shù)T1可寫為

其中，dt=Δt/2n。為了提高計(jì)算精度，一般取n=20，精細(xì)區(qū)段dt就已經(jīng)非常小的區(qū)段了[17]。當(dāng)dt非常小時(shí)，可采用Taylor級(jí)數(shù)展開近似矩陣指數(shù)，即

將式（30）代入式（29），矩陣指數(shù)T1近似解為

執(zhí)行以下運(yùn)算，

由式（26）可以構(gòu)造如下離散映射

其中，

然后，系統(tǒng)在單個(gè)時(shí)間周期上的狀態(tài)轉(zhuǎn)移矩陣Ф定義為

最后，根據(jù)Floquet理論，若系統(tǒng)轉(zhuǎn)移矩陣Ф特征值都小于1，則系統(tǒng)穩(wěn)定；否則，系統(tǒng)發(fā)生顫振。

3 數(shù)值仿真結(jié)果分析

3.1 收斂性分析

局部離散誤差分析可以反映所提方法近似解收斂速度的快慢。將局部離散誤差||μ|–|μ0||表示為周期離散數(shù)m的函數(shù)，其中參考值|μ0|是m=500時(shí)采用1stFDM得到的特征值的模。以單自由度銑削動(dòng)力學(xué)模型為例，對(duì)本文所提出方法與1stFDM和3rdUFDM的收斂性進(jìn)行分析。本文中數(shù)值仿真所采用的仿真參數(shù)均與文獻(xiàn)[14]相同。切削參數(shù)為：切向力系數(shù)Kt=6×108N/m2，法向力系數(shù)Kn=2×108N/m2，徑向切深比a/D=1，主軸轉(zhuǎn)速Ω=5000r/min，軸向切削深度w=0.2mm、0.5mm、0.7mm、1mm。模態(tài)參數(shù)：順銑，刀具齒數(shù)N=2，模態(tài)質(zhì)量mt=0.03993kg，模態(tài)阻尼ζ=0.011，固有頻率ωn=922×2πrad/s。本文所有仿真均在臺(tái)式機(jī)電腦（Inter(R) Core(TM) i5-7400 CPU, 8GB）上運(yùn)行，運(yùn)行平臺(tái)為MATLAB 9.2，操作系統(tǒng)為Windows 10。

圖1給出了1stFDM、3rdFDM和文中方法在4種不同軸向切削深度時(shí)的局部離散誤差變化趨勢(shì)。由圖1可知，隨著離散數(shù)m的增加，1stFDM、3rdUFDM和文中方法的局部離散誤差逐漸接近于零。在離散數(shù)相同的情況下，文中方法的局部離散誤差最小。比如主軸轉(zhuǎn)速Ω=5000r/min，軸向切削深度w=0.7mm，離散數(shù)m=45時(shí)，1stFDM、3rdUFDM和文中方法的局部離散誤差分別為0.043，0.014和0.003，文中方法的計(jì)算精度比1stFDM和3rdUFDM的計(jì)算精度分別提高約93%和79%。在相同的局部離散誤差下，文中方法所需的周期離散數(shù)更小。因此，文中方法的計(jì)算效率遠(yuǎn)遠(yuǎn)高于1stFDM和3rdUFDM。

文中方法兩自由度的局部離散誤差變化趨勢(shì)如圖2所示。兩自由度銑削過程的系統(tǒng)參數(shù)取值的大小與單自由度的參數(shù)相同。由圖2可知，兩自由度模型的收斂速度和單自由度模型的收斂速度一樣快，且當(dāng)離散誤差m=40時(shí)都接近于0。

3.2 穩(wěn)定性葉瓣圖

3.2.1 單自由度銑削模型

在單自由度銑削模型中，為了比較穩(wěn)定性邊界的計(jì)算精度，分別使用1stFDM、3rdUFDM和文中方法獲得穩(wěn)定性葉瓣圖。將模態(tài)參數(shù)設(shè)置與3.1節(jié)相同，主軸轉(zhuǎn)速5000≤Ω≤10000r/min，軸向切削深度0≤w≤10mm，徑向比a/D=1，順銑，主軸轉(zhuǎn)速與軸向切削深度平面的網(wǎng)格數(shù)為200×100。當(dāng)離散數(shù)m足夠大時(shí)，1stFDM獲取的穩(wěn)定性邊界接近于真實(shí)值[10]。因此，取離散數(shù)m=200，1stFDM獲取的穩(wěn)定性葉瓣圖作為參考并在圖中以紅色給出。

圖1 1stFDM、3rdUFDM和本文方法收斂速度的比較Fig.1 Convergence rate comparison of 1stFDM, 3rdUFDM and proposed method

圖3給出了在不同離散數(shù)下3種方法所得的穩(wěn)定性葉瓣圖與參考曲線的對(duì)比結(jié)果，同時(shí)給出了計(jì)算時(shí)間?？芍S著離散數(shù)m的增加，3種方法的計(jì)算精度逐漸提高。顯然當(dāng)離散數(shù)m相同時(shí)，文中方法的計(jì)算精度比1stFDM和3rdUFDM的計(jì)算精度都高，而且計(jì)算時(shí)間都是最短的。所以，文中方法是一種具有較高的計(jì)算精度和計(jì)算效率的銑削穩(wěn)定性預(yù)測(cè)方法。

3.2.2 雙自由度銑削模型

雙自由度銑削過程的系統(tǒng)參數(shù)取值的大小與單自由度的參數(shù)相同。周期離散數(shù)m取為40，主軸轉(zhuǎn)速與軸向切削深度平面劃分為200×100網(wǎng)格點(diǎn)，徑向切深比a/D=0.05，0.1，0.5順銑。在雙自由度銑削模型中，1stFDM、3rdUFDM和文中方法的銑削穩(wěn)定性葉瓣圖計(jì)算結(jié)果和時(shí)間列于圖4中。可知，當(dāng)a/D=0.5時(shí)，在沒有損失任何計(jì)算精度的前提下，文中方法的計(jì)算時(shí)間較1stFDM縮短21%，較3rdUFDM縮短9%。故文中方法是一種更加高效的銑削穩(wěn)定性預(yù)測(cè)方法。

4 試驗(yàn)驗(yàn)證

針對(duì)本文提出的用于預(yù)測(cè)銑削穩(wěn)定性的算法，本節(jié)通過試驗(yàn)對(duì)銑削穩(wěn)定性預(yù)測(cè)結(jié)果進(jìn)行驗(yàn)證。Gradisek等[18]已經(jīng)通過試驗(yàn)驗(yàn)證了半離散法的有效性，所以本次驗(yàn)證試驗(yàn)所用的模態(tài)參數(shù)和切削試驗(yàn)數(shù)據(jù)均與文獻(xiàn)[18]相同，仿真結(jié)果和試驗(yàn)結(jié)果如圖5所示。

可以看出，當(dāng)徑向比a/D取1、0.5和0.1時(shí)，文中方法所獲得的穩(wěn)定性葉瓣圖很好地預(yù)測(cè)出了銑削系統(tǒng)的穩(wěn)定區(qū)域。對(duì)于小徑向銑削，即a/D=0.05，文中方法所獲得的穩(wěn)定性葉瓣圖和試驗(yàn)結(jié)果還是吻合的，但是在主軸轉(zhuǎn)速11000≤Ω≤17000r/min時(shí)預(yù)測(cè)誤差較大，這可能是由于本文使用的銑削動(dòng)力學(xué)模型沒有考慮實(shí)際銑削過程中存在的某些非線性因素，導(dǎo)致上述穩(wěn)定性葉瓣圖與實(shí)際試驗(yàn)之間存在誤差。綜上可見，本文方法是一種高效準(zhǔn)確的預(yù)測(cè)銑削穩(wěn)定性的方法。

5 結(jié)論

針對(duì)考慮再生顫振的銑削動(dòng)力學(xué)模型，本文提出了三階牛頓–埃爾米特插值全離散法預(yù)測(cè)銑削穩(wěn)定性。

圖2 本文方法兩自由度收斂速度Fig.2 Convergence rate of the proposed method for two degree-of-freedom

圖3 1stFDM、3rdUFDM和文中方法對(duì)于單自由度模型的穩(wěn)定性葉瓣圖Fig.3 Stability lobes of 1stFDM, 3rdUFDM and proposed method in this paper for a single degree of freedom model

圖4 1stFDM、3rdUFDM和文中方法對(duì)于雙自由度模型的穩(wěn)定性葉瓣圖Fig.4 Stability lobes of 1stFDM, 3rdUFDM and proposed method in this paper for two-degree-of-freedom model

圖5 文中方法的穩(wěn)定性預(yù)測(cè)結(jié)果和試驗(yàn)驗(yàn)證結(jié)果Fig.5 Predicted stability boundaries and stability obtained from experiment

（1）通過對(duì)局部離散誤差的分析可知，文中方法的計(jì)算精度比1stFDM和3rdUFDM的計(jì)算精度分別提高約93%和79%。隨著軸向切削深度增加，文中方法在較少的離散數(shù)時(shí)能達(dá)到1stFDM和3rdUFDM較多離散數(shù)時(shí)的局部離散誤差。

（2）在單自由度銑削模型中，隨著離散數(shù)的增加，以上3種方法的計(jì)算精度都逐漸提高。當(dāng)離散數(shù)相同時(shí)，文中方法的計(jì)算精度高于1stFDM和3rdUFDM，并且計(jì)算效率也是最高的。在雙自由度銑削模型中，當(dāng)離散數(shù)相同時(shí)，3種方法的計(jì)算精度相當(dāng)，但文中方法的計(jì)算效率比1stFDM提高約21%，比3rdUFDM提高約9%。

（3）試驗(yàn)表明，文中方法可以高效準(zhǔn)確地預(yù)測(cè)銑削穩(wěn)定性。