微細立銑刀動力學分析的Chebyshev譜方法

黃意新， 尹青峰， 趙　陽

(1.哈爾濱工業大學 航天學院，哈爾濱　150001； 2. 中國工程物理研究院 機械制造工藝研究所，四川 綿陽　621900)

?

微細立銑刀動力學分析的Chebyshev譜方法

黃意新1， 尹青峰2， 趙陽1

(1.哈爾濱工業大學 航天學院，哈爾濱150001； 2. 中國工程物理研究院 機械制造工藝研究所，四川 綿陽621900)

采用 Chebyshev 譜方法對微細立銑刀動力學特性進行研究。考慮刀具剪切變形與轉動慣量效應，建立刀具各段均勻或變截面 Timoshenko 梁動力學模型，利用 Chebyshev譜方法對動力學方程進行離散，通過動態子結構法綜合得到刀具整體動力學模型，編制了刀具動力學特性分析軟件;仿真得到了微細立銑刀固有頻率及模態振型并與分段階梯梁方法及有限元方法進行了比較，驗證了方法的可行性;研究了刀柄長度、刀頸半錐角及刀頭長徑比等刀具結構參數對刀具固有頻率的影響。算例的數值計算結果表明：當 Chebyshev 多項式階數>16時，前四階計算結果偏差收斂至 0.04% 內，計算精度較高;該方法可用于微細立銑刀參數優化設計。

微細立銑刀；切比雪夫譜方法；動力學特性；動態子結構法

微銑削是加工微小零件和高精密零件的一種新型加工技術[1]，是微細加工研究領域中，由電加工跨入非電加工、硅微工藝跨入非硅工藝、二維加工跨入三維加工的一項重要的先進制造技術[2]，在制造航空航天等領域所用的三維復雜微結構零件方面具有明顯優勢[3]。由于加工尺度的縮小，微細銑削刀具尺寸小、剛性弱，加工質量與穩定性對刀具性能敏感[4-5]，刀具結構參數將影響刀具動力學特性進而影響加工過程中的刀具動態響應[6]。

目前用于銑刀動力學分析的模型主要有：整體懸臂梁模型、分段階梯梁模型及有限元模型等[7]。張俊等[8]基于等質量、等截面積以及等剛度的原則將刀齒等效為均勻直徑梁，主要用于快速獲取銑刀頻響特性。方澤平等[9-11]應用分段思想建立微細銑刀懸伸部分的動力學模型，將刀具分為三部分，刀柄、刀頭部分為等截面圓柱梁，刀徑部分則簡化為等梯度變化的階梯梁，分析了刀具結構參數對刀具動力學特性的影響。YANG等[12-13]將銑刀分為：裝夾段、懸伸段及刀刃段三部分，采用 Euler-Bernoulli 梁和 Timoshenko梁理論建立了銑刀動力學模型用于刀尖頻響函數分析。MUSTAPHA等[14]應用廣義Hamilton變分原理推導了微細銑刀動力學模型，采用譜元法計算了刀具固有頻率隨刀具結構參數的變化規律。譜方法是繼差分法和有限元法之后又一種重要的數值方法，它用在整個區間都非零的連續函數的線性組合來逼近精確解，其精度可直接由級數展開式的項數來決定，具有指數收斂特性，常被視為具有“無窮階”收斂性，當真解足夠光滑時，采用譜方法可以得到很好的效果。譜方法可應用于微分方程數值求解[15]，梁動力學問題分析等[16]。從函數近似的角度看，譜方法可以分為Fourier方法，Chebyshev或Legendre方法。前者適用于周期性問題，后者適用于非周期性問題[17-18]。

為獲得微細立銑刀的精確動力學特性，基于Timoshenko梁理論，采用Chebyshev譜方法及動態子結構法建立微細銑刀的變截面Timoshenko梁動力學模型，分析刀具結構參數對刀具動力學特性的影響，為微細銑刀結構優化設計提供參考。

1　微細立銑刀結構動力學模型

微細立銑刀見圖1(a) ，包括刀柄、刀頸、刀頭三部分，并可簡化為圖1(b)的變截面懸臂梁模型。刀柄為直徑d1，長L1的均勻圓柱梁；刀頸為半錐角φ的錐形梁，其長為L2；刀頭為直徑d3，長L3的均勻圓柱梁。由于截面非均勻，應用動態子結構法，將刀具分為(1)、(2)、(3) 三個子結構，①、④為邊界，②、③為子結構交界面。首先獲得子結構運動方程，再根據邊界條件、子結構交界面的位移和力的協調條件，將離散的子結構運動方程組裝成整體結構運動方程并計算其動力學特性。

圖1　微細立銑刀及其幾何模型Fig.1　Micro-endmill and its geometric model

對于子結構 (2)，變截面Timoshenko梁運動方程為：

(1)

式中：κ為剪切系數;G為材料剪切模量;E為材料楊氏模量;橫截面積A=A(x)、慣性矩I=I(x)沿軸向發生變化。令：χ=ρ/κG，σ=E/κG，F=-f/κG，Ψ=-ψ/κG，式 (1) 可表示為：

(2)

對于子結構(1)、(3)，橫截面積A、慣性矩I為常數，運動方程式(2)可表示為：

(3)

2　Chebyshev 譜方法

2.1Chebyshev多項式理論

文獻 [19] 給出了Chebyshev譜方法理論的詳細論述。此處給出Chebyshev多項式逼近函數及其導數以及函數內積的方法，作為推導Timoshenko梁離散運動方程的基礎。第一類Chebyshev多項式是一類遞歸正交多項式，其第k項可表示為：

Tk(x)=cos(karccosx)

(4)

式中：k為非負整數，多項式滿足遞推關系：

Tk(x)=2xTk-1(x)-Tk-2(x)

(5)

式中：x∈[-1,1]。Chebyshev也可以定義在其它區間范圍，如l∈[0,L]，需對l進行變換x(l)=(2l-L)/L。函數y(x)的N-1階Chebyshev多項式逼近為：

(6)

式中:ak為多項式系數，組成N維系數向量a，采用Gauss-Lobattoe采樣點xk=cos((k-1)π/(N-1))，k=1,2,…,N，插值點函數值yN(xk)組成N維函數值向量y，則有：

y=ΓBa

(7)

式中：ΓB為函數值向量與多項式系數向量間的變換矩陣，在已知N個采樣點及其函數值時可求得系數向量a，進而可以獲得函數y(x)在任意點的函數值，由式(7)可得：

(8)

Chebyshev多項式導數：

T0′(x)=0,T1′(x)=T0(x),

T2′(x)=4T1(x)

(9)

(10)

(11)

上述各式中k為大于1的正整數，由式(6)、式(9)、式(10)、式(11)可得函數y(x)導數的N階Chebyshev多項式逼近：

(12)

多項式系數bk構成N維系數向量b:

b=Da

(13)

式中：D可由式(10)、式(11)、式(13)獲得:

(14)

綜合式 (8)、式(12)可得函數y(x)的N階導數值向量：

(15)

由式 (6)、式(7)、式 (15)可知，在已知函數y(x)N個采樣點上函數值的情況下，可以采用Chebyshev多項式逼近y(x)的各階導數值。

Chebyshev多項式積分：

(16)

式中：n為正整數，由式(6)、式(16)可得：

(17)

式中：v為含N維列向量，稱作Chebyshev積分向量，其元素可由式 (16) 得到。

函數內積采用Chebyshev多項式逼近可表示為：

(18)

V、Vh為N×NChebyshev內積矩陣。

由Chebyshev多項式理論可知，在已知函數N個采樣點函數值時，可以采用Chebyshev多項式逼近函數及其各階導數或積分值，計算其與其它函數的內積。

2.2Timoshenko梁離散運動方程

對于子結構 (2)，采用N階Chebyshev多項式逼近其運動方程式(2) 的解：y=y(x,t)，α=α(x,t):

(19)

由于有限項Chebyshev多項式逼近函數時存在截斷誤差，不能精確滿足微分方程式(2)，將產生殘差：

(20)

則方程式(2)的等效積分形式為：

(21)

式中：θ1(x)、θ1(x)為權函數，采用Galerkin加權余量法，θ1(x)、θ1(x)亦為Chebyshev多項式，設x為Gauss-Lobattoe節點向量，y、α分別為函數y=y(x,t)、α=α(x,t)某時刻在節點處的函數值向量，A、I分別為A(x)、I(x)的在節點處的函數值向量，由式(15)、式(18)可知：

y′=Q1y, α′=Q1α

(22)

(23)

(24)

則式(21)的離散形式為：

θ1T[(VAQ2+VA′Q1)y-

(25)

θ2T[σ(VIQ2+VI′Q1)α+

(26)

令：

(27)

式(25)、式(26)可寫成：

(28)

(29)

(30)

(31)

式中：0為N×N矩陣塊。

對于子結構(1)、結構(3)，橫截面積A、慣性矩I為常數，上述四式可簡化為：

(32)

(33)

(34)

(35)

(36)

式中:I為N階單位矩陣。

2.3邊界條件及系統綜合方程

根據式(28)～式(31)、式(32)～式(36)可分別得到3個子結構的未施加約束條件的離散運動方程，現綜合各方程得到系統運動方程，令：

qG=[y(1)α(1)y(2)α(2)y(3)α(3)]T

(37)

式中下標(i)代表子結構i，由式(28)、式(32)可得：

(38)

(39)

(40)

(41)

各子結構交界面位移與力的協調條件為：

y(1)(1)=y(2)(N)，y(2)(1)=y(3)(N)

(42)

α(1)(1)=α(2)(N)，α(2)(1)=α(3)(N)

(43)

α′(1)(1)=α′(2)(N)，α′(2)(1)=α′(3)(N)

(44)

y′(1)(1)-α(1)(1)=y′(2)(N)-α(2)(1)

(45)

y′(2)(1)-α(2)(1)=y′(3)(N)-α(3)(1)

(46)

邊界條件：y(0,t)=0，α(0,t)=0,α′(L，t)=0,y′(L,t)-α(L，t)=0,離散形式為：

y(1)(N)=0, α(1)(N)=0

(47)

(48)

式 (32)～式 (48) 等12個線性齊次邊界條件可表示為：

(49)

其離散形式可表示為：

βq=b

(50)

式中:β為12×6N邊界條件矩陣，其每行代表一個邊界條件，b為12維零向量，對于式 (42)～式(48) 中的邊界條件有：

(51)

式中:e1、eN為N維行向量，其第1個、第N個分量分別為1，其余均為0；對于任意滿足式 (50) 的解q可表示為：

(52)

(53)

(54)

(55)

3　算例

3.1固有頻率及模態

選取刀具結構參數見表1，材料密度ρ=14.45×103kg/m3，楊氏模量E=580 GPa，剪切模量G=242 GPa，截面剪切系數κ=0.9。得到刀具前四階固有頻率收斂性見圖2，當N等于20時，前四階固有頻率分別為：5 089.714 Hz、29 024.152 Hz、71 461.987 Hz、108 043.023 Hz，與文獻 [9] 中吻合，結果對比如表2所示。當N繼續增大，計算結果相對偏差收斂至<0.000 1%、0.001%、0.002%、0.04%。若只計算一、二階固有頻率，當N=10時，相對偏差便可收斂至<0.05%，映證了Chebyshev譜方法具有的快速收斂性與高精度特性，相較于分段階梯梁方法具有一定優勢。

表1　刀具結構參數

表2　刀具各階固有頻率比較

得到刀具前四階振型見圖3，圖3(a)為振動位移y前四階模態，圖3(b)為彎曲角α前四階模態，各圖中豎直線代表刀具各段交界面②、③。從圖中可知,刀具變形以刀頭部分最大，高階振型更為顯著，三階、四階位移模態中，從交界面③處相對位移迅速增大，二階、三階、四階彎曲模態中，從交界面③處相對彎曲角亦迅速增加，說明影響刀具變形量的主要因素是刀頭部分參數。

圖2　模型收斂性Fig.2　The convergence of the first four natural frequencies of the micro end-mill

圖3　刀具前四階模態振型，N = 20Fig.3　The first four mode shapes and natural frequencies of the micro end-mill when N =20

3.2結構參數影響分析

為合理設計微細立釬刀使之滿足高精度的動力學特性要求，需要分析各刀具參數對其固有頻率、模態振型的影響，基于所建立的刀具動力學模型，分析了刀柄長度、刀具半錐角、刀頭長徑比等參數對刀具動力學特性的影響。

圖 4為刀柄長度對固有頻率的影響，可知其它參數不變時，前四階固有頻率均隨刀柄懸長的增加而減小，以一階固有頻率為例，當刀柄長度L1=15 mm時，其值為8 259.901 Hz，而當刀柄長度增加至25 mm時，其值減小至3 440.102 Hz，減幅達58.352%。而在高速微細銑削中，由于刀徑減小，為獲得較大的切削速度就需要很高的主軸轉速，可達150 000 r/min，對于雙刃微徑銑刀，其切削力頻率可達5 000 Hz，因此當刀具一階固有頻率低于切削力頻率時極可能造成刀具共振，導致刀具、加工件損壞。因此在滿足裝卡條件的前提下，應盡可能縮短刀柄懸伸段長度。

圖4　刀柄長度對固有頻率的影響Fig.4　Change of natural frequencies with length of shank for micro-endmill

圖5為刀頸段半錐角對刀具固有頻率的影響，可知其它參數不變時，刀具前四階固有頻率均隨半錐角的增大而增大，增幅逐漸減小。以一階固有頻率為例，當半錐角從6°到10°時，其值由4 086.832 Hz增大到4 944.045 Hz，增加857.213 Hz，從10°到14°時，增加428.826 Hz至5 372.871 Hz。因此，適當增加刀頸半錐角有利于提高刀具固有頻率。

圖5　刀具半錐角對固有頻率的影響Fig.5　Change of natural frequencies with taper angle for a micro-endmill

圖6為刀頭長徑比對刀具固有頻率的影響，可知刀頭長徑比影響非常小。只有第四階固有頻率隨刀頭長徑比增加而減小，其它各階頻率基本保持不變。這說明在微細銑刀設計時，刀頭長徑比應更多關注如何增加刀頭靜、動剛度，不必過于關注其對刀具固有頻率的影響。

圖6　刀頭長徑比對固有頻率的影響Fig.6　Change of natural frequencies with ratio of length to diameter of tool for a micro-endmill

4　結　論

(1) 對于微細立銑刀動力學特性分析，Chebyshev譜方法具有良好收斂性和較高的計算精度，在選取較小Chebyshev多項式階數時，仍然獲得低階固有頻率、模態振型較高精度的解。相關結果與分段階梯梁方法及有限元方法吻和，證明了方法的可行性。

(2) 分析刀具前四階模態振型可知，刀柄、刀頸及刀頭中刀頭變形量最大，刀具振動位移、彎曲角度在各段交界面處迅速增大，振型階數越高越顯著。

(3) 刀具各結構參數中刀柄長度對固有頻率影響最大，前四階固有頻率隨刀柄長度增加而迅速減??；隨前刀頸半錐角的增大，前四階固有頻率逐漸增大，但增大幅度逐漸減??；刀頭長徑比對刀具固有頻率的影響很小。

[1] 肖才偉. 振動微銑削動力學建模及其分析[D]. 哈爾濱：哈爾濱工業大學, 2008.

[2] 李迎. 微銑削加工技術研究現狀及發展趨勢[J]. 電子機械工程, 2008, 24(6):26-32.

LI Ying. Recent adavances in micro milling and its future development trend[J]. Electro-Mechanical Engineering, 2008, 24(6): 26-32.

[3] 陳明君, 陳妮, 何寧, 等. 微銑削加工機理研究新進展[J]. 機械工程學報, 2014, 50(5):161-172.

CHEN Mingjun, CHEN Ni, HE Ning, et al. The research progress of micromilling in machining mechanism[J]. Journal of Mechanical Engineering, 2014, 50(5): 161-172.

[4] MIAO J C, CHEN G L, LAI X M, et al. Review of dynamic issues in micro-end-milling[J]. International Journal of Advanced Manufacturing Technology,2007,31(9):897-904.

[5] SHUNMUGAM M S. Machining challenges: macro to micro cutting[J]. Journal of the Institution of Engineers(India): Series C, 2015, 2(22):1-19.

[6] SHI Y, MAHR F, WAGNER U V, et al. Chatter frequencies of micromilling processes: influencing factors and online detection via piezoactuators[J]. International Journal of Machine Tools & Manufacture, 2012, 56(2):10-16.

[7] 萬敏, 張衛紅. 銑削過程中誤差預測與補償技術研究進展[J]. 航空學報, 2008, 29(5):1340-1349.

WAN Min, ZHANG Weihong. Overviews of technique research progress of form error prediction and error compensation in milling process[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(5): 1340-1349.

[8] 張俊, 黃保華, 趙萬華,等. 面向動態特性快速求解的銑刀等效建模方法[J]. 振動工程學報, 2013, 26(3):351-356.ZHANG Jun, HUANG Baohua, ZHAO Wanhua, et al. Equivalent modeling of endmills for rapid dynamics prediction[J]. Journal of Vibration Engineering, 2013, 26(3): 351-356.

[9] 方澤平, 王西彬, 劉志兵, 等. 基于Timoshenko梁理論的微細銑刀動力學建模[J]. 振動與沖擊, 2014, 33(23): 111-115.

FANG Zeping, WANG Xibin, LIU Zhibing, et al. Dynamic modeling of a micro-endmill based on Timoshenko beam theory[J]. Journal of Vibration and Shock, 2014, 33(23): 111-115.

[10] 方澤平. 微細立銑刀的設計基礎理論與刃磨技術研究[D]. 北京：北京理工大學, 2014.

[11] RAHNAMA R, SAJJADI M, PARK S S. Chatter suppression in micro end milling with process damping[J]. Journal of Materials Processing Technology,2009,209(17):5766-5776.

[12] YANG Y, WAN M, MA Y C, et al. An improved method for tool point dynamics analysis using a bi-distributed joint interface model[J]. International Journal of Mechanical Sciences, 2015, 105(2016): 239-252.

[13] YANG Y, ZHANG W H, MA Y C, et al. Generalized method for the analysis of bending, torsional and axial receptances oftool-holder-spindle assembly[J]. International Journal of Machine Tools and Manufacture,2015,99(2015): 48-67.

[14] MUSTAPHA K B, ZHONG Z W. A new modeling approach for the dynamics of a micro end mill in high-speed micro-cutting[J]. Journal of Vibration & Control, 2013, 19(6):901-923.

[15] 李雪慧. 截斷Chebyshev多項式的譜方法求解數值微分問題[D]. 蘭州：蘭州大學, 2010.

[16] YAGCI B, FILIZ S, ROMERO L L, et al. A spectral-Tchebychev technique for solving linear and nonlinear beam equations[J]. Journal of Sound & Vibration,2009,321(1):375-404.

[17] TREFETHEN L N. Spectral methods in MATLAB[M]. Philadelphia, Pennsylvania: Society for Industrial & Applied Mathematics , 2001.

[18] 王健平. 譜方法的基本問題與有限譜法[J]. 空氣動力學學報, 2001, 19(2):161-171.

WANG Jianping. Fundamental problems in spectral methods and finite spectral method[J]. Acta Aerodynamica Sinca, 2001, 19(2): 161-171.

[19] BOYO J P. Chebyshev and fourier spectral methods[M]. 2nd ed. New York: Dover Publications, 2000.

Dynamic analysis of micro-endmill by Chebyshev spectral method

HUANG Yixin1, YIN Qingfeng2, ZHAO Yang2

(1. School of Astronautics, Harbin Institute of Technology, Harbin 150001, China；2. Institute of Machinery Manufacturing Technology, China Academy of Engineering Physics, Mianyang 621900, China)

A Chebyshev spectral method for studying the dynamic characteristics of micro-endmill was presented. Considering the shear deformation and rotary inertia of micro-endmill, the uniform and tapered Timoshenko beam models were used in the dynamic analysis of the segments of endmill. The discrete forms of the dynamic models were gained by Chebyshev spectral method and synthesized by dynamic substructure technique to get the equation of motion of the whole system. A software package for solving the natural frequencies and modal shapes was developed. To validate the method, the impacts of structural parameters on dynamic characteristics were analysed. The simulation results were compared with those obtained by the stepped beam method and finite element method. It is shown that the calculated result converges into a range of 0.04% when the order of polynomials is larger than 16. The method is of good precision and can be applied to optimize structural parameters of micro-endmill.

micro-endmill; Chebyshev spectral method; dynamic characteristics; dynamic substructure technique

國家“973”計劃(2013CB733004)

2015-11-20修改稿收到日期：2016-02-04

黃意新 男，博士生，1987年生

趙陽 男，教授，博士生導師，1968年生

TH122；TG714

A DOI:10.13465/j.cnki.jvs.2016.14.005