（西北工業(yè)大學(xué)機電學(xué)院，西安 710072）

金屬銑削加工過程中的影響因素很多，比如切削參數(shù)、刀具幾何參數(shù)以及工件、刀具材料等，這也使得新加工工藝的開發(fā)和切削條件的優(yōu)化變得異常復(fù)雜。傳統(tǒng)上，銑削過程中常采用試錯法選取合理的工藝參數(shù)以保證加工精度。然而，試錯法成本高、周期長，缺乏可操作性與量化分析，已經(jīng)遠(yuǎn)遠(yuǎn)不能滿足市場競爭的要求。開展基于切削物理機理的建模仿真技術(shù)研究，從理論上給出指導(dǎo)工藝參數(shù)選擇的方法，能有效解決上述問題。因此，關(guān)于切削過程機理及建模仿真研究工作也引起了國內(nèi)外學(xué)者的廣泛關(guān)注。加拿大英屬哥倫比亞大學(xué)Altintas教授團隊[1]圍繞該主題進(jìn)行了系統(tǒng)而深入的研究工作，取得的許多原創(chuàng)性成果得到了國內(nèi)外學(xué)者廣泛采用。國內(nèi)西北工業(yè)大學(xué)、北京航空航天大學(xué)、南京航空航天大學(xué)、山東大學(xué)、大連理工大學(xué)等高校也在相關(guān)領(lǐng)域進(jìn)行了深入研究，取得了許多積極有益的成果。本文將對銑削力建模、動力學(xué)建模與穩(wěn)定性分析及變形誤差預(yù)測與控制等工作進(jìn)行介紹，以期使讀者對銑削加工建模仿真技術(shù)有較系統(tǒng)的認(rèn)識。

銑削力建模方法

銑削力建模是所有銑削加工物理仿真的核心基礎(chǔ)。現(xiàn)有銑削力模型大致可以分為兩類：經(jīng)驗?zāi)Ｐ团c理論模型。經(jīng)驗?zāi)Ｐ褪峭ㄟ^大量切削試驗，對所測得的數(shù)據(jù)用數(shù)學(xué)方法進(jìn)行處理得到的。如Fuh[2]、Abou-El-Hossein[3]借助大量的試驗，利用統(tǒng)計學(xué)原理建立了切削力經(jīng)驗?zāi)Ｐ?。由于?jīng)驗?zāi)Ｐ托枰M(jìn)行大量試驗，且通用性較差，因此理論模型得到了更多的關(guān)注。理論模型是基于切削工藝幾何分析，經(jīng)過一系列力學(xué)推導(dǎo)得到的切削力函數(shù)。20世紀(jì)60年代，Koenigsberger和Sabberwal[4]建立了一種切削力理論模型，該模型認(rèn)為整個切削物理過程發(fā)生在剪切區(qū)域，工件材料在剪切區(qū)域發(fā)生剪切滑移產(chǎn)生的力即為整個切削力，因此切削力表達(dá)為與切屑載荷的比例函數(shù)，而比例系數(shù)則命名為切削力系數(shù)。之后，許多學(xué)者對該模型進(jìn)行了發(fā)展和延拓，形成了經(jīng)典的銑削力學(xué)模型（Mechanical Model）[1]，這里簡述如下:

圓柱坐標(biāo)系下的切削力可以表達(dá)為切削力系數(shù)（KT、KR、KA）與切削載荷的乘積，如下式：

式中，F(xiàn)T、FR和FA分別代表刀具切向、徑向和法向切削力；h代表瞬時未變形切屑厚度；b代表軸向切深。顯然，h、b取決于工藝參數(shù)，容易計算得到。因此，要計算切削力，關(guān)鍵是確定切削力系數(shù)。目前，關(guān)于切削力系數(shù)的確定主要有兩組方法，一種是直接標(biāo)定切削力系數(shù)，如Altintas[1]提出的用固定徑向、軸向切深與轉(zhuǎn)速，改變進(jìn)給率的方法標(biāo)定切削力系數(shù)，但該方法通用性不強，當(dāng)?shù)毒邘缀螀?shù)改變時，需要進(jìn)行重新標(biāo)定；另一種是通過正交與斜角切削幾何關(guān)系，建立切削力系數(shù)與刀具幾何和材料屬性參數(shù)表達(dá)式，通過計算得到切削力系數(shù)，即通過幾何關(guān)系將切削力系數(shù)轉(zhuǎn)化成剪切角（φn）、摩擦角（βn）、剪切應(yīng)力（τs）、刀具法向前角（αn）、刀具螺旋角刀（β）和切屑流動角（η）的函數(shù)：

上式中的刀具幾何參數(shù)可以測得，切屑流動角近似等于刀具螺旋角，因此只有剪切角（φn）、摩擦角（βn）、剪切應(yīng)力（τs）需要進(jìn)行標(biāo)定。在Altintas[1]提出的模型中，需要180組車削試驗來確定以上3個參數(shù)。本節(jié)簡要給出作者提出的快速銑削標(biāo)定算法[5]。

圖1 銑削過程示意圖Fig.1 Diagram of milling process

圖1給出了銑削過程示意圖[5]，在建模求解銑削力時，首先將銑刀沿軸向離散為許多微分單元，然后通過積分得到切削合力。圖1中i和j分別表示刀齒編號和微分單元編號。由斜角切削幾何關(guān)系有：

式中，F(xiàn)Rn,i(φ)為刀具第i個刀齒切削刃的法平面內(nèi)切向力，F(xiàn)Tn,i(φ)為刀具第i個刀齒的切削刃的法平面內(nèi)的進(jìn)給力，顯然，通過試驗測定相關(guān)切削力即可按式（3）標(biāo)定得到摩擦角。需要說明的是測力儀測定切削力一般為笛卡爾坐標(biāo)系下（即X、Y、Z方向）測力值，還需進(jìn)行坐標(biāo)變換得到FR和FT，詳細(xì)方法參考文獻(xiàn)[5]。剪切角的標(biāo)定基于最大剪應(yīng)力理論，即剪切速度方向和合力方向成45°角。結(jié)合斜角切削幾何關(guān)系，有：

上式相關(guān)參數(shù)定義同文獻(xiàn)[5]，解式（4）即可得到剪切角。確定了θi、θn和φi等角度，依據(jù)計算公式求解剪切應(yīng)力如下：

該模型由于標(biāo)定了求解切削力系數(shù)所需的物理參數(shù)，在刀具變化時，只需將相關(guān)刀具幾何參數(shù)代入式（2）計算即可，因此，該模型對于任何幾何刀具具有普遍適用性。

也有學(xué)者將上述銑削力模型稱為集成模型或一元模型。相應(yīng)地，將同時考慮了切削過程中的剪切效應(yīng)和犁切效應(yīng)（即將切削力表達(dá)為犁切力和剪切力兩項和的形式）稱為二元模型。在此基礎(chǔ)上，作者還提出了考慮底刃在切削過程中對切削力的作用的模型，即三元模型[6]。此外，作者還首次提出了銑削力可以表達(dá)為名義量與偏心引起的擾動量的疊加，并可以利用名義量來標(biāo)定切削力系數(shù)，擾動量標(biāo)定偏心參數(shù)。更多關(guān)于銑削力建模方面的工作，讀者可以參考文獻(xiàn)[7-10]。

銑削過程動力學(xué)建模

顫振是金屬切削加工過程中的一種較易發(fā)生、由切屑厚度動態(tài)再生引起的自激振動，也是導(dǎo)致切削失穩(wěn)的最主要因素。銑削過程顫振極易誘發(fā)刀具快速磨損、低精度波紋加工表面、噪音污染，還會嚴(yán)重影響機床壽命。因此，銑削顫振的動力學(xué)機理研究具有十分重要的學(xué)術(shù)和應(yīng)用價值。要對銑削工藝系統(tǒng)進(jìn)行動力學(xué)建模分析，首先要獲得工藝系統(tǒng)動力學(xué)特性，然后才能求解在切削力的激勵下，系統(tǒng)動力學(xué)響應(yīng)規(guī)律。因此，本節(jié)將闡述兩方面內(nèi)容，一是刀尖頻響建模預(yù)測方法，即工藝系統(tǒng)動力學(xué)特性的確定；另一方面，還將對顫振穩(wěn)定性預(yù)測方法進(jìn)行介紹，即系統(tǒng)動力學(xué)響應(yīng)規(guī)律。

1 刀尖頻響建模方法

目前，工藝系統(tǒng)動力學(xué)特性的確定主要有模態(tài)試驗和建模計算兩條主要途徑。由于實際生產(chǎn)中，一個零件加工過程往往選用多種刀具和刀柄組合，顯然試驗法更適用于實驗室，生產(chǎn)單位則更需要一套高效計算的方法。下文將介紹作者提出的一種適用于各種旋轉(zhuǎn)刀具的通用刀尖頻響建模思想[11]。

如圖2所示[11]，刀具-刀柄-主軸組合結(jié)構(gòu)被分解成為4個子結(jié)構(gòu)，最終刀尖點頻響是通過各子結(jié)構(gòu)動態(tài)響應(yīng)組合來獲得的。其中，SHS代表主軸-刀柄子結(jié)構(gòu)，STC代表刀具無刃桿部單元，F(xiàn)PTC代表刀具刃部單元，THJIC代表刀具-刀柄連接接觸單元。這里，刀具-刀柄被當(dāng)成零厚度阻尼彈性層。對于主軸-刀柄子結(jié)構(gòu)，首先將其分為NH個節(jié)點連接的（NH-1）個等長單元，則對該子結(jié)構(gòu)有：QH(ω)=HH(ω)FH-TH(ω)；其中，HH(ω)可以通過響應(yīng)耦合子結(jié)構(gòu)法求得。在獲得主軸-刀柄子結(jié)構(gòu)頻響后，依次對其他幾個單元進(jìn)行求解，最終組合得到刀尖點頻響結(jié)果。詳細(xì)求解信息，請參考文獻(xiàn)[11]。

2 顫振穩(wěn)定性預(yù)測方法

Tobias[12]等率先提出了應(yīng)用穩(wěn)定性葉瓣圖（Stability Lobe）來表征車削過程穩(wěn)定切削域，并為后續(xù)研究者廣泛采用以開展穩(wěn)定性分析工作。Altintas[1]指出連續(xù)的切削表面之間的相位移動引發(fā)指數(shù)增長的切屑厚度是誘發(fā)切屑再生顫振的根本原因，并提出了經(jīng)典的頻域求解穩(wěn)定域零階算法，同時對該算法進(jìn)行了許多后續(xù)完善。除頻域法外，許多學(xué)者還從時域研究了穩(wěn)定性預(yù)測問題，其中比較具有代表性的是匈牙利學(xué)者Insperger提出的半離散算法（Semi-Discretization Method）[13]。但是，該方法在提出時是基于理想工藝系統(tǒng)的，無法應(yīng)用于刀具偏心等引起的多延時問題。對此，Wan等[14]基于半離散法提出了一種通用多延時銑削穩(wěn)定性預(yù)測方法。隨后，又將該算法推廣到預(yù)測多模態(tài)銑削系統(tǒng)[15]。下文給出該方法的基本思想。

圖2 刀具-刀柄-主軸裝配建模原理圖Fig.2 Principle view of tool-holder-spindle assembly modeling

在銑削過程中，根據(jù)Altintas專著[1]中的定義，將關(guān)聯(lián)第j階模態(tài)在t時刻X、Y方向的刀具動態(tài)變形記作xj(t)和yj(t)。相應(yīng)地，第v個刀齒、第u個軸向單元徑向位移（或沿切屑厚度方向的位移可以表達(dá)為：

時間延時量τl對應(yīng)的動態(tài)切屑厚度為：

時間延時量的定義及確定方法參考文獻(xiàn)[14]。動態(tài)切屑厚度導(dǎo)致的動態(tài)切削力（FdX,j,uv(t),FdY,j,uv(t)），可以表達(dá)為：

其中，

進(jìn)而，基于模態(tài)疊加原理，構(gòu)建多模態(tài)多延時穩(wěn)定性控制方程，

其中，

式中，diag（*）代表對角矩陣運算，Nm表示主導(dǎo)模態(tài)數(shù)，ma,j、ca,j和ka,j分別表示a方向（a=X或Y）第j階模態(tài)下的模態(tài)質(zhì)量、阻尼系數(shù)和模態(tài)剛度。當(dāng)Nm等于1時，上述問題轉(zhuǎn)化為單模態(tài)問題，文獻(xiàn)[14]對于該問題求解進(jìn)行了詳細(xì)論述，并對刀具幾何、進(jìn)給等對穩(wěn)定性影響規(guī)律進(jìn)行了深入研究討論。對上式，采用半離散算法求解，即可獲得穩(wěn)定性葉瓣圖，更多信息可參考文獻(xiàn)[15]。由于該問題在求解時計算時間成本過高，作者還對該方法進(jìn)行了改進(jìn)，提出了一種最低包絡(luò)線的快速算法，并理論證明了該方法，同時對其適用性進(jìn)行討論，這里不再贅述。圖3給出了多延時系統(tǒng)預(yù)測的應(yīng)用實例[14]。

圖3 多延時銑削系統(tǒng)穩(wěn)定性結(jié)果Fig.3 Stability results for milling system with multiple delays

變形誤差預(yù)測與控制

1 變形誤差預(yù)測

加工過程中，由于切削力的作用，刀具和工件均產(chǎn)生彈性變形。走刀過后變形彈性恢復(fù)，致使部分材料殘留而未被切除，造成了零件表面的加工誤差。Kline指出，沿被加工面法向的分力是決定工件表面誤差的主要因素[16]。目前，對于無顫振的靜態(tài)周銑加工過程表面誤差的研究，主要有3類模型：（1）剛性模型，直接用名義切削量預(yù)測切削力，不考慮刀具及工件變形對加工過程的反饋影響，也不考慮材料去除引起的工件剛度變化的影響；（2）定剛度柔性模型，此模型相對剛性模型，在預(yù)測切削力時考慮刀具及工件變形的反饋影響，但不考慮材料去除引起的工件剛度變化的影響；（3）變剛度柔性模型，顧名思義，其預(yù)測切削力不僅考慮刀具/工件變形的反饋影響，而且還要考慮材料去除引起的工件剛度變化的影響。這里簡要介紹作者提出的一種基于商用軟件SAMCEF，同時考慮工件與刀具變形影響和材料去除引起剛度變化的算法[17]。

圖4 薄壁零件銑削加工示意圖Fig.4 Diagram of thin-walled workpiece milling process

如圖4所示[18]，使用平頭立銑刀銑削薄壁面AB，應(yīng)該切除陰影部分ABDCA，但由于切削力的作用，使刀具和工件均產(chǎn)生彈性變形，A、C和D點分別移到A'、C'和D'點，致使加工過程中的實際切削量（如徑向切削量）小于名義值，走刀過后變形彈性恢復(fù)，從而殘留部分材料未被切除，造成板厚方向的加工誤差。

變形誤差ep可表達(dá)為：

式中，δw為工件變形，δc為刀具變形。

刀具變形按懸臂梁理論計算，同時假設(shè)剛性主軸與夾頭中刀柄裝夾部分間為線彈性關(guān)系，即存在裝卡剛度。在某加工瞬態(tài)，刀齒節(jié)點（i，j）受作用于刀齒節(jié)點（k，m）的切削力而引起的靜態(tài)彈性變形可計算為：

式中，Ec表示刀具楊氏模量，Ic表示銑刀截面對其形心軸的慣性矩，kc表示刀具與夾頭間的裝卡剛度，表示銑刀裝卡后的懸臂長度。沿著刀具軸向?qū)λ械洱X變形求和即可得到所有刀齒節(jié)點受力在刀齒節(jié)點（i，j）引起的變形。

圖5 工件變形計算Fig.5 Calculation of workpiece deformation

工件加工過程中，影響加工精度的主要是工件彈性變形引起的，對工件變形的計算采用單位載荷法。如圖5所示[8]，若在A點、B點和C點分別加X方向、Y方向和Z方向的單位載荷，它們在D點X方向引起的變形分別為δ(AX,D,x)和δ(BY,D,X)和δ(CZ,D,X)，那么當(dāng)作用于A點X方向、B點Y方向和C點z方向的實際載荷為別為F(AX)、F(BY)、F(CZ) 時，D點X方向的變形大小為

刀具/工件變形的反饋影響按定剛度柔性模型中介紹的方法進(jìn)行考慮。對材料去除引起的工件剛度變化，借鑒結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計中的密度變量思想，考慮在恒定網(wǎng)格下修改被去除網(wǎng)格單元剛度的方法，實現(xiàn)對工件剛度的修正。

對于工件某網(wǎng)格單元i，如果有材料已被切除，采用下式對其剛度進(jìn)行修正：

式中表示工件單元i的名義剛度，ΔVi表示該單元i被切削后的剩余體積，Vi表示該單元i被切削前的實際體積。如果ΔVi=0，則取ηi=ε,取10-6=ε≤ηi是為了避免單元剛度矩陣的奇異性。更多詳細(xì)內(nèi)容，讀者可參閱文獻(xiàn)[17-19]。

圖6 雙參數(shù)誤差控制方法計算流程Fig.6 Optimization algorithm for surface form error control

2 變形誤差控制方法

當(dāng)預(yù)測變形誤差超過預(yù)設(shè)公差要求時，需對變形誤差進(jìn)行控制。鏡像補償法，即采用鏡像對稱原理對徑向切削深度進(jìn)行偏置補償?shù)姆椒ǎ瑥V泛用于變形誤差控制，其主要原理為：將名義徑向切削深度沿加工變形反方向偏移一個大小為變形量e的距離，可以實現(xiàn)減小變形誤差的目的。需要指出，這種方法雖然可以減小變形誤差,但不能保證補償后的變形誤差落在公差帶范圍內(nèi)。

為解決此問題，作者提出了一種雙參數(shù)優(yōu)化方法，即通過同時優(yōu)化單齒進(jìn)給量f與徑向切削深度ae，在維持生產(chǎn)率的同時實現(xiàn)最終誤差控制在公差帶范圍的目的，其數(shù)學(xué)模型可表達(dá)如下：

式中，emax和emin為給定單齒進(jìn)給量與徑向切削深度下變形誤差的最大值和最小值，cmax和cmin為給定變形誤差上下限，fmax和fmin為單齒進(jìn)給量上下限，ae,max和ae,min為徑向切削深度上下限。優(yōu)化求解式（14）的流程如圖6所示[20]，可以看出，求解的關(guān)鍵步驟之一是建立emax和emin與f和ae的關(guān)聯(lián)關(guān)系。

關(guān)于求解式（14）以及如何建立emax和emin與f和ae的關(guān)聯(lián)關(guān)系的詳細(xì)內(nèi)容，請參見文獻(xiàn)[20]。

結(jié)論

本文系統(tǒng)地介紹了銑削加工過程物理仿真的幾大關(guān)鍵技術(shù)問題，涵蓋了從切削力建模、變形誤差預(yù)測到工藝系統(tǒng)動力學(xué)建模與顫振穩(wěn)定性分析等多方面內(nèi)容，以期使讀者能夠較全面了解基于切削工藝力學(xué)機理的建模仿真技術(shù)。

[1]ALTINTAS Y.Manufacturing automation[M].Cambridgeshire: Cambridge University Press, 2012.

[2]FUH K-H, HWANG R-M.A predicted milling force model for high-speed end milling operation[J].International Journal of Machine Tools and Manufacture, 1997, 37(7):969-979.

[3]A B O U-E L-H O S S E I N K,KADIRGAMA K, HAMDI M, et al.Prediction of cutting force in end-milling operation of modified AISI P20 tool steel[J].Journal of Materials Processing Technology, 2007, 182(1): 241-247.

[4]KOENIGSBERGER F, SABBERWAL A J P.An investigation into the cutting force pulsations during milling operations[J].International Journal of Machine Tool Design and Research, 1961, 1: 15-33.

[5]WAN M, PAN W J, ZHANG W H, et al.A unified instantaneous cutting force model for flat end mills with variable geometries[J].Journal of Materials Processing Technology, 2014, 214(3):641-650.

[6]WAN M, LU M S, ZHANG W H, et al.A new ternary-mechanism model for the prediction of cutting forces in flat end milling[J].International Journal of Machine Tools and Manufacture, 2012, 57: 34-45.

[7]WAN M, ZHANG W H, DANG J W, et al.New procedures for calibration of instantaneous cutting force coefficients and cutter runout parameters in peripheral milling[J].International Journal of Machine Tools and Manufacture, 2009, 49:1144-1151.

[8]WAN M, ZHANG W H.Systematic study on cutting force modelling methods for peripheral milling[J].International Journal of Machine Tools and Manufacture, 2009, 49:424-432.

[9]WAN M, ZHANG W H, TAN G, et al.An in-depth analysis of the synchronization between the measured and predicted cutting forces for developing instantaneous milling force model[J].International Journal of Machine Tools and Manufacture, 2007,47: 2018-2030.

[10]WAN M, ZHANG W H, QIN G, et al.Efficient calibration of instantaneous cutting force coefficients and runout parameters for general end mills[J].International Journal of Machine Tools and Manufacture, 2007,47:1767-1776.

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

[12]TOBIAS S A.Machine tool vibration[M].New York: Blackie and Sons Ltd,1965.

[13]INSPERGER T, STEPAN G.Updated semi-discretization method for periodic delaydifferential equations with discrete delay[J].International Journal for Numerical Methods in Engineering, 2004, 61: 117-141.

[14]WAN M, ZHANG W H, DANG J W, et al.A unified stability prediction method for milling process with multiplt delays[J].International Journal of Machine Tools &Manufacture, 2010, 50(1): 29-41.

[15]WAN M, MA Y C, ZHANG W H,et al.Study on the construction mechanism of stability lobes in milling process with multiple modes[J].International Journal of Advanced Manufacturing Technology, 2015,79:589-603.

[16]KLINE W A, DEVOR R E,SHAREEF I A.The prediction of surface accuracy in end milling[J].Journal of Engineering for Industry-Transactions of the ASME, 1982,104: 272-278.

[17]WAN M, ZHANG W H, QIU K P,et al.Numerical prediction of static form errors in peripheral milling of thin-walled workpieces with irregular meshes[J].Transactions of the ASME Journal of Manufacturing Science and Engineering, 2005, 127: 13-22.

[18]萬敏.薄壁件周銑加工過程中表面靜態(tài)誤差預(yù)測關(guān)鍵技術(shù)研究[D].西安: 西北工業(yè)大學(xué), 2005.

WAN Min.Numerical prediction of static form errors in the peripheral milling of thinwalled workpiece[D].Xi’an: Northwestern Polytechnical University, 2005.

[19]WAN M, ZHANG W H.Efficient algorithms for calculations of static form errors in peripheral milling[J].Journal of Materials Processing Technology, 2006,171:156-165.

[20]WAN M, ZHANG W H , QIN G H, et al.Strategies for error prediction and error control in peripheral milling of thin-walled workpiece[J].International Journal of Machine Tools and Manufacture, 2008, 48: 1366-1374.