楊赫然，耿超緒，孫興偉，董祉序

(沈陽工業(yè)大學 機械工程學院，沈陽 110870)

加工過程的振動及抑制策略直接影響螺桿轉子銑削加工質量和加工精度，并可對機床及刀具產生嚴重損害[1]。加工過程中的振動形式主要有再生顫振和強迫振動兩種[2]。再生顫振是一種強烈的自激振動，亦是高性能機床在加工過程中必須控制的主要振動形式[3]。按照產生機理[4]，顫振可以分為再生型顫振、振型耦合型顫振、摩擦型顫振、力-熱耦合型顫振等。

其中，再生型顫振往往先于其他形式發(fā)生，也是金屬切削過程中主要的顫振形式[5]。

為了防止螺旋曲面銑削加工過程中再生顫振對加工的影響，加工前有必要對螺旋曲面銑削穩(wěn)定性進行預測，以保證加工時能夠延長刀具的使用壽命并且獲得較高的表面質量[6-7]。

針對再生型顫振對精工的影響，通常用線性時滯微分方程來表示再生型顫振影響的銑削加工動力學模型[8]。所以如何能夠獲取時滯微分方程精確的解，且獲得穩(wěn)定性葉瓣圖，是預測螺旋曲面銑削系統(tǒng)穩(wěn)定性的關鍵。

到目前為止，預測銑削穩(wěn)定性的方法有很多種，主要分為時域法、頻域法和試驗方法。

時域法包括半離散法、全離散法和時間有限元方法等。Insperger等[9-11]提出半離散法，該方法對時滯微分方程中的時滯項進行離散，其算法的精度與離散步長呈正相關，基于半離散法，又提出改進半離散法以及高階半離散法。Niu等[12]提出基于Runge-Kutta的分析方法。Li等[13]提出了另一類全離散法，該方法離散了所有時間相關項，利用數(shù)值迭代方法得到迭代公式，從而得到Floquet轉移矩陣。Ding等[14]基于時間有限元法，以積分方程技術求解銑削動態(tài)系統(tǒng)響應，提出數(shù)值積分法，并將其發(fā)展為具有指數(shù)收斂階的譜方法和針對多時滯工況的變步長法[15]。Ding等[16-17]提出全離散方法，包括一階全離散法和二階全離散法，一階全離散是通過線性插值對狀態(tài)項和時滯項進行離散，二階全離散法是通過Lagrange插值對狀態(tài)項進行離散，通過線性插值對時滯項進行離散。Guo等[18]通過牛頓插值提出了三階全離散法，能夠得到更高的計算精度和更快的計算效率。Zhang等[19]基于Simpson公式提出了一種全新的方法預測銑削穩(wěn)定性。

頻域法一般為將顫振模型的時滯微分方程組利用傅里葉變換轉到頻域表示，基于控制理論，解析計算銑削穩(wěn)定邊界。Budak等[20-21]提出了頻域法求解銑削動力學方程，但這種方法的計算精度無法滿足小徑向切深。Merdol等[22]提出了多頻率法，該方法考慮了定向因子的高次諧波，然而在計算過程中需要迭代搜索顫振頻率，需求解多個特征值。Bachrathy等[23]將多頻率法擴展到適用于所有刀具結構的穩(wěn)定性預測，包括可引入分布時滯的復雜刀具結構，將擴展多頻率法和多維二分法結合以提高計算效率，并證明了在所測頻響函數(shù)質量較差的情況下，擴展多頻率法依然可以得到可靠的穩(wěn)定預測結果。

綜上所述，眾多學者針對銑削加工過程的穩(wěn)定性預測進行了較深入的研究，然而隨著新算法的涌現(xiàn)，在預測精度及速度上仍有較大的提升空間，且針對不同動力學系統(tǒng)的適應性有進一步探究的必要。

Adams作為典型的線性多步法代表，當插值區(qū)間較大，計算步數(shù)較多時，由于其參與運算的函數(shù)個數(shù)較其他算法少，因此在計算效率方面具有明顯優(yōu)勢。智紅英等[24]針對單自由度及雙自由度高速銑削系統(tǒng)進行預測研究，也表明該方法在效率及精度方面優(yōu)勢明顯。螺桿轉子銑削動力學模型為三自由度動力學模型，需要在保證預測精度的同時提高預測速度，因此，本文基于隱式Adams方法對螺旋曲面銑削系統(tǒng)的穩(wěn)定性進行預測，探究隱式Adams方法對螺桿轉子銑削動力學模型的適用性，為螺旋曲面銑削加工穩(wěn)定性預測提供理論依據。

1 動力學模型建立

對于多自由度系統(tǒng)，預測銑削穩(wěn)定性的再生顫振動力學模型，可以用時滯微分方程進行表示，如式(1)所示

(1)

式中：M為模態(tài)質量矩陣；C為模態(tài)阻尼矩陣；K為模態(tài)質剛度矩陣；q(t)為模態(tài)坐標矩陣；H(t)為周期系數(shù)矩陣；Tc為刀齒周期(時滯量)，Tc=60/(Nn)；N為刀具齒數(shù)；n為主軸轉速；ap為切削深度。

1.1 三維銑削力模型

在銑削螺桿轉子曲面時，由于數(shù)控系統(tǒng)中給出的進給量沿著工件Zw軸方向，而刀具沿著Xc軸上下移動，且刀具存在擺角，所以總的進給量由Zw向及Xc軸復合而成。因此，首先要建立銑削力模型，如圖1所示[25]。

圖1 銑削力模型Fig.1 Milling force model

以瞬時剛性切削理論為依據，銑削過程中的瞬時切削厚度、切向銑削力，徑向銑削力、軸向銑削力的表達式如式(2)和式(3)所示

(2)

(3)

式中：ΔXc，ΔYc和ΔZc分別為Xc，Yc和Zc方向上盤銑刀銑削螺桿轉子曲面時產生的動態(tài)位移；ft為橫向進給率；fT為軸向進給率。由于決定銑削穩(wěn)定性的是動態(tài)銑削力，所以忽略靜態(tài)銑削力，將切向銑削力Ftl、徑向銑削力Frl和軸向銑削力Fal分別映射到刀具坐標的Xc，Yc，Zc，可以得到三個方向的力如式(4)所示

(4)

1.2 三自由度動力學模型

隨后在建立三自由度力學模型的-基礎上建立動力學模型，如圖2所示。

圖2 三自由度動態(tài)銑削動力學模型Fig.2 Three-degree-of-freedom dynamic milling dynamic model

則銑削動力學模型可以由下列方程進行表示

(5)

式中：mtx，mty和mtz為刀具Xc，Yc和Zc方向的模態(tài)質量；ξx，ξy和ξz為Xc，Yc和Zc方向的阻尼比；ωnx，ωny和ωnz為Xc，Yc和Zc方向的固有圓頻率。

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

式中：Kt為切向切削力系數(shù)；Kr為徑向切削力系數(shù)；Ka為軸向切削力系數(shù)；ωl(t)為第l個刀齒的位置角。

(15)

式中：ωst為刀具的切入角；ωex為刀具的切出角。對于順銑ωst=arccos(2b/D-1)，ωex=π；對于逆銑，ωst=0，ωex=arccos(1-2b/D)，其中b/D為徑向切深和直徑的比值。

(16)

其中，

(17)

(18)

2 參數(shù)識別

若要精確的求解動力學模型，需要獲取盤銑刀的模態(tài)參數(shù)，介于盤銑刀結構特點，僅通過試驗不易獲得三個方向的模態(tài)參數(shù)，因此本章采用模態(tài)試驗與仿真結合的方法，獲取盤銑刀的模態(tài)參數(shù)。對裝配狀態(tài)下的盤銑刀進行模態(tài)分析，獲得各階振型。

2.1 盤銑刀裝配狀態(tài)下的模態(tài)仿真分析

為精確分析盤銑刀動態(tài)特性，采用有限元分析軟件對裝配狀態(tài)下的盤銑刀部件進行模態(tài)分析，得到其前六階振型，其中前兩階如圖3所示。前六階振型的變形方向均圍繞Z軸，且X、Y、Z向模態(tài)參數(shù)相同。

圖3 裝配狀態(tài)下的盤銑刀前兩階模態(tài)振型Fig.3 Mode shape of disc milling cutter in assembly state

2.2 盤銑刀裝配狀態(tài)的諧響應分析

在模態(tài)仿真分析的基礎上進行諧響應分析可獲取其位移響應所集中的位置，同時為模態(tài)試驗提供支撐依據。根據振型圖，在盤銑刀形變位移最大的位置，選取三個位置分別施加300 N的簡諧作用力，方向沿著刀具Z軸負方向。將諧響應分析求解的最大頻率范圍設置為10 000 Hz，求解間隔設置為400，基于2.1節(jié)模態(tài)分析進行求解，得到盤銑刀的時間位移響應，如圖4所示。

通過圖4可知，Xc方向的響應頻率峰值出現(xiàn)在1 725 Hz，3 925 Hz和4 400 Hz附近，Yc方向的響應頻率峰值出現(xiàn)在1 725 Hz，1 800 Hz，3 925 Hz和4 400 Hz附近，Zc方向的響應頻率峰值出現(xiàn)在1 725 Hz，3 925 Hz和4 400 Hz附近。由模態(tài)分析得到的盤銑刀裝配狀態(tài)下1階～6階的固有頻率分別為：1 724.9 Hz，1 725.9 Hz，1 788.7 Hz，3 919.6 Hz，4 430.5 Hz和4 432.2 Hz，由于盤銑刀存在重根模態(tài)的影響，故會出現(xiàn)重復頻率。因此在盤銑刀裝配狀態(tài)下諧響應分析0～10 000 Hz頻率范圍內，盤銑刀裝配狀態(tài)下1階～6階固有頻率容易激發(fā)。

圖4 Xc，Yc，Zc位移-頻率響應曲線Fig.4 Frequency response curve in Xc,Yc and Zc directions of the assembled state of disc milling cutter

2.3 模態(tài)參數(shù)獲取試驗

本次試驗采用錘擊法獲取盤銑刀的刀尖頻響函數(shù)曲線，試驗采用DEWEsoft公司的SIRIUS-ACC數(shù)據采集器，力錘和加速度傳感器采用IEPE型，試驗過程采用移動力錘法。根據模態(tài)仿真分析可知，盤銑刀在裝配狀態(tài)下前6階的模態(tài)振型不同且出現(xiàn)重根模態(tài)，根據模態(tài)振型位移最大位置布置傳感器。傳感器在盤銑刀的布置點以及試驗臺布局分別如圖5和圖6所示。圖5中，按順時針方向對傳感器進行編號分別為1號、2號、3號和4號傳感器，根據傳感器布點位置，力錘敲擊產生振動后即可獲得所關心頻率范圍內的傳遞函數(shù)。

圖5 傳感器布置示意圖Fig.5 Schematic diagram of sensor arrangement

圖6 試驗裝置示意圖Fig.6 Schematic diagram of experimental device

由頻響函數(shù)參數(shù)估計的基本方程中留數(shù)和模態(tài)振型之間的關系和模態(tài)振型矩陣的互易性可知，矩陣的每一列元素，都包含模態(tài)振型的信息，因此只要模態(tài)振型在Xc，Yc和Zc方向存在與參考自由度相關的振型，則測量數(shù)據就能觀察到模態(tài)信息。

每隔一個刀齒選取一個錘擊點，盡可能多的選取錘擊點數(shù)，試驗采樣頻率設置為10 000 Hz，譜線數(shù)設置為8 192條，根據Nyquist理論，所能獲取的頻率窗口長度為采樣頻率的一半，經計算得采樣頻率的分辨率為0.61 Hz。根據傳感器布點，試驗獲取的傳遞函數(shù)曲線及相干性曲線如圖7和圖8所示。通過圖7的幅頻特性曲線可以獲取模態(tài)的固有頻率以及阻尼比等模態(tài)參數(shù)信息。模態(tài)試驗獲取的相干性曲線，反應了響應是否和激勵相關，當響應是由激勵引起時，相干性曲線的值接近1，則說明兩者相關。反之，則不相關。

圖7 傳遞函數(shù)幅頻特性曲線Fig.7 Amplitude-frequency characteristic curve of transfer function

圖8 相干性曲線Fig.8 Coherence curve

根據模態(tài)仿真分析結果、模態(tài)試驗數(shù)據以及傳感器布點位置綜合分析可知，試驗數(shù)據中的一階頻率對應模態(tài)仿真中的第三階固有頻率，其中的二階頻率對應模態(tài)仿真分析中的第四階固有頻率。

除1號和2號傳感器測量的數(shù)據外，3號和4號傳感器的數(shù)據分別與1號和2號的數(shù)據非常接近，說明盤銑刀具有良好的對稱性，但動態(tài)特性又稍有區(qū)別，并非完全對稱。將試驗數(shù)據與模態(tài)仿真數(shù)據進行對比分析可知，兩者的誤差均在10%以內，如表1所示。由此可知仿真結果具有較高可信度。

表1 盤銑刀裝配狀態(tài)模態(tài)試驗與模態(tài)仿真結果對比Tab.1 Comparison of assembly state modal experiment and simulation results of disc milling cutter

2.4 基于導納圓擬合法的模態(tài)參數(shù)識別

導納圓擬合法是一種比較直觀的方法，可以在數(shù)據質量稍低時，在一定數(shù)據范圍內搜索出峰值，擬合出峰值頻率。

根據傳遞函數(shù)幅頻特性曲線以及相干性曲線，可以判斷出本次試驗數(shù)據良好。根據模態(tài)試驗數(shù)據以及導納圓擬合方法分析得到的數(shù)據如表2所示，表2中固有頻率、阻尼比、模態(tài)質量和模態(tài)剛度對應仿真數(shù)據中的第一、第二和第四階模態(tài)信息。由模態(tài)仿真分析可知，盤銑刀X、Y、Z三向的模態(tài)質量、模態(tài)阻尼及模態(tài)剛度分別對應相等，因此以表2中的一階模態(tài)參數(shù)作為后續(xù)穩(wěn)定性預測的計算數(shù)據。

表2 盤銑刀試驗模態(tài)參數(shù)Tab.2 Test modal parameters of disk milling cutter

2.5 平均切削力系數(shù)識別

刀齒在一個銑削周期的平均銑削力可以表示為

(19)

由于顫振和表面粗糙度等因素的影響，使刀具的每個旋轉周期測得的力都存在一定差值，因此選取刀具切削螺桿轉子1/5的切削力，即5頭螺桿切削過程的1/5部分所測得的瞬時銑削力數(shù)據的平均值作為平均銑削力，使用定積分計算不同的銑削加工參數(shù)下的平均銑削力數(shù)據，50 mm棒料作為試驗材料，試驗加工參數(shù)由正交試驗法獲得，如表3所示。

表3 螺桿轉子曲面銑削試驗方案Tab.3 Experimental plan for screw rotor milling

由正交試驗表1可知，需進行10組螺桿轉子銑削試驗。其中：a為盤銑刀每加工一個截面沿工件Zw方向的位移量；n為盤銑刀轉速；F為加工倍率。根據試驗方案，測得不同銑削參數(shù)下的銑削力試驗數(shù)據，使用坐標變換式將銑削力投影到刀具坐標下，得到刀具坐標系Oc-XcYcZc下三個方向所受的平均銑削力，如表4所示。

表4 刀具平均切削力Tab.4 Average cutting force of tool

將表4所示的平均銑削力帶代入式(19)中，分別計算出中各組銑削參數(shù)下的平均銑削力系數(shù)如表5所示,系數(shù)中的負號表示該方向的銑削力與文中規(guī)定的正方向相反。

表5 螺桿轉子平均銑削力系數(shù)計算結果Tab.5 Calculation result of average milling force coefficient of screw rotor

通過多元非線性回歸求得切向、徑向和軸向的模型，如式(20)、式(21)和式(22)所示

Kt=82 488.861+6 057.801a-745.850n-

1 552.520F-1 040.347a2+2.161n2+28.378F2

(20)

Kr=164 592.802+3 842.936a-1 717.104n-

1 184.012F-722.357a2+4.932n2+22.163F2

(21)

Ka=-84 180.568+995.286a+782.904n+

667.366F-5.972a2-2.212n2-11.324F2

(22)

3 數(shù)值求解

3.1 基于隱式Adams方法的數(shù)值分析

基于隱式Adams方法的穩(wěn)定性預測判斷流程如圖9所示。

圖9 穩(wěn)定性預測流程圖Fig.9 Flow chart of stability prediction

將式(1)進行狀態(tài)空間變換，令初始時刻為t0，可以通過求解獲得

(23)

刀齒周期Tc可以分為兩個階段，自由振動時間間隔tz和強迫振動時間間隔Tc-tz。當?shù)毒咛幱谧杂烧駝訒r刻t，即t∈[t0,t0+tz]，則可以得到

x(t)=e(A(t-t0))x(t0)

(24)

加工過程中，當?shù)毒咛幱趶娖日駝訒r刻t時，即t∈[t0+tz,Tc]，將強迫振動的切削時間Tc-tz進行離散，平均分成m個時間間隔，則每個時間間隔可以表示為h=(Tc-tz)/m，則相應的離散點可以表示為

ti=t0+tz+(i-1)h,i=1,2,…,m+1

(25)

當t∈[ti,ti+1]時，代入式(23)可以得到

(26)

當i=1時，將式(25)代入式(24)可以得到狀態(tài)量x(t1)和時滯量x(tm+1-Tc)的關系，如式(27)所示

x(t1)=x(t0+tz)=e(Atz)x(t0)=

e(Atz)x(tm+1-Tc)

(27)

通過隱式Adams方法表示其他離散點，為了簡化表達式，將B(ti)表示為Bi，將x(ti)表示為xi，x(ti-Tc)表示為xi-T。

19e(Ah)Bi(xi-xi-T)-5e(2Ah)Bi-1(xi-1-xi-1-T)+

e(3Ah)Bi-2(xi-2-xi-2-T)]

(28)

簡化式(28)得到

Hi-2xi-2+Hi-1xi-1+(-e(Ah)+Hi)xi+(I+Hi+1)=
Hi-2xi-2-T+Hi-1xi-1-T+Hixi-T+Hi+1xi+1-T

(29)

聯(lián)立式(27)和式(29)可得

(30)

其中，

(31)

(32)

通過IAM法求得系統(tǒng)的狀態(tài)傳遞矩陣Φ為

Φ=M-1N

(33)

根據Floquet理論，當傳遞矩陣Φ的特征值等于1時，系統(tǒng)處于穩(wěn)定邊界；當傳遞矩陣Φ的特征值小于1時，系統(tǒng)處于漸進穩(wěn)定狀態(tài)；當傳遞矩陣Φ的特征值大于1時，系統(tǒng)會進入顫振狀態(tài)。

3.2 穩(wěn)定性分析

切向切削力系數(shù)、徑向切削力系數(shù)和軸向切削力系數(shù)可通過回歸模型式(20)、式(21)和式(22)計算得出。設置切向切削力系數(shù)Kt=10 499.441 N/m2，徑向切削力系數(shù)Kr=7 988.942 N/m2，軸向切削力系數(shù)Ka=-2 332.002 N/m2。將模態(tài)參數(shù)設置為表1中第一行參數(shù)，從安全性范圍考慮將主軸轉速范圍設置在常用的加工參數(shù)，主軸轉速n從150～230 r/min，并設置200個等間距的轉速；切削深度ap從0～4 mm，并設置200個等間距切削深度，穩(wěn)定性葉瓣圖構成了200×200的網格。設置離散數(shù)m為40，使用計算機軟件進行數(shù)值仿真結果如圖10所示。

圖10 螺桿轉子銑削穩(wěn)定性葉瓣圖Fig.10 Milling stability lobe diagram of screw rotor

如圖8所示，在穩(wěn)定性葉瓣圖中，曲線以下部分代表穩(wěn)定區(qū)域，即區(qū)域Ⅱ，在曲線以上的部分代表不穩(wěn)定區(qū)域，即區(qū)域Ⅰ。穩(wěn)定性葉瓣圖中存在多個葉瓣，所有葉瓣的最低點連成線所構成的穩(wěn)定區(qū)域為絕對穩(wěn)定區(qū)域?？梢詾橄乱徊降那邢髟囼炋峁┮罁?。

4 銑削試驗

為驗證本文提出的穩(wěn)定性預測方法，進行螺桿轉子銑削試驗。試驗所用銑床為LXK300G螺旋槽數(shù)控銑床，盤銑刀直徑為290 mm。被加工工件材料為45號鋼，直徑為50 mm，工件通過頂針和三爪卡盤固定，采用三個加速度傳感器多點布置，以保證采集振動信號的質量，如圖11所示。

圖11 試驗現(xiàn)場Fig.11 Test site

試驗的加工參數(shù)為：螺桿為左旋，導程為400 mm，螺旋角為18.511°，出于安全性考慮，主軸轉速設定為160 r/min，170 r/min，180 r/min，190 r/min，200 r/min，210 r/min，Zc方向的間歇位移量為常用的加工范圍，將其設置為2.0 mm，2.5 mm，3.0 mm，3.5 mm，4.0 mm，4.5 mm，5.0 mm。將間歇位移量經螺旋角轉換后得到切削深度如表6所示，將其反映到穩(wěn)定性葉瓣圖中，如圖12所示。

表6 銑削穩(wěn)定性試驗參數(shù)Tab.6 Parameters for milling stability experiment

圖12 銑削穩(wěn)定性驗證試驗參數(shù)Fig.12 Experimental parameters for milling stability verification

選取不同的主軸轉速和切削深度進行切削，利用DEWE Soft公司的SIRIUS-ACC數(shù)據采集器搜集振動信號，并根據統(tǒng)計理論獲得不穩(wěn)定點和穩(wěn)定點。穩(wěn)定點A的時域圖及頻譜圖如圖13所示。不穩(wěn)定點B的時域圖及頻譜圖如圖14所示。

圖13 穩(wěn)定點A的信號圖Fig.13 Signal diagram of stable point A

圖14 不穩(wěn)定點B的信號圖Fig.14 Signal diagram of unstable point B

由圖13和圖14可以看出，穩(wěn)定加工狀態(tài)和不穩(wěn)定加工狀態(tài)時域的振動信號圖中幅值相差近200倍，點A和點B的頻域特征均存在振動，該振動信號是由于測點布置、工件之間結合部和加工過程盤銑刀徑向切入量的變化導致的，除穩(wěn)定點A和不穩(wěn)定點B的共同特征外，在不穩(wěn)定點B的頻域信號中存相似明顯的顫振頻率，表明在螺桿轉子銑削過程中存在顫振現(xiàn)象，與解析結果吻合良好。對其他組測試數(shù)據分析可以得到類似結果。隱式Adams方法的計算機求解結果和試驗結果吻合良好，且適用于三自由度螺旋曲面銑削系統(tǒng)，能夠準確的預測切削穩(wěn)定性。在計算效率方面，與全離散法進行了對比，采用CPU為I5-4200U、內存為12 G的便攜式計算機進行計算，取離散度為40，全離散法耗時1 621.010 s，而隱式Adams方法耗時為682.128 s，耗時約為全離散方法的42%。

5 結 論

本文針對螺旋曲面銑削系統(tǒng)在加工過程中產生的振動現(xiàn)象，基于隱式Adams方法來預測螺旋曲面銑削加工過程的穩(wěn)定性，得到如下結論：

(1)通過分析螺旋曲面加工系統(tǒng)特性，建立了三自由度螺旋曲面銑削系統(tǒng)動力學模型。采用模態(tài)試驗，利用錘擊法獲取傳遞函數(shù)幅頻特性曲線，并應用導納圓擬合法提取穩(wěn)定性預測所需要的模態(tài)參數(shù)。

(2)提出基于隱式Adams的三自由度螺旋曲面銑削動力學系統(tǒng)穩(wěn)定性預測方法，并通過銑削試驗進行驗證。試驗結果表明，數(shù)值求解結果與試驗結果相吻合，即該方法適用于螺旋曲面銑削系統(tǒng)的穩(wěn)定性預測。