基于Kriging模型的數(shù)控機床空間切削穩(wěn)定性*

鄧聰穎， 苗建國， 楊閃閃， 殷國富， 魏 博

(1.重慶郵電大學(xué)先進(jìn)制造工程學(xué)院 重慶，400065) (2.四川大學(xué)空天科學(xué)與工程學(xué)院 成都,610065)

引 言

切削過程中出現(xiàn)的顫振失穩(wěn)現(xiàn)象是惡化工件加工質(zhì)量和加工效率、加速機床磨損和精度喪失的重要原因[1-2]。通常情況下，一般的鋁合金等輕金屬的高速切削穩(wěn)定性主要受到刀具系統(tǒng)動力學(xué)特性的影響，而鑄鐵、鋼等材料在粗加工與半精加工時的切削穩(wěn)定性主要與機床主體結(jié)構(gòu)動力學(xué)特性相關(guān)[3]。在加工過程中，刀具點空間位置的改變使機床整體結(jié)構(gòu)隨之變化，形成動態(tài)的整機質(zhì)量矩陣、剛度矩陣和阻尼矩陣，引起機床動力學(xué)特性表征指標(biāo)變化，進(jìn)而使切削穩(wěn)定性在整個加工空間發(fā)生演變[4]。因此，研究機床切削穩(wěn)定性與加工位置的關(guān)聯(lián)關(guān)系，有利于提高切削穩(wěn)定性預(yù)測精度。

近年來，國內(nèi)外許多學(xué)者研究了機床結(jié)構(gòu)動力學(xué)特性與加工位置的關(guān)系[4-6]，但是機床用戶更關(guān)心面向加工過程的切削穩(wěn)定性。Law等[7]采用模態(tài)選擇準(zhǔn)則對機床子部件進(jìn)行簡化，耦合各子部件構(gòu)建機床整機動力學(xué)模型，通過更新結(jié)合面約束方程計算不同位置處刀尖點的頻響函數(shù)，進(jìn)而計算各位置的切削穩(wěn)定性。楊毅青等[8]基于實驗研究了機床沿不同進(jìn)給方向以及機床主軸處于不同位置的顫振穩(wěn)定域圖預(yù)測，但是針對每個加工位置進(jìn)行實驗或仿真計算，需要較大的工作量。文獻(xiàn)[3，9]結(jié)合實驗設(shè)計、有限元仿真以及響應(yīng)面模型，可計算機床任意位置的最小極限切削深度，但是在建立響應(yīng)面模型時未考慮樣本數(shù)據(jù)的空間相關(guān)性，也沒有進(jìn)一步探討加工位置的優(yōu)選方法。

筆者以機床加工點空間位置變化為研究基點，考慮切削穩(wěn)定性表征指標(biāo)的空間相關(guān)性，建立最小極限切削深度的Kriging預(yù)測模型，并采用改進(jìn)粒子群算法確定具有最小極限切削深度最大值的加工位置。在此基礎(chǔ)上，采用切削實驗確定機床易顫振模態(tài)，進(jìn)而基于能量分布理論確定機床薄弱結(jié)合部，通過提出結(jié)合部動剛度優(yōu)化配置方案進(jìn)一步提高最小極限切削深度值，以期更加精確地預(yù)測切削穩(wěn)定性空間分異特性，并擴大穩(wěn)定區(qū)域選擇范圍。

1 機床切削穩(wěn)定性空間信息預(yù)測建模

1.1 機床空間切削穩(wěn)定性概念

加工過程中切削力的動態(tài)變化量[10]可表示為

(1)

其中：ΔFx和ΔFy分別為x和y方向的切削力變化量；ap為軸向切削深度；Ktc為切向切削力系數(shù)；A為方向系數(shù)矩陣。

(2)

其中：Krt為徑向切削力系數(shù)和切向切削力系數(shù)的比值；φst和φex分別為切入角和切出角；φjl為當(dāng)前刀齒瞬時接觸角。

在頻域范圍內(nèi)分析式(1)可得

(3)

其中：F為切削力矩陣；G(iωc)為工件-刀具接觸區(qū)域的傳遞函數(shù)矩陣。

當(dāng)工件系統(tǒng)剛性較大時，通常采用刀尖點頻率響應(yīng)函數(shù)矩陣表示。式(3)的特征方程及特征值Λ為

(4)

根據(jù)式(4)，極限切削深度aplim及主軸轉(zhuǎn)速n為

(5)

采用圖1所示的顫振穩(wěn)定域葉瓣圖描述式(5)中切削深度和主軸轉(zhuǎn)速的關(guān)系，直線表示最小極限切削深度，其以下為絕對穩(wěn)定區(qū)，在任意轉(zhuǎn)速下選擇該區(qū)域內(nèi)的切削深度加工都不會發(fā)生顫振。針對一般較硬材料的粗加工、半精加工，切削過程中選擇的主軸轉(zhuǎn)速較低，對應(yīng)的條件穩(wěn)定區(qū)相對狹窄，因此通常選擇在最小極限切削深度以下進(jìn)行加工[3]。傳遞函數(shù)矩陣G(iωc)是繪制切削顫振穩(wěn)定域圖的關(guān)鍵。研究表明，G(iωc)受加工位置的影響，因此針對較硬材料在粗加工、半精加工階段的切削穩(wěn)定性準(zhǔn)確預(yù)測問題，以最小極限切削深度aplimmin在加工空間的動態(tài)變化為研究對象，提出更為全面的廣義最小極限切削深度函數(shù)

aplimmin=F(x,y,z)

(6)

其中：x，y，z分別為各向運動軸的位移。

基于式(6)，進(jìn)一步將aplimmin與加工位置的關(guān)聯(lián)關(guān)系顯性化，即可求解aplimmin在加工空間的演化規(guī)律。

圖1 顫振穩(wěn)定域葉瓣圖Fig.1 Chatter stability chart

1.2 廣義最小極限切削深度Kriging預(yù)測模型

Kriging方法是基于已知樣本信息的動態(tài)構(gòu)造，從變量的空間相關(guān)性和變異性角度，構(gòu)建對象問題的近似函數(shù)關(guān)系式，對非樣本點信息進(jìn)行無偏、最優(yōu)估計[6,11]。該方法包括線性回歸部分和非參數(shù)部分，通常Kriging近似模型具有如下表達(dá)

(7)

協(xié)方差矩陣可表示為

cov(z(xl),z(xj))=σ2R(xl,xj)

(8)

R(xl，xj)為樣本空間中任意兩個樣本點xl和xj的空間相關(guān)函數(shù)，目前廣泛采用高斯函數(shù)和指數(shù)函數(shù)表達(dá)，它決定了模擬的精度，依次具有如下形式

(9)

其中：參數(shù)θk為待定參數(shù)；xlk，xjk為樣本點xl，xj的kth分量。

根據(jù)實驗樣本點及其對應(yīng)的響應(yīng)值，形成一個相關(guān)矩陣

(10)

采用加權(quán)最小二乘法，得到Kriging模型的回歸系數(shù)矩陣

(11)

基于樣本點xs的響應(yīng)值Y的線性加權(quán)差值來預(yù)估待測點x的響應(yīng)值

(12)

近似值與真實值之間的誤差為

(13)

將式(9)、式(10)的函數(shù)代入式(13)，可得

基于無偏差預(yù)測要求

(15)

式(15)的均方差為

(16)

結(jié)合上式與方差最小條件、拉格朗日乘子法

(17)

綜上，Kriging模型可最終表達(dá)為

(18)

其中：γ*=R-1(YT-Fsβ*)。

為檢驗Kriging模型的預(yù)測精度，采用設(shè)計空間任意樣本點的相對誤差(RE)與平均相對誤差(ERR)兩個標(biāo)準(zhǔn)[14]來評價

(19)

如果RE與ERR的值越接近0，則表明Kriging模型的精度越高。

構(gòu)建廣義最小極限切削深度Kriging預(yù)測模型的主要流程如下：

1) 采用實驗設(shè)計方法離散化機床工作空間，確定樣本點的空間位置坐標(biāo)(xi,yi,zi)，并通過實驗獲取各樣本點位置對應(yīng)的頻率響應(yīng)函數(shù)G(iωc)i；

2) 繪制各樣本點位置對應(yīng)的顫振穩(wěn)定域葉瓣圖，進(jìn)而確定各樣本點位置的最小極限切削深度；

3) 基于式(7)～式(18)構(gòu)建廣義最小極限切削深度的Kriging預(yù)測模型

aplimmin(x,y,z)=fT(x,y,z)β+rT(x,y,z)γ*

(20)

4) 求解新樣本點處的最小極限切削深度，結(jié)合式(19)、式(20)檢驗步驟3中Kriging預(yù)測模型的精度，根據(jù)驗證結(jié)果調(diào)整模型，提高擬合精度。

2 基于Kriging預(yù)測模型的廣義空間切削穩(wěn)定性優(yōu)化問題描述

本研究中廣義空間切削穩(wěn)定性優(yōu)化問題是在給定切削寬度和進(jìn)給量的基礎(chǔ)上，基于式(20)建立的最小極限切削深度Kriging預(yù)測模型，尋求具有較高最小極限切削深度值的加工位置，并在此基礎(chǔ)上進(jìn)一步提高最小極限切削深度值，擴大切削參數(shù)選擇范圍。

2.1 改進(jìn)粒子群優(yōu)化算法

粒子群算法作為一種機器學(xué)習(xí)的群體智能優(yōu)化算法，能有效實現(xiàn)非線性和多峰值復(fù)雜目標(biāo)的優(yōu)化[12]。在D維空間中，假設(shè)n個粒子組成一個群落，每個粒子的特性通過位置Xi=(xi1,xi2,...,xiD)、速度Vi=(vi1,vi2,...,viD)以及適應(yīng)值fi來表征。每一個粒子位置對應(yīng)該優(yōu)化的一個潛在解，速度表征粒子當(dāng)前的飛行方向和距離，適應(yīng)度值則采用適應(yīng)度函數(shù)計算，表示粒子在當(dāng)前位置的優(yōu)劣性。每個粒子的位置更新取決于個體極值Pbset和群體極值Gbest的變化，同時更新適應(yīng)度值，通過分析新粒子、個體極值和群體極值的適應(yīng)度值，進(jìn)一步更新個體極值與群體極值。位置與速度按下式更新

(21)

其中：ω為慣性權(quán)重；i為第i個粒子，i=1,2,…,n；d為維數(shù)，d=1,2,…,D；c1和c2表示速度加權(quán)因子為0或正的常數(shù)；r1和r2為[0,1]之間的隨機數(shù)。

慣性權(quán)重ω的設(shè)置對粒子群算法的求解極為重要，為更好地滿足實際要求，迭代求解過程中基于粒子不同的特性，動態(tài)調(diào)整慣性權(quán)重ω。假設(shè)全部粒子平均適應(yīng)度以favg表征，fmax為本輪迭代中的最大適應(yīng)度值，fhavg代表適應(yīng)度值大于favg的粒子的平均適應(yīng)度值，則迭代過程中按如下方式動態(tài)調(diào)整ω[12]：

1) 當(dāng)fi≤favg時，優(yōu)先考慮提高該類粒子的全局搜索能力，應(yīng)設(shè)置較大的ω

(22)

2) 當(dāng)favg

(23)

3) 當(dāng)fhavg

(24)

式(22)～式(24)中，右側(cè)ω值遵循線性遞減規(guī)律，k為當(dāng)前迭代次數(shù)，kmax為初始化最大迭代次數(shù)

(25)

為避免粒子群算法出現(xiàn)早熟收斂，采用變異算子的形式對粒子群算法中某些變量以一定概率重新初始化

(26)

以加工點各向位移x，y，z為變量，最小極限切削深度的最大值為目標(biāo)，結(jié)合式(23)～(26)可建立刀具點加工位置尋優(yōu)模型，并按圖2所示流程求解。

圖2 改進(jìn)粒子群算法流程Fig.2 Improved particle swarm optimization algorithm

2.2 基于能量分布的切削穩(wěn)定性優(yōu)化

在2.1節(jié)確定的機床最優(yōu)加工點位置附近進(jìn)行切削時，為擴大切削參數(shù)選擇范圍，提高切削效率，從機床結(jié)構(gòu)動力學(xué)特性和切削穩(wěn)定性關(guān)聯(lián)關(guān)系角度進(jìn)行切削穩(wěn)定性優(yōu)化?？紤]機床振動問題60%以上來源于結(jié)合部，且結(jié)合部動剛度特性研究比阻尼特性研究更易實現(xiàn)，筆者結(jié)合切削實驗和能量分布理論，確定機床易顫振模態(tài)及對應(yīng)的薄弱結(jié)合部，通過優(yōu)化薄弱結(jié)合部動剛度來提高機床加工點動剛度，進(jìn)而提高最小極限切削深度。

首先以2.1節(jié)中求解的最小極限切削深度的最大值為基準(zhǔn)，線性遞增作為切削深度并進(jìn)行切削實驗，采用麥克風(fēng)拾取加工過程的聲音信號，根據(jù)聲音信號的頻譜分析確定機床易顫振模態(tài)。在此基礎(chǔ)上，采用能量分布理論分析機床在易顫振階模態(tài)振動時彈性能在整機中的分布情況。假設(shè)系統(tǒng)由N個模塊構(gòu)成，其中第m個模塊在i階模態(tài)振動時，系統(tǒng)與模塊m的彈性能[13]可表達(dá)為

(27)

其中：VAi，Vmi分別為系統(tǒng)在第i階模態(tài)振動時系統(tǒng)與模塊m的彈性能；Ai為i階模態(tài)下模塊m所有自由度的幅值向量；Ki為模塊m的剛度矩陣。

模塊m在系統(tǒng)中的彈性能分布率Ri定義為

(28)

機床在易顫振階模態(tài)振動時，如果計算的結(jié)合部Ri值越高，則表明該結(jié)合部剛度越低，需要對其進(jìn)行優(yōu)化設(shè)計。

2.3 廣義空間切削穩(wěn)定性優(yōu)化的實現(xiàn)

圖3為廣義空間切削穩(wěn)定性優(yōu)化流程，首先，基于構(gòu)造的樣本信息在Matlab環(huán)境中建立最小極限切削深度Kriging模型，并結(jié)合改進(jìn)粒子群算法求解最優(yōu)加工位置及最小極限切削深度；其次，采用切削實驗和彈性能分布率確定機床易顫振模態(tài)下的薄弱結(jié)合部，提出其動剛度優(yōu)化方案，進(jìn)而計算整機模態(tài)和切削穩(wěn)定性，獲取各薄弱結(jié)合部動剛度最優(yōu)配置以及優(yōu)化后的最小極限切削深度值。

圖3 廣義空間切削穩(wěn)定性優(yōu)化的實現(xiàn)Fig.3 Optimization of cutting stability in generalized space

3 基于Kriging預(yù)測模型的廣義空間切削穩(wěn)定性優(yōu)化實例

將提出的廣義空間切削穩(wěn)定性優(yōu)化方法應(yīng)用于普瑞斯數(shù)控機床廠的一臺三軸立式加工中心，如圖4所示。針對該機床耗時最多的某鑄件的粗加工、半精加工階段，采用的切削深度僅為2.42 mm。為提高加工效率，結(jié)合理論分析與圖3，對該工序的最小極限切削深度進(jìn)行優(yōu)化。

3.1 廣義最小極限切削深度Kriging模型構(gòu)建

工作臺、床鞍以及主軸箱的行程依次為x-400 mm，y-400 mm和z-550 mm。筆者結(jié)合中心復(fù)合表面設(shè)計和Box-Behnken實驗設(shè)計法規(guī)劃樣本點空間位置，如表1和圖4所示。各樣本點位置對應(yīng)的刀尖點頻率響應(yīng)函數(shù)(frequency response functions,簡稱FRF)由錘擊法模態(tài)實驗獲取，進(jìn)而通過式(1)～(5)計算對應(yīng)的最小極限切削深度aplimmin。

采用圖4中實心標(biāo)記的27個空間點作為位置變量輸入，關(guān)聯(lián)的aplimmin作為輸出，選擇隨機函數(shù)z(x,y,z)協(xié)方差計算模型為高斯函數(shù)，獲取正則化θ參數(shù)(θx=0.825,θy=0.625,θz=0.825)，完成Kriging預(yù)測模型構(gòu)建。為驗證該模型的準(zhǔn)確性，首先，采用該模型預(yù)測圖4中空心標(biāo)記的6個加工位置的aplimmin；其次，在該6個位置進(jìn)行切削實驗，刀具為直徑20 mm、螺旋角45°的4齒硬質(zhì)合金立銑刀，工件材料為45鋼，主軸轉(zhuǎn)速為3kr/min，徑向切削寬度為14 mm，進(jìn)給速度為0.06 mm/z。通過磁力座將PCB公司測量范圍為4 Hz～20 kHz的探頭型傳聲器固定在工作臺上，拾取切削過程的聲音信號并進(jìn)行頻譜分析，確定顫振及對應(yīng)的aplimmin，其采樣頻率為6 kHz。根據(jù)式(19)獲取由表2給出的樣本點實驗值和預(yù)測值的相對誤差(RE)與平均相對誤差(ERR)，其中RE最大值為3.12%，ERR最大值為2.03%，可見建立的Kriging預(yù)測模型擬合精度較高。

/mmxyz/mm/mmxyz/mm/mmxyz/mm13805303802.991020202002.7919380530201.9323802753801.8711202752003.2220380275202.673380203801.9812205302002.632138020203.324200203801.64132005302002.312220020203.0952002753801.73142002752002.1923200275201.8462005303802.6115200202002.6824200530201.537205303803.0716380202002.812520530201.978202753802.09173802752003.172620275202.64920203801.82183805302002.79272020203.46

表2 精度檢驗樣本點相對誤差

圖5描述了最小極限切削深度aplimmin在加工空間的演化規(guī)律，由圖可以看出，aplimmin隨加工位置變化而改變。圖5(b)中，y-z軸聯(lián)動時，aplimmin最低為1.47 mm，最高為3.10 mm，變化率達(dá)111.7%；圖5(e)中，x-z軸聯(lián)動時，aplimmin最低為2.15 mm，最高為3.27 mm，變化率達(dá)52.0%；圖5(h)中，x-y軸聯(lián)動時，aplimmin最低為1.74 mm，最高為3.30 mm，變化率達(dá)89.4%。分析圖5(a～c)發(fā)現(xiàn)，當(dāng)y-z軸坐標(biāo)確定時，工作臺x向位置的變化對aplimmin的影響相對較小。上述研究表明，在規(guī)劃機床切削工藝時，應(yīng)盡量使機床處于具有較高切削深度且其變化率較小的加工區(qū)域，并且應(yīng)優(yōu)先使用極限切削深度變化較小的運動方式來保障機床加工的平穩(wěn)性。

圖5 最小極限切削深度在加工空間的演化規(guī)律Fig.5 Evolution law of minimum limit cutting depth in machining space

3.2 廣義最小極限切削深度動態(tài)優(yōu)化

采用改進(jìn)粒子群算法獲取使最小極限切削深度值最大化的加工位置，首先設(shè)置基本參數(shù)如下：初始化種群規(guī)模為100，以加工點各向位移x，y，z為變量，種群中每個粒子的維數(shù)為3，粒子的位置和初始速度隨機給定，更新過程按式(22)～(25)進(jìn)行，根據(jù)式(20)建立的Kriging預(yù)測模型計算粒子的適應(yīng)度值，每次更新后采用式(26)進(jìn)行變異操作，并重新初始化粒子。按照圖2的求解流程，計算的加工空間aplimmin最大值為3.47 mm，對應(yīng)的最優(yōu)加工位置為30.77,14.12,36.78 mm。

在最優(yōu)加工位置附近進(jìn)行切削實驗，通過切削中聲音信號的頻譜分析確定顫振頻率分布范圍為406～443 Hz(如圖6所示的切削深度為4.00 mm時聲音信號的頻譜)。由于顫振常常出現(xiàn)在系統(tǒng)某階固有頻率附近，可根據(jù)機床整機模態(tài)信息判斷機床易顫振模態(tài)由第10階模態(tài)(423.7 Hz)引起，對應(yīng)模態(tài)振型表現(xiàn)為主軸-立柱系統(tǒng)繞y軸的擺動。因此，以該系統(tǒng)中z向主軸軸承結(jié)合部、z向滾珠絲杠結(jié)合部、z向?qū)к壔瑝K結(jié)合部以及床身-立柱結(jié)合部為研究對象，根據(jù)式(27)、式(28)計算出各結(jié)合部彈性能分布依次為26.04%，6.58%，32.80%和3.91%，從而判定主軸軸承結(jié)合部、導(dǎo)軌滑塊結(jié)合部為較薄弱結(jié)合部，并以該類結(jié)合部動剛度為優(yōu)化變量。

圖6 切削過程中聲音信號頻譜圖Fig.6 Spectrogram of sound signal in cutting process

主軸與直線導(dǎo)軌系統(tǒng)均為外購件，并且其待優(yōu)化結(jié)合部均為輕預(yù)載，因此根據(jù)產(chǎn)品技術(shù)手冊和課題組結(jié)合部動力學(xué)特性辨識成果，基于表3中的中預(yù)載和高預(yù)載情況下的結(jié)合部動剛度，按照表4的優(yōu)化方案依次將其寫入處于最優(yōu)工位的機床整機有限元模型中，以進(jìn)行諧響應(yīng)分析獲取各組刀尖頻響函數(shù)矩陣，采用式(1)～(5)計算各aplimmin。由表4看出，aplimmin由初始的2.42 mm提高到4.86 mm，增幅為100.8%，驗證了筆者提出的廣義空間最小極限切削深度優(yōu)化方法的有效性。

表3 薄弱結(jié)合部剛度數(shù)據(jù)

/(N·μm-1)k1k2/mm1110550.07673.952110550.09334.373157706.57674.694157706.59334.86

4 結(jié) 論

1) 基于Kriging方法理論建立最小極限切削深度的空間信息預(yù)測模型，研究其在加工空間的演化規(guī)律，融合了機床運動部件位姿的隨機組合對切削穩(wěn)定性的影響，有利于在完整加工空間中準(zhǔn)確描述切削穩(wěn)定性，提高了切削穩(wěn)定性預(yù)測精度。

2) 結(jié)合最小極限切削深度Kriging預(yù)測模型與改進(jìn)粒子群算法，可有效確定使最小極限切削深度最大化的加工位置，并在此基礎(chǔ)上引入切削實驗和彈性能分布理論，優(yōu)化機床在易顫振模態(tài)下薄弱結(jié)合部的動剛度配置，進(jìn)一步提高了最小極限切削深度，也為機床廠家外購件的選擇提供了技術(shù)支持。

3) 采用提出的機床空間切削穩(wěn)定性研究方法，可進(jìn)一步探討軸向極限切削深度、主軸轉(zhuǎn)速、切削寬度等切削穩(wěn)定性關(guān)聯(lián)指標(biāo)與加工空間的關(guān)聯(lián)關(guān)系，為更全面和更準(zhǔn)確地進(jìn)行切削加工參數(shù)選擇和工藝規(guī)劃提供理論依據(jù)和數(shù)據(jù)支持。