基于銑削系統動力學響應的球頭銑刀銑削表面形貌建模

董永亨， 李淑娟， 張倩， 李鵬陽， 李旗， 賈禎， 李言

(1.西安理工大學 機械與精密儀器工程學院, 陜西 西安 710048; 2.桂林電子科技大學 機電工程學院, 廣西 桂林 541004)

0 引言

零件的表面形貌對結合部接觸狀態、表面磨損、潤滑狀態、摩擦、振動等均有很大影響，快速預測加工表面的形貌對零件表面質量的提高至關重要。球頭銑刀因加工對象適應性強等優點而被廣泛應用于航空航天、模具和汽車等行業中重要零件的加工中。然而，對于高精度表面，大部分采用球頭銑刀銑削后再拋光的工藝路線，不但耗時耗力，而且影響零件表面質量的一致性。隨著高速加工技術、多軸控制技術和精密制造技術不斷進步和日趨完善，通過球頭銑刀精細銑削達到或接近零件的最終表面質量要求已經成為可能。然而，球頭銑刀的刀齒形狀復雜，切削時刀齒和工件之間的接觸點在不斷地發生變化，銑削過程中的幾何和物理因素對表面形貌形成的影響機理復雜。特別是在弱剛度銑削系統中，加工過程中由振動誘導的工藝系統動態位移對表面形貌的影響嚴重，直接決定著表面的使用效果。因此，探尋考慮工藝系統動態響應的球頭銑刀銑削表面形貌的建模與仿真方法具有重要的意義。

目前，銑削表面形貌建模與仿真的常用方法主要是Z-MAP法和數值計算法。

Z-MAP法是將工件進行離散，轉換成-網格節點，再根據切削加工中刀齒的位置求取相應每個網格節點對應的坐標，最后利用--矢量來表述已加工表面。Altintas等給出一種通用螺旋銑刀切削刃幾何模型，為復雜刀齒形狀銑刀加工表面形貌的建模奠定了基礎；閻兵等基于該方法建立了球頭銑刀的表面創成模型，具有良好的可擴展性；Quinsat等基于該方法分析了不同加工參數集合對零件表面五軸加工表面形貌的影響；Liu等基于該方法開發出一個球頭銑刀銑削表面特征和粗糙度的綜合仿真系統，在該系統里可以研究切削參數、刀具跳動、刀具磨損和刀軸傾斜對表面的影響；范思敏等預測了球頭銑刀銑削球面的表面形貌，為一般曲面加工表面形貌的建模奠定了基礎。傳統Z-MAP法雖然可以綜合考慮多種因素對表面形成的影響，但仿真精度和效率依賴于工件網格精度、刀齒離散精度和加工時間步長，當仿真精度要求高時，仿真效率就會降低，為此，一些研究者對傳統Z-MAP法進行了改進。趙厚偉等將具有高計算效率的實體建模法引入Z-MAP法中仿真預測了球頭銑刀加工表面形貌，但是，實體建模法中忽略進給方向殘留面積的缺點必然會影響結果的準確性。一些研究者用離散平面近似刀齒掃掠面，以提高Z-MAP計算效率。王仁偉等和Antoniadis等基于三角面和四邊面的刀齒微元掃掠面假設較高效地仿真了球頭銑刀銑削表面形貌，然而，用離散平面近似刀具掃掠面時，存在一定的近似誤差，并且采樣時間間隔越大，近似誤差越大。

數值計算法是在建立刀齒相對于工件的運動軌跡的基礎上，通過數值迭代的方法求解方程組，求出工件某坐標平面上選定點對應的輸入參數解，并將這些解回代方程組求出相應的軸向或法向高度，以最小的高度值作為坐標平面選定點高度，進而構建出整體加工表面形貌的一種建模方法。Gao等基于牛頓迭代法仿真了立銑刀周銑表面形貌，為球頭銑刀銑削表面形貌仿真提供了思路；譚剛等和Zhang等提出一種球頭銑刀多軸銑削表面形貌仿真的廣義迭代算法，無需離散刀齒，也無需對工件劃分網格，具有通用性好的特點；Chen等用牛頓- 辛普森公式求解了球頭銑刀間歇進給方向表面凹坑的高度坐標，但對進給方向表面凹坑的忽略容易在高速加工中引起較大的誤差；Mizugaki等提出一種球頭銑刀等高加工球面表面形貌的預測方法，用牛頓- 辛普森法求解出球面上對應位置的法向高度，但需借助CAD系統繪制出表面形貌；Arizmendi等基于切比雪夫展開式轉換提出了預測球頭銑刀加工表面形貌的模型，雖然無需設置初始點，但超越方程組的求解難度較大。數值法可以看作是運動學求逆解的過程，該方法在計算精度和計算效率上，都擺脫了對離散精度的依賴。然而，數值計算法主要從運動學角度計算最終成形表面的高度坐標，較難將諸如振動誘導下刀具和工件的動態位移等因素直接考慮進去，同時，數值計算的精度和效率與迭代精度的設置和初始點的選擇有很大關系。

文獻[20-21]將Z-MAP法和數值計算法有機地結合起來完成了球頭銑刀銑削表面形貌的建模與仿真，較好地兼顧了計算精度和效率，但未考慮工藝系統動力學響應的影響。梁鑫光等也綜合應用這兩種方法，提出考慮系統動態位移的球頭銑刀銑削表面形貌的仿真方法，然而，在計算表面殘留高度坐標時僅通過線性疊加的辦法考慮振動位移所引起的高度坐標變換，而忽略了在振動作用下作用于工件網格點的刀齒位置點變化。

本文以弱剛度銑削系統中球頭銑刀銑削表面形貌建模與仿真為研究對象，在建立刀具運動學模型和工藝系統動力學響應模型的基礎上，綜合應用Z-MAP法和泰勒公式，建立加工軌跡驅動下的球頭銑刀銑削表面形貌的幾何仿真模型，并記錄工件網格點的刀齒作用點和作用時間等信息。在此基礎上，以切削時間為紐帶，根據工藝系統的動態位移使用泰勒公式和線性插值法修正工件網格點的相關信息，求出考慮工藝系統動態位移的工件網格點高度，進而獲得考慮動力學相應的表面形貌。

1 刀齒幾何運動軌跡方程的建立

為了清楚地表達球頭銑刀銑削過程中刀齒上任意切削點的軌跡方程，建立如圖1所示的坐標系。圖1中：(簡稱{})為固連在工件上的全局坐標系；(簡稱{})為刀具瞬時進給坐標系，坐標軸矢量與進給速度方向平行且同向，為理想的被加工表面法線方向，指向實體外，為與的叉乘，也為和的叉乘；(簡稱{})為主軸隨動坐標系，該坐標系固連在機床主軸上，與主軸軸線重合，且使刀具遠離工件方向為正向，當與完全重合時，該坐標系的另外兩個坐標軸及其方向與{}的完全重合，但是，實際工況中當刀具姿態調整時，與之間存在夾角，該情況是使{}通過相對于和的旋轉實現主軸姿態的調整，進而實現刀具姿態的調整，從而獲得不同的銑削方式；(簡稱{})為刀具坐標系，坐標原點固連在刀具的球頭中心，與刀具的理論軸線重合，且與始終保持平行，二者之間的距離為，與基準刀齒(第1個刀齒)刃線在坐標平面上投影線起點的切線方向重合，該坐標系繞主軸坐標系中的坐標軸以角速度旋轉，坐標軸與之間的夾角為+(為主軸未開始旋轉的初始狀態下二者之間的夾角，為時刻主軸旋轉過的角度，=)；(簡稱{})為刀齒的局部坐標系，坐標原點也固連在刀具的球頭中心，坐標軸與完全一致，坐標軸與刀齒的刃線在坐標平面上投影線起點的切線方向重合。

圖1 球頭銑刀銑削運動的參考坐標系Fig.1 Reference coordinate system for ball-end milling

圖2 考慮刀具跳動的坐標系Fig.2 Coordinate system considering cutter’s jump

由于制造和裝夾誤差等因素的影響，刀具的中心軸線與主軸的中心軸線之間總存在偏心，如圖2所示，為點的軸向位置角(°)。假定刀具中心和主軸中心之間的偏心距離為(mm)，矢量相對于坐標軸的夾角為(°)，且規定繞坐標軸順時針旋轉方向為正；主軸順時針方向旋轉，其轉速為(r/min)，則角速度=π30(rad/s)，則(s)時刻旋轉過的角度=180π(°)，{}相對于{}的變換矩陣為

(1)

式中：=+，為初始狀態下與的初始夾角，順時針旋轉方向為正，為了簡化研究，本文設定=0°。

球頭銑刀銑削中往往使刀具相對于工件被加工表面傾斜，形成如圖3所示的側傾角(°)和前傾角(°)。先使{}繞旋轉角度′(°)，使′=arctan (tancos)，再使{}繞旋轉角度，且規定繞坐標軸矢量正向逆時針旋轉方向為正，反之為負。則{}相對于{}的變換矩陣

(2)

本文以單向直線進給銑削平面為研究對象，如圖3所示，此時

(3)

圖3 刀具的姿態調整及走刀Fig.3 Gesture adjustment and feed of the cutter

式中：(,)為首次進給時在{}中的起始位置；為刀具進給次數，=1,2,3…；為每齒進給量(mm/r)；為進給行距(mm)；為單次走刀長度(mm)；為毛坯高度(mm)；為刀具半徑(mm);為吃刀深度(mm)。

以定導程螺旋刃球頭銑刀為切削刀具，通過齊次坐標矩陣變換可得到加工過程中刀齒上任意點在工件坐標系{}下的軌跡方程：

(4)

式中：為該點所對應的螺旋滯后角(°)，=180tan(1-cos)π，為圓柱面上刀齒刃口曲線的螺旋角(°)；為刀齒與基準刀齒的夾角，=360(-1)(°)，為刀齒總數，=1～。

2 銑削系統動力學響應的建模

在弱剛度銑削系統中，由于剛度、阻尼和慣性的作用，刀具和工件會在銑削力的激勵下產生振動，進而誘導產生動態位移，該動態位移會疊加到刀齒幾何運動軌跡上，從而對加工表面形貌產生直接影響。

2.1 瞬時銑削力的求解

按照切削力機械建模法的思想，首先需要對刀齒微元進行微元切削力的建模，為此，將刀齒分割成等刀齒軸向位置角增量的諸多刀齒微元，如圖4所示，(,,)為與刀齒上離散點在時刻所在位置的連線在平面上的投影相對于坐標軸矢量順時針轉過的角度(°)，稱為實際徑向位置角。以刀齒離散點的特征信息來代表刀齒上(-1)～點之間的刀齒微元信息，刀齒上的微元在時刻所受的切削力可以分解為切向單元力切削力d(,,)、徑向單元力切削力d(,,)和軸向單元力切削力d(,,)，結合切削力的機械建模法，考慮到球頭銑刀刀齒的不同軸向位置角處切削狀態的不同，綜合考慮剪切和耕犁的共同作用，將其表示為

(5)

圖4 刀齒微元受力圖Fig.4 Force diagram of the discrete cutter teeth element

圖5 刀- 工切觸狀態的識別Fig.5 Identification of the cutter-workpiece engagement state

圖6 未變形切屑厚度求解示意圖Fig.6 Schematic diagram of the solution for undeformed chip thickness

將刀齒微元瞬時所受的切向力、徑向力和軸向力通過(6)式轉化至主軸隨動坐標系{}下，即

(6)

則刀具在時刻所受的瞬時切削力在主軸隨動坐標系{}下可表示為

(7)

式中：為刀齒微元總數。

通過齊次坐標變換原理得，刀具在時刻所受的瞬時切削力在工件坐標系{}下可表示為

(8)

2.2 工藝系統動態位移的求解

弱剛度銑削系統中工件和刀具振動的激勵包含靜態銑削力，而且包含動態銑削力，該動態銑削力主要由切削過程的再生效應產生。如圖7所示的弱剛度銑削系統，實際的切屑厚度不但包含名義厚度，還包含了振動誘導下刀具和工件的動態位移所引起的再生切屑厚度。

圖7 柔性刀具- 柔性工件動力學銑削系統Fig.7 Dynamic milling system of flexiblecutter-flexible workpiece

在銑削過程中，可以將刀具和工件簡化為如圖7所示的兩個垂直方向彈簧- 阻尼系統，則其動力學方程為

(9)

式中：、、分別為系統的模態質量矩陣、模態阻尼矩陣、模態剛度矩陣，

(10)

(11)

(12)

受刀齒結構和刀具跳動等因素的影響，球頭銑刀在銑削過程中的時滯量隨著刀齒、刀齒上的微元和切削時間的變化而變化，因此，穩定性分析中的時滯時間不再是刀齒的切削周期，而是主軸的旋轉周期，=60。根據文獻[26]，為了保證計算精度，在設定時間步長時需保證不能大于刀具- 工件系統最高振動模態周期的14倍，按照此要求將均勻地分為份，則時間步長單元

(13)

在Ding等提出的全離散法基礎上考慮切削過程的變時滯特性，求解(11)式就可以獲得球頭銑刀切削系統中工件和刀具的動態位移()、()、()、()。

2.3 考慮刀具和工件動態位移的刀齒運動學建模

(14)

(15)

因此，綜合考慮第1節所述的加工軌跡驅動下的刀具相對于工件的純幾何運動和振動誘導下刀具和工件的動態位移，則刀具坐標系{}相對于{}以及瞬時進給坐標系{}相對于{}的齊次變換矩陣分別如(16)式和(17)式所示：

(16)

(17)

結合(4)式，可以得出在考慮刀具和工件的動態位移情況下球頭銑刀加工過程中刀齒上任意點在工件坐標系{}下的軌跡方程：

(18)

因此，振動誘導下刀具和工件的總動態位移在工件坐標系{}下可以表達為

(19)

由(4)式、(18)式和(19)式，可得

(20)

3 表面形貌的仿真

圖8 考慮動力學響應的刀齒掃掠工件網格點高度求解示意圖Fig.8 Schematic diagram for calculating grid point height of workpiece swept by cutter teeth considering dynamic response

考慮振動誘導下刀具和工件動態位移的球頭銑刀銑削表面形貌的仿真過程如下：

1)參數設置。參數主要包括：

①毛坯尺寸及其Z-MAP模型表示，設置工件毛坯尺寸為××，用×個網格離散工件，保證圖8所示的工件在軸和軸方向的網格間距d和d滿足條件：max (d,d)≤15 min(,)；為了給工件的每個網格點賦予更多的信息，引入結構體數組[,](其中，=1，2，…，+1；=，2，…，+1)表示工件網格點的相關信息，同樣，用結構體數組的元素[,]·表示工件網格點(,)的高度坐標，仿真初始，為每一個結構體數組元素(,)·賦值，網格點(,)的其它信息如表1所示。

表1 工件網格關聯信息

②切削參數設置，主要有主軸轉速、每齒進給量、吃刀深度和切削行距等。

③刀具參數設置，設置球頭銑刀半徑、刀齒總數和螺旋角等，離散刀齒，設置刀齒上相鄰兩離散點的軸向位置角增量=180min (d,d)(π)，以保證刀齒微元在平面的有投影長度不超過d或d。

④刀具姿態設置，主要有側傾角和前傾角。

⑤表面形貌仿真時間步長的設置，為了保證計算精度，在離散切削時間時需要保證單位時間步長內刀齒離散點的掃掠軌跡曲線的弧長不大于工件網格間距，如(21)式所示：

(21)

式中：為最大軸向位置角，=arccos((-))。

⑥按(13)式設置用于切削系統動態位移計算的時間單元。

2)靜態銑削力計算。 在不考慮誘導產生銑削系統動態位移引起未變形切屑厚度變化的情況下，按照(5)式～(8)式計算球頭銑刀銑削過程中的靜態銑削力。

3)刀具和工件動態位移的計算。 按照22節求解在靜態力和再生效應所引起的動態力的共同作用下刀具和工件的動態位移，并按(19)式獲得工件坐標系{}下的任意切削時刻的總動態位移坐標((),(),())。

4)表面形貌的幾何仿真。 基于(4)式所建立的加工軌跡驅動下的刀齒運動軌跡方程，根據Z-MAP思想完成球頭銑刀銑削表面形貌的幾何仿真，為了在保證精度的前提下提高仿真效率，表面形貌的幾何仿真分兩步完成：表面形貌的預備幾何仿真和建立在其上的正式幾何仿真。

圖9 球頭銑刀銑削表面形貌幾何仿真示意圖Fig.9 Schematic diagram of geometric simulation of ball-end milling surface topography

(22)

上述用接近工件網格點的刀齒掃掠點高度坐標近似計算工件網格點高度坐標的方法雖然能夠提高仿真效率，但是存在一定的近似誤差，為此，在通過預備性幾何仿真獲得工件網格點相關的近似切削信息后，考慮到軸和軸方向坐標的近似誤差不會超過各自方向網格間距的一半，因此，可以基于刀齒運動方程的泰勒展開式實現誤差修正，進而獲得準確的表面形貌幾何仿真效果。由泰勒公式和(4)式，可得

(23)

解方程組可得Δ和Δ，結合(4)式和(23)式，工件網格點對應的當前時間段內刀齒掃掠面上的點的高度坐標為

(24)

用更新結構數組元素(,)·，并用(,)·+Δ和(,)·+Δ更新(,)·和(,)·，依次類推，按照此方法修正工件網格點的信息，通過各網格點修正后的(,)·值可以實現準確的表面形貌幾何仿真。

5)考慮振動誘導下工件和刀具動態位移的表面形貌的物理仿真。通過幾何仿真獲取了工件各網格點由幾何運動驅動的刀齒掃掠信息，但是，在振動誘導下刀具和工件動態位移的影響下，由幾何仿真所得到的刀齒掃掠信息會發生變化，從而導致工件網格點的掃掠高度也發生改變，因此，在考慮該動態位移的表面形貌物理仿真中需要重新計算。由于在實際加工中振動誘導的刀具和工件動態位移較小，因此，在獲取表面形貌的幾何仿真結果及刀具和工件動態位移結果基礎上，基于泰勒公式計算振動作用下加工表面的實際高度，進而獲得表面形貌的物理仿真結果。首先，需要求取每個工件網格點所對應的刀具和工件動態位移坐標，由步驟4獲取幾何運動驅動下每個工件網格點(,)的刀齒作用時間(,)·，通過其與動態位移計算時間單元的比值(,)·可以獲取該時間在動態位移時間歷程中的位置，由于用于表面形貌仿真計算的時間步長與的尺度不同，因此，(,)·在動態位移時間歷程中往往處于某一個區間，假定(,)·∈(,(+1))(=0,1,2,…)，而由(19)式可以獲得在某時間點=,(+1)處刀具和工件的動態位移坐標((),(),())和((+),(+),(+))，則可以由(25)式插值求得(,)·所對應的動態位移坐標值

(25)

由泰勒公式和(20)式，可得

(26)

按照(20)式將(26)式展開，可得

(27)

解方程組可得Δ和Δ，由此可得，在振動誘導的刀具和工件動態位移作用下，工件網格點(,)的刀具作用時間和作用刀齒的軸向位置角分別如(28)式和(29)式所示：

(,)·=(,)·+Δ

(28)

(,)·=(,)·+Δ

(29)

設=?(,)·」+1，則在振動誘導的刀具和工件動態位移作用下，工件網格點(,)的高度坐標更新為

(30)

4 實驗驗證

采用秦川機床廠生產的加工中心MV510作為實驗機床，該機床的最高主軸轉速為5 000 r/min，定位精度為0.018 mm，重復定位精度為0.008 mm；實驗刀具是整體長度為120 mm的Y330整體式硬質合金螺旋刃球頭銑刀，其直徑為10 mm，齒數為2，螺旋角為35°，彈性模量為210 GPa，刀具的裝夾懸伸長度為85 mm，安裝之后經測得徑向跳動量為0.028 mm；工件材料為航空鋁合金7050-T6，工件的長×寬×高為42 mm×11 mm×180 mm；為了在三軸機床上使刀具相對于工件的姿態調整，采用組合夾具的基礎板和圓柱支撐部件，按照正弦規的工作原理調整角度。

銑削系統的動態性能測試系統由力錘、加速度傳感器以及對應的數據采集設備和信號分析軟件LMS Test. Lab8B組成，如圖10所示，測試中采用力錘激勵，假定切削力集中在刀頭處，在刀頭上安裝加速度計，測量刀頭處的響應加速度，數據采集器將力錘激勵信號和加速度信號數字化采集到計算機上，然后由信號分析軟件LMS Test. Lab8B進行分析和處理，得到系統動態性能參數，結果如表2所示。

圖10 弱剛度銑削系統動態性能參數測試及切削實驗Fig.10 Dynamic property test and milling experiments for the weak-stiffness milling system

表2 刀具和工件的模態參數

在表3所示的切削條件下進行加工實驗，加工結束后用酒精擦拭加工表面，接著用德國產的徠卡激光共聚焦顯微鏡DCM-3D測量表面形貌，該設備的分辨率為0.1 nm，如圖11所示。

表3 實驗的切削參數

圖11 表面形貌測量裝置Fig.11 Surface topography measuring device

圖12 弱剛度銑削系統加工表面形貌的仿真和實驗結果對比Fig.12 Comparison of simulation and experimental results of surface topography machined by a weak-stiffness milling system

球頭銑刀銑削表面形貌的實測和仿真結果如圖12所示，二者的表征參數表面幅度算術平均偏差及輪廓幅度算術平均偏差對比如表4所示，無論從整體形狀，還是的對比結果上可以看出，仿真與實測結果都比較接近，二者在間歇進給方向的相應截面輪廓的吻合度也很高，主要是因為該方向截面輪廓是刀齒廓形的“復印”，實驗所用刀具的刀齒鋒利，且主軸轉速高，加工過程中不易產生積屑瘤，金屬在該方向的塑性流動小。然而，進給方向上的截面輪廓雖然在變化趨勢上一致，但有一定的差異，主要有以下5方面原因：

1)在強迫振動的沖擊作用下，弱剛度刀具的中心軸線較大幅地偏離主軸中心軸線，可能會出現副后角切削的情況，此時，刀齒后刀面與工件的摩擦增大，會使殘留凸體頂部的金屬略微向圖示右邊偏坍，使得峰谷截面曲線出現了圖示的偏態分布特點；

2)在沖擊力和摩擦力的共同作用下，接近刀齒刃帶的金屬會發生較大的塑性流動，流動方向沿著進給方向，接近塑性流動層的金屬會發生塑性變形和彈性變形，當刀齒切離工件時，又有彈性變形的恢復，從而使得實測輪廓的波動與理論仿真的有一定差異；

3)刀齒在工件表面形貌的谷底換向切出時可能會使金屬從基體撕裂，從而導致輪廓谷深較大；

4)由于更強的摩擦和擠壓作用，使得輪廓峰頂略向左偏移，造成了毛刺翻邊，使得粗糙度值增大；

5)力錘敲擊法測量模態參數存在一定誤差，會對工藝系統動態位移的計算精度產生一定的影響，進而對表面形貌預測產生直接影響。

表4 弱剛度銑削系統中實測和仿真表面Sa的對比

若改用文獻[21]所提出的方法實現本例的表面形貌幾何仿真環節，最終仿真所獲得表面形貌的為6.49 μm，在精度上與本文所提方法相當，而仿真計算所需的時間為695.4 s，相對于本文方法(551.2 s)效率顯低。

從整體上來講，仿真結果與實測結果有較好的一致性，可以用于弱剛度球頭銑刀銑削表面形貌的預測之中。

5 結論

1)按照齊次坐標變換原理建立了加工軌跡驅動下的球頭銑刀刀齒的運動方程，綜合Z-MAP法和泰勒公式，能夠實現球頭銑刀銑削表面形貌的預備性幾何仿真和正式性幾何仿真。

2)提出了考慮變時滯效應的弱剛度工藝系統動態位移的全離散求解方法，在修正刀具運動軌跡的基礎上，通過改進Z-MAP法完成考慮振動誘導下工件和刀具動態位移的表面形貌物理建模，實現了弱剛度球頭銑刀銑削系統銑削表面形貌的仿真。

3)在相同條件下，實驗和仿真表面形貌的表面幅度算術平均偏差的絕對誤差僅為0.02 μm，表明本文所所提出的考慮系統動力學響應的球頭銑刀銑削表面形貌建模方法是可靠的，可以為實際加工中參數的選擇和優化提供理論依據。