黃 超，楊文安，蔡旭林，鄭方志，穆英娟，郭國強(qiáng)

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

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

Altintas 等[5]提出了一種零階半解析法來預(yù)測銑削過程中的穩(wěn)定性，但該方法無法預(yù)測小徑向切深下的穩(wěn)定切削區(qū)域。Merdol 等[6]提出了一種考慮高次諧波影響的多頻率穩(wěn)定性預(yù)測方法，該方法能滿足小徑向切深時的預(yù)測精度，但是計算效率相對較低。Insperger 等[7]首先提出了具有高精度的時域數(shù)值穩(wěn)定性預(yù)測方法，該方法僅對時滯項(xiàng)進(jìn)行離散處理?；谥苯臃e分思想，Ding 等[8]提出了一種全離散法穩(wěn)定性預(yù)測方法，全離散法利用線性插值同時逼近時滯微分方程中的狀態(tài)項(xiàng)和時滯項(xiàng)，極大地提高了該方法的計算效率。張偉等[9]將變齒距銑刀銑削作為顫振的抑制策略，通過試驗(yàn)驗(yàn)證變齒距銑刀具有較好的銑削顫振控制能力。Wan 等[10]提出了由于變齒距引發(fā)的多時滯銑削穩(wěn)定性預(yù)測模型，并基于半離散法求解變齒距銑削的穩(wěn)定性邊界。Zhang等[11]提出了一種改進(jìn)的全離散來求解由變齒距引起的多時滯問題的銑削穩(wěn)定性葉瓣圖，并指出銑削過程中刀具跳動會影響銑削穩(wěn)定性邊界。Jin 等[12]提出了一種考慮螺旋角的變齒距銑削改進(jìn)半離散法，并研究了不同變齒距銑刀對銑削顫振穩(wěn)定性的影響。劉寬等[13]利用AdvantEdge 有限元仿真和Matlab 軟件進(jìn)行了銑削力建模和分析，驗(yàn)證了變齒距銑刀具有降低振幅和減小切削力的效果。Jin 等[14]建立了變齒距和變主軸轉(zhuǎn)速的時滯微分運(yùn)動方程并運(yùn)用改進(jìn)半離散法求解，獲得了銑削穩(wěn)定性葉瓣圖。

針對目前變齒距銑削穩(wěn)定性預(yù)測算法計算效率低且通用性差的問題，在不降低預(yù)測精度的情況下，本文提出了一種基于精細(xì)積分的變齒距銑削穩(wěn)定性預(yù)測的全離散法，以快速精確地獲得變齒距銑削系統(tǒng)的穩(wěn)定性葉瓣圖。

1 數(shù)學(xué)模型

如圖1所示，銑削過程可以簡化為兩自由度振動“質(zhì)量–彈簧–阻尼”模型。假定刀具相對于工件是彈性體。根據(jù)牛頓運(yùn)動定律，銑削運(yùn)動方程可表示為

圖1 兩自由度銑削系統(tǒng)示意圖Fig.1 Schematic diagram of two-degree-of-freedom milling system

式中，M、C和K分別為模態(tài)質(zhì)量、阻尼和剛度矩陣；q(t)=[x(t)y(t)]T是位移向量；F(t)=[Fx(t)Fy(t)]T是切削力向量。

由于刀具考慮螺旋角效應(yīng)，實(shí)際參與切削的刀齒的高度時刻發(fā)生變化。如圖2所示，第j齒在微單元高度dz上沿切向和法向所受到的微元切削力表示為

圖2 變齒距螺旋刀具幾何圖Fig.2 Geometry of variable pitch spiral cutter

式中，Kt和Kr是切向和法向切削力系數(shù)；hj（t,z）是動態(tài)切削厚度，可定義為

式中，τj是時滯；fz是進(jìn)給量；?j（t,z）是第j齒隨時間變化的角度位置，定義為

式中，Ω、β和R分別是主軸轉(zhuǎn)速、刀具螺旋角和半徑；ψi是第j齒和j–1 齒之間的齒距角。

在t時刻，第j齒在x和y方向的切削力為

則刀具作用在x和y方向上總切削力為

式中，窗函數(shù)g（?j（t,z））為

式中，?st和?cx是切入和切出角。

聯(lián)立式（2）、（3）、（5）和（6），則可以得到式（1）切削力向量表達(dá)式

將式（8）帶入式（1），可得到考慮螺旋角的變齒距銑削動力學(xué)方程為

2 精細(xì)積分全離散法

令x(t)=[q(t),Mq(t)+(t)/2]，略去對銑削穩(wěn)定性無影響的靜態(tài)切削力部分，通過矩陣變換，則式（9）可改寫為狀態(tài)空間表達(dá)形式，即

將主軸旋轉(zhuǎn)周期T分為m等份，時間步長記為Δt=T/m，則第j–1 齒與j齒之間的時滯τj所對應(yīng)的離散間隔數(shù)mj可近似表示為[12,14]

式中，int(*) 是將正數(shù)四舍五入到0 的函數(shù)，比如int（3.76）=3。為了簡化運(yùn)算，令v（t）=A（t）x（t）和θj（t–τj）=Bj（t）x（t–τj）。在時間間隔[kτ,（k+1）τ]（k=0，1，…，m），對式（12）進(jìn)行常微分方程求解，初始條件為xk=x（tk），得到xk+1表達(dá)式為

式中，δ=ξ–kΔt，δ∈[0，Δt]。

其次，求解式（14）的積分項(xiàng)。在間隔[kτ,（k+1）Δt]中，式（14）狀態(tài)項(xiàng)v（δ）和時滯項(xiàng)θj（δ–τj（t））都利用線性插值多項(xiàng)式來近似表示，即

將式（15）和（16）代入式（14），得到

式中，

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

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

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

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

最后，將vk=Akxk，vk+1=Ak+1xk+1，θk–mj=Bj，k+1xk+1–mj和θk+1–mj=Bj，k+1xk+1–mj代入式（17），得

式（23）可獲得原方程的離散映射形式為

其中，

這里Zk+1=（I–T1hk+1Ak+1）–1。近似求得系統(tǒng)的狀態(tài)傳遞矩陣為Φ=Dm-1，Dm-2，…，D0。

依據(jù)Floquet 理論，假如傳遞矩陣Φ特征值的模小于1，系統(tǒng)穩(wěn)定；否則，系統(tǒng)不穩(wěn)定。

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

對于考慮螺旋角的變齒距銑刀銑削，采用文獻(xiàn)[12]提出的改進(jìn)半離散法對文中提出的方法開展理論驗(yàn)證研究。本文中數(shù)值模擬所有采用的仿真參數(shù)見表1[12]。所有運(yùn)算均基于Matlab9.2 平臺并在同一臺電腦上（Inter(R) Core(TM) i5-7400 CPU，8GB）完成。

在不同的徑向切深率（a/D）下，文中方法和改進(jìn)半離散法在表1所示切削條件下所獲得的穩(wěn)定性預(yù)測結(jié)果對比如圖3所示[12]?？芍?，本文方法與文獻(xiàn)[12]中所提的改進(jìn)半離散法所預(yù)測結(jié)果基本吻合。由于改進(jìn)半離散法的計算結(jié)果得到了很好的試驗(yàn)驗(yàn)證，也間接地驗(yàn)證了本文方法預(yù)測考慮螺旋角變齒距穩(wěn)定性預(yù)測的準(zhǔn)確性。本文方法在4 種不同的徑向切深的計算時間約為35s，而改進(jìn)半離散法約為320s，本文方法計算效率提高約90%。由此可知，本文方法具有更高的計算效率。

圖3 本文方法與改進(jìn)半離散法的葉瓣圖[12]Fig.3 Stability lobe diagrams of proposed method and updated semi-discretization method[12]

本文方法效率提升的原因：（1）當(dāng)銑削穩(wěn)定性預(yù)測時，式（19）中的矩陣T1只與主軸轉(zhuǎn)速有關(guān)，而與軸向切深無關(guān)，因此，在通過掃描軸向切深時，確定式（25）中轉(zhuǎn)移傳遞矩陣Dk的過程中，無須在不同的切削深度計算矩陣指數(shù)。然而，對于改進(jìn)半離散法，在掃秒軸向切深時，矩陣指數(shù)的計算是必需的，并且是計算4 個類似矩陣指數(shù)。如果將主軸轉(zhuǎn)速和切削深度組成的參數(shù)平面分為Ns×Nd網(wǎng)格點(diǎn)，通過改進(jìn)半離散法獲得穩(wěn)定性葉瓣圖所需要的矩陣指數(shù)計算次數(shù)為Ns×Nd×k×4，但是用本文方法，該次數(shù)只為Ns次，因此計算效率極大提高。（2）本文方法通過精細(xì)積分計算矩陣指數(shù)，以替代改進(jìn)半離散法中通過迭代公式計算矩陣指數(shù)，并且避免了求解矩陣指數(shù)時的求逆運(yùn)算，從而有效地減少了計算消耗時間。

4 時域仿真

時域仿真可以預(yù)測銑削過程中刀具切削力和位移，廣泛應(yīng)用于銑削穩(wěn)定性預(yù)測和模型驗(yàn)證。本節(jié)使用一種考慮跳刀情況下的變齒距銑削的時域仿真來驗(yàn)證所提出的算法[16]。所采用的仿真參數(shù)如表1所示[12]。通過繪制每轉(zhuǎn)一次采樣數(shù)據(jù) （“+”符號）來研究模擬響應(yīng)的穩(wěn)定性。

表1 刀具和系統(tǒng)模態(tài)參數(shù)[12]Table 1 Parameters of tool and system modal[12]

圖4[12]為文中提出方法和改進(jìn)半離散法[12]在徑向切深率a/D=0.5 和離散間隔數(shù)m=72 銑削條件下獲得的穩(wěn)定性葉瓣圖?？梢钥闯?，文中提出方法與文獻(xiàn)[12]中的改進(jìn)半離散法所得預(yù)測結(jié)果有一定差異。因此，為了驗(yàn)證哪一種方法預(yù)測結(jié)果更加精確，選取圖4中A（6500r/min，3mm）和B（8500r/min，5.8mm）2 個點(diǎn)進(jìn)一步研究。從圖5和6 可以看出，A點(diǎn)處于穩(wěn)定銑削狀態(tài)，因?yàn)槊哭D(zhuǎn)一次采樣數(shù)據(jù)接近一個固定點(diǎn)值。相反，點(diǎn)B處于不穩(wěn)定的銑削過程（圖7），因?yàn)橥讲蓸拥臄?shù)據(jù)分布在x和y方向的位移圖中（圖8）。因此本文方法具有更高的銑削穩(wěn)定性預(yù)測精度。

圖4 本文方法與改進(jìn)半離散法預(yù)測結(jié)果比較[12]Fig.4 Stability lobe diagrams of proposed method and updated semi-discretization method[12]

圖5 圖4中A點(diǎn)x方向（進(jìn)給）和y方向的位移時程Fig.5 Time history for the x direction (feed) and y direction displacements obtained from the point A in Fig.4

圖6 A點(diǎn)x和y方向的位移圖Fig.6 Plot of x and y direction displacements obtained from the point A

圖7 圖4中B點(diǎn)x方向（進(jìn)給）和y方向的位移時程Fig.7 Time history for the x direction (feed) and y direction displacements obtained from the point B in Fig.4

圖8 B點(diǎn)x和y方向的位移圖Fig.8 Plot of the x and y direction displacements obtained from the point B

5 結(jié)論

針對變齒距銑削過程中變齒距引起的多時滯問題，提出了一種基于精細(xì)積分全離散法的變齒距銑削穩(wěn)定性的預(yù)測方法?；诙杂啥瓤紤]螺旋角變齒距銑削模型，通過與改進(jìn)半離散法對比及時域仿真，驗(yàn)證了本文方法的有效性。而且比較結(jié)果顯示本文方法具有更高的計算效率和預(yù)測精度，在同等仿真條件下，本文方法計算效率可以提高約90%，擁有更高的使用價值。