梁志強 石貴紅 杜宇超 葉玉玲 籍永建 陳司晨仇天陽 劉志兵 周天豐 王西彬

1.北京理工大學(xué)先進加工技術(shù)國防重點學(xué)科實驗室,北京,100081 2.北京神工科技有限公司,北京,100013 3.北京信息科技大學(xué)機電系統(tǒng)測控北京市重點實驗室,北京,100192

0 引言

工業(yè)機器人由于其靈活性、智能性、經(jīng)濟性等優(yōu)勢而被廣泛應(yīng)用于大型復(fù)雜結(jié)構(gòu)件的加工。然而，由于機器人開鏈式多自由度的結(jié)構(gòu)特點，導(dǎo)致其剛性與機床相比相差將近兩個數(shù)量級[1]，銑削過程極易產(chǎn)生顫振，嚴重降低了加工表面質(zhì)量與精度。理論推導(dǎo)機器人銑削穩(wěn)定性葉瓣圖是對機器人銑削穩(wěn)定性控制的一種有效方式[2]，其關(guān)鍵是對機器人銑削系統(tǒng)刀尖頻響函數(shù)的準確獲取[3-5]。由于機器人刀尖頻響函數(shù)的位姿依賴性，如何在變位姿下準確獲取機器人刀尖頻響函數(shù)，成為機器人銑削領(lǐng)域研究的熱點和難點。

目前用于獲取加工系統(tǒng)頻響函數(shù)的方法主要有模態(tài)實驗法、子結(jié)構(gòu)響應(yīng)耦合法(receptance coupling substructure analysis, RCSA)[5-6]與基于實驗的有限元分析法[7-8]。2000年，SCHMITZ等[9]推導(dǎo)出了裝配體頻響與各子結(jié)構(gòu)頻響的函數(shù)關(guān)系，首次提出了子結(jié)構(gòu)響應(yīng)耦合方法。隨后的研究者將該方法用于機床頻響計算、結(jié)合部建模與刀尖頻響預(yù)測。李宇庭[3]運用子結(jié)構(gòu)耦合方法結(jié)合模態(tài)實驗建立了機器人刀尖頻響預(yù)測模型，并通過實驗驗證了其準確性。不過加工系統(tǒng)子結(jié)構(gòu)之間的連接狀態(tài)模型和連接參數(shù)的辨識一直是該方法的難題，辨識方法還處于研究階段，距離工程應(yīng)用還有一定差距[10]，此外，該方法中的刀柄、刀具以及主軸端的頻響多依賴有限元仿真與實驗?zāi)B(tài)分析，預(yù)測模型建立的工作量較大。2013 年，加拿大的LAW等[11]利用有限元方法，提出了降階模型機床位置依賴的多體動力學(xué)建模方法，提高了機床刀尖頻響預(yù)測效率。MOUSAVI等[12-13]基于矩陣結(jié)構(gòu)分析(matrix structural analysis，MSA)方法的思想，利用ANSYS有限元模態(tài)仿真分析對MATLAB數(shù)值模型進行參數(shù)初步標定，再通過模態(tài)實驗對未知參數(shù)進行最后標定。有限元法雖適用于對不同結(jié)構(gòu)進行模態(tài)分析，但對于大型加工系統(tǒng)而言，隨著網(wǎng)格數(shù)增加，計算效率顯著降低，因此預(yù)測精確性與計算效率不可兼具。此外，國內(nèi)外學(xué)者也提出了其他方法用于機器人銑削系統(tǒng)動力學(xué)分析。CHEN等[4]利用已知一定空間范圍內(nèi)機器人銑削位姿的實驗?zāi)B(tài)參數(shù)，提出逆距離加權(quán)法，計算出任意位姿的頻響函數(shù)，并驗證了該方法在一定加工范圍內(nèi)的準確性。MOHAMMADI等[14]與PAN等[15]在研究機器人銑削中軸向振動對加工穩(wěn)定性的影響，以及機器人銑削中的模態(tài)耦合顫振產(chǎn)生機理與抑制方法時，將機器人本體簡化為質(zhì)量-彈簧-阻尼模型，分別計算模型中的質(zhì)量、剛度與阻尼矩陣?？傮w而言，對加工系統(tǒng)頻響預(yù)測方法的研究普遍存在效率與準確性無法兼具的問題。

對于機器人銑削系統(tǒng)的頻響預(yù)測，大都將末端主軸系統(tǒng)結(jié)合面考慮為剛性連接[3-5，12-15]，ALTINTAS等[16]指出，將主軸-刀柄結(jié)合面看作柔性，系統(tǒng)的頻響函數(shù)會受到影響，因此，在建模過程中不應(yīng)簡單地將主軸系統(tǒng)整體簡化為剛性連接。趙萬華等[17]提出一種可以考慮刀具-夾套、夾套-刀柄以及刀柄-主軸結(jié)合面接觸特性的主軸轉(zhuǎn)子系統(tǒng)動力學(xué)解析建模方法，解決了傳統(tǒng)主軸轉(zhuǎn)子系統(tǒng)動力學(xué)模型中的結(jié)合面接觸剛度通過實驗識別的方法獲取時通用性差的問題。但目前在機器人銑削系統(tǒng)的頻響預(yù)測研究中，很少對結(jié)合面接觸剛度的因素進行考慮。

綜上，本文提出了一種考慮主軸-刀柄結(jié)合面接觸剛度的機器人銑削系統(tǒng)刀尖頻響預(yù)測方法。將機器人本體動力學(xué)模型簡化為質(zhì)量-彈簧-阻尼模型，并結(jié)合歐拉-拉格朗日法、關(guān)節(jié)剛度辨識與比例阻尼模型依次確定了機器人的質(zhì)量、剛度和阻尼矩陣。在此基礎(chǔ)上，將末端執(zhí)行器的質(zhì)量與幾何參數(shù)補償?shù)奖倔w動力學(xué)模型中，并建立了主軸-刀柄結(jié)合面接觸剛度解析模型，基于有限元主副自由度理論對本體模型進行二次補償，從而構(gòu)建了機器人銑削系統(tǒng)刀尖頻響預(yù)測模型。

1 機器人銑削系統(tǒng)頻響預(yù)測模型建立

1.1 機器人本體動力學(xué)建模

建立機器人銑削系統(tǒng)刀尖頻響預(yù)測模型首先需要對機器人本體動力學(xué)進行建模。機器人銑削加工中進給速度遠低于切削速度，因此本體動力學(xué)方程中的慣性力(即科里奧利力和離心力)影響可以忽略；同時忽略非線性慣性力的影響，考慮機器人關(guān)節(jié)柔性，將本體動力學(xué)模型簡化為圖1所示的質(zhì)量-彈簧-阻尼線性系統(tǒng)[14]。同時假設(shè)工件相對于機器人系統(tǒng)為剛性，因此可忽略加工中的工件變形。在以上簡化的基礎(chǔ)上，機器人銑削加工系統(tǒng)的動力學(xué)模型可由微分方程表示：

圖1 機器人動力學(xué)模型簡化示意圖Fig.1 Simplified diagram of robot dynamics model

(1)

式中，M、C、K分別為機器人基坐標系中的質(zhì)量、阻尼和剛度矩陣；x(t)為位移矢量；F(t)為切削力矢量。

1.1.1機器人質(zhì)量矩陣計算

機器人本體動力學(xué)模型中的質(zhì)量矩陣M通過將關(guān)節(jié)坐標系下的質(zhì)量矩陣利用雅可比矩陣變換到基坐標系下得到：

M=J-TMqJ-1

(2)

(3)

關(guān)節(jié)坐標系下的質(zhì)量矩陣Mq通過KUKA KR210工業(yè)機器人慣性參數(shù)結(jié)合機器人改進D-H模型參數(shù)，利用歐拉-拉格朗日法推導(dǎo)得出。圖2所示為實驗用機器人的CAD模型，示出了該機器人通過模型導(dǎo)出的連桿2的質(zhì)量、初始位姿下的慣性張量和基坐標系中的質(zhì)心坐標。通過機器人的工作空間圖提取桿長數(shù)據(jù)，建立了機器人的改進D-H模型，其D-H參數(shù)如表1所示。各個連桿的慣性參數(shù)與D-H參數(shù)作為數(shù)值模型的初始輸入?yún)?shù)。

圖2 KUKA KR210工業(yè)機器人CAD模型與連桿慣性參數(shù)提取Fig.2 CAD model of KUKA KR210 industrial robot and inertia parameter extraction of links

表1 KUKA KR210 型機器人的改進D-H參數(shù)表Tab.1 Modified D-H parameter table of KUKA KR210 robot

根據(jù)歐拉-拉格朗日方程[18]，關(guān)節(jié)空間中的質(zhì)量矩陣如下所示：

(4)

式中，mi為連桿i的質(zhì)量；Jvi(q)和Jωi(q)分別為連桿i雅可比矩陣的上半部和下半部；Ri(q)為連桿i的定向變換矩陣；Ici為連桿i質(zhì)心的慣性張量矩陣。

1.1.2機器人剛度矩陣計算

機器人剛度矩陣通過將關(guān)節(jié)剛度矩陣利用雅可比矩陣變換到機器人基坐標系下得到：

K=J-TKqJ-1

(5)

式中，Kq為關(guān)節(jié)坐標系下的對角剛度矩陣。

機器人各關(guān)節(jié)剛度值通過載荷實驗辨識獲得。在機器人法蘭邊緣一點施加載荷，對比加載與空載狀態(tài)下的位移偏差，改變載荷質(zhì)量與施力方向重復(fù)多次，根據(jù)下式：

Fo=J-TKqJ-1D

(6)

式中，F(xiàn)o為外力；D為由于外力產(chǎn)生的位移變量。

利用最小二乘法擬合出機器人各關(guān)節(jié)剛度值。由于力矩值較小，此處只考慮測量點的位移變量。

測量對象為KUKA KR210工業(yè)機器人，采用Leica AT960/930型號的激光跟蹤儀進行位移測量，實驗原理和實驗現(xiàn)場如圖3所示。

(a)實驗原理

(b)實驗現(xiàn)場圖3 機器人關(guān)節(jié)剛度測試實驗Fig.3 Robot joint stiffness test experiment

首先將激光跟蹤儀的坐標系與機器人基坐標系重合，這樣激光跟蹤儀即可測量機器人基坐標系中某點坐標值。隨后用繩子和滑輪裝置在機器人末端分別懸掛10 kg、20 kg、30 kg的砝碼；分別將靶球安裝在圖3b中A、B位置，利用激光跟蹤儀測量向量AB的值作為力矢量F的方向(繩子方向)，將砝碼所施加力按照向量AB的方向余弦進行分解，即可得到基坐標系下x、y、z方向的分力。再用激光跟蹤儀測得機器人末端位置空載與加載的位移變量，即可擬合關(guān)節(jié)剛度值[19]。更換滑輪C1到C2改變載荷方向重復(fù)實驗。

擬合結(jié)果如下：

Kq=diag(3.3479×106，3.1318×106，4.2895×106，

2.4787×105，6.4292×105，1.8847×105) N·m/rad

(7)

1.1.3機器人阻尼矩陣計算

阻尼矩陣采用比例阻尼模型計算[20]，可將其簡化為已獲得的質(zhì)量矩陣和剛度矩陣的線性組合：

C=αK+βM

(8)

其中，α、β由刀尖頻響實驗標定得到。

1.2 主軸-刀柄結(jié)合面接觸剛度建模

1.2.1模型建立

由于主軸系統(tǒng)靠近刀尖，對刀尖動力學(xué)影響較大，將主軸-刀柄結(jié)合面剛性處理會對頻響函數(shù)預(yù)測精度造成誤差，因此建立主軸-刀柄結(jié)合面剛度模型。

在主軸-刀柄結(jié)合面剛度建模中，同時考慮法向、切向與扭轉(zhuǎn)接觸剛度，采用均布的彈簧單元模擬主軸-刀柄結(jié)合面剛度特性[17]，如圖4a所示。將主軸與刀柄的接觸區(qū)域沿軸向劃分為N個單元，每個單元簡化為以中徑為直徑的梁，每個單元的主軸面與刀柄面通過相對應(yīng)節(jié)點處產(chǎn)生的耦合力與力矩相互作用[21]，如圖4b所示。

(a)結(jié)合面整體彈簧單元分布

(b)結(jié)合面單元彈簧分布圖4 主軸-刀柄結(jié)合面接觸剛度模型Fig.4 Contact stiffness model of spindle-toolholder interface

當(dāng)兩個對應(yīng)節(jié)點產(chǎn)生相對運動時，其耦合力與力矩可由下式表示：

(9)

式中，ui1、vi1、θi1分別為部件A第i個單元沿x、y和繞z方向的位移量；ui2、vi2、θi2分別為部件B第i個單元沿x、y和繞z方向的位移量；kix、kiy、kiθ分別為結(jié)合面第i個單元沿x、y和繞z方向的接觸剛度；cix、ciy、ciθ分別為結(jié)合面第i個單元沿x、y和繞z方向的阻尼；Fx、Fy、Mz分別為結(jié)合面第i個單元沿x、y和繞z方向產(chǎn)生的耦合力與力矩。

相應(yīng)的每對節(jié)點的接觸剛度矩陣表示如下：

(10)

將每對節(jié)點的接觸剛度矩陣kij按照有限元法組裝即可得到整個結(jié)合面的總剛度矩陣KJ。

式(10)矩陣中接觸剛度kix、kiy和kiθ需要通過對主軸-刀柄結(jié)合錐面的法向、切向與角接觸剛度進行計算進而變換到機器人的基坐標系下獲得。

圖5為刀柄的受力示意圖，靜態(tài)下，主軸拉刀力使結(jié)合面產(chǎn)生接觸剛度。

圖5 刀柄結(jié)合錐面受力示意圖Fig.5 Force diagram of toolholder joint

平衡方程為

(11)

式中，F(xiàn)為主軸拉刀力；Fn為結(jié)合錐面所受法向正壓力；Fμ為結(jié)合錐面所受切向摩擦力；θ為錐角；μ為錐面摩擦因數(shù)。

刀柄接觸面的法向正壓力為

(12)

主軸-刀柄結(jié)合面接觸面積為

(13)

式中，Di為刀柄結(jié)合錐面第i個單元的大徑尺寸；di為刀柄結(jié)合錐面第i個單元的小徑尺寸；Li為刀柄結(jié)合錐面第i個單元的軸向長度。

主軸-刀柄結(jié)合面第i個單元所受壓力為

(14)

根據(jù)吉村允孝的單位面積法[22]，第i個單元的法向接觸剛度為

(15)

式中，α′、β′為結(jié)合面基礎(chǔ)特性參數(shù)，由實驗擬合得到。

結(jié)合面的切向剛度系數(shù)一般為法向剛度系數(shù)的1/3～2/3[22]，本文取每個單元的切向接觸剛度為其法向接觸剛度的1/2。

主軸-刀柄結(jié)合錐面第i個單元的角接觸剛度為

(16)

(17)

(18)

式中，kit為主軸-刀柄結(jié)合錐面的切向接觸剛度；li1與li2為幾何相關(guān)系數(shù)；Dim為主軸-刀柄結(jié)合面第i個單元中徑。

角接觸剛度的詳細計算原理見文獻[17]。

1.2.2結(jié)合面剛度矩陣的降階變換

有限元法將連續(xù)體(無限多自由度的系統(tǒng))離散為有限的N階自由度系統(tǒng)[23]，而本模型只考慮較低階模態(tài)，為最終簡化為二自由度系統(tǒng)，需要將總剛度矩陣KJ的維數(shù)進一步減小。利用有限元方法中的主副自由度理論，按照加速動力縮聚法，將矩陣左上角一定維數(shù)的元素定為主自由度，其他部分作為副自由度[24]，進行矩陣變換。

將主自由度記為δa，副自由度記為δb，總剛度矩陣按照所選取的主副自由度個數(shù)分塊表示：

(19)

只考慮主副自由度間的彈性聯(lián)系，則主副自由度間滿足以下約束方程：

Kbaδa+Kbbδb=0

(20)

由此可得

(21)

(22)

(23)

式中，T為變換方程，是一個長方形矩陣。

將總剛度矩陣KJ變換為所需的低階剛度矩陣：

(24)

最后，綜合考慮機器人本體與末端執(zhí)行機構(gòu)的動力學(xué)以及主軸-刀柄結(jié)合面接觸剛度，將式(24)代入式(1)中，加工系統(tǒng)的整體動力學(xué)方程如下：

(25)

2 模型標定與實驗驗證

2.1 模型標定

由于機器人執(zhí)行末端結(jié)構(gòu)較復(fù)雜，幾何、材料等相關(guān)數(shù)據(jù)無法準確獲取，另外阻尼系數(shù)未知，因此需要對這些參數(shù)進行標定。在將末端執(zhí)行器與機器人第六軸看作剛性連接的情況下，選擇預(yù)測模型中第六軸的質(zhì)量m6與桿長d6作為調(diào)整參數(shù)來標定仿真頻響函數(shù)的固有頻率，選擇阻尼系數(shù)α、β作為調(diào)整參數(shù)來標定仿真頻響函數(shù)的幅值。本文利用模態(tài)實驗獲得不同位姿下的多組系統(tǒng)刀尖頻響函數(shù)，基于建立的考慮主軸-刀柄結(jié)合面接觸剛度的機器人銑削系統(tǒng)動力學(xué)模型，通過使實驗與數(shù)值仿真得到的頻響函數(shù)的固有頻率及其峰值差值最小化，來識別數(shù)值模型中以上的未知參數(shù)[13]。根據(jù)下面的公式通過數(shù)值模型計算得到頻響函數(shù)：

(26)

式中，ωnk、ηk分別為數(shù)值模型的第k階模態(tài)的固有頻率和阻尼比；pik、pjk為數(shù)值模型的第k階模態(tài)的振型矢量；Fj(ω)和Xi(ω)分別為在數(shù)值模型自由度j上的激勵力和自由度i上的位移量。

將1.1.1節(jié)與1.1.3節(jié)得到的質(zhì)量矩陣與阻尼矩陣轉(zhuǎn)化到模態(tài)空間后，可得到第k階模態(tài)質(zhì)量mk和模態(tài)阻尼ck。

加工系統(tǒng)安裝的是FIACHER SD5084型號的主軸和TRIBOS-S HSK-E32型號的刀柄。實驗設(shè)備使用安裝有YD-5T型石英傳感器的力錘敲擊刀尖產(chǎn)生激勵信號，使用INV9822型加速度傳感器采集刀尖處的響應(yīng)信號，加速度傳感器靈敏度為10.355 mV/(m·s-2)，采用INV3062T型4通道數(shù)據(jù)采集卡，采樣頻率設(shè)置為5120Hz。最后利用此套設(shè)備配備的模態(tài)分析模塊計算得到實驗頻響函數(shù)。機器人銑削系統(tǒng)與信號采集設(shè)備和實驗示意圖見圖6。

(a)系統(tǒng)與采集設(shè)備

(b)實驗示意圖圖6 實驗設(shè)備與測試示意圖Fig.6 Experimental equipment and test diagram

標定過程如圖7所示。首先進行頻率標定，調(diào)整m6、d6兩參數(shù)使仿真頻響的固有頻率與實驗重合，如圖7b所示；再進行阻尼標定，調(diào)整α、β兩參數(shù)使仿真頻響的峰值與實驗重合，如圖7c所示；此過程改變位姿重復(fù)進行多次，優(yōu)化所選參數(shù)的標定值。

(a)未標定

(b)頻率標定

(c)阻尼標定圖7 數(shù)值模型標定過程示意圖Fig.7 Numerical model calibration process diagram

2.2 實驗驗證

為驗證機器人銑削系統(tǒng)刀尖頻響函數(shù)預(yù)測模型的準確性，開展刀尖頻響函數(shù)測試實驗。實驗隨機選取機器人的5個位姿,如表2所示，采用錘擊法對刀尖頻響函數(shù)進行測試，并與上述數(shù)值模型計算得到的頻響結(jié)果對比，以驗證本文提出的數(shù)值模型的準確性。

表2 所選機器人位姿的各關(guān)節(jié)角Tab.2 Joint angles of the selected robot postures

采用式(26)的頻響預(yù)測數(shù)值模型對以上所選的機器人位姿進行仿真頻響計算。在主軸-刀柄結(jié)合面接觸剛度矩陣計算中，將刀柄與主軸的結(jié)合面沿軸向劃分為6個單元，每個單元建立一組彈簧模型。本文所用主軸與刀柄接觸部分材料為鋼，在無潤滑接觸情況下，主軸與刀柄的摩擦因數(shù)μ取為0.15，結(jié)合面基礎(chǔ)特性參數(shù)α′、β′依據(jù)文獻[22]分別選取為3 262 540與0.604。

所選位姿下刀尖頻響的實驗與仿真對比結(jié)果如圖8所示，表3給出各位姿下仿真與實驗得到的幅值與固有頻率以及兩者的相對誤差。實驗所得頻響函數(shù)隨位姿變化，這是不同位姿下加工系統(tǒng)的動力學(xué)特性有所改變導(dǎo)致的。提出的預(yù)測模型將機器人各關(guān)節(jié)角度作為變量，考慮了位姿變化的影響，結(jié)果表明，提出的機器人銑削系統(tǒng)刀尖頻響預(yù)測方法計算出的頻響函數(shù)與實驗得到的頻響函數(shù)相比，固有頻率誤差最大不超過6.63%，幅值誤差最大不超過9.80%，驗證了考慮主軸-刀柄結(jié)合面接觸剛度的機器人銑削系統(tǒng)刀尖頻響預(yù)測方法的有效性。造成仿真結(jié)果與實驗結(jié)果之間誤差的原因如下：①所用加速度傳感器的自身質(zhì)量在測量時會對刀尖動力學(xué)造成一定影響；②實驗所用模態(tài)分析軟件在將加速度信號進行兩次積分算得位移信號時不可避免地產(chǎn)生計算誤差；③對模型的簡化造成一定誤差，如將末端執(zhí)行器與機器人法蘭看作完全剛性連接等。

(a)P1

(b)P2

(c)P3

(d)P4

(e)P5圖8 所選位姿下刀尖頻響實驗與仿真對比結(jié)果Fig.8 Comparison of experimental and simulation results of tool tip frequency response under selected postures

表3 由模態(tài)實驗與數(shù)值模型得到的刀尖頻響結(jié)果對比Tab.3 The comparison of tool-tip frequency responses obtained by modal experiment and numerical model

3 結(jié)論

通過機器人本體動力學(xué)模型與末端結(jié)合面動力學(xué)模型疊加，提出一種考慮主軸-刀柄結(jié)合面接觸剛度的機器人銑削系統(tǒng)刀尖頻響預(yù)測方法，并通過加工系統(tǒng)的刀尖模態(tài)實驗進行驗證，結(jié)果表明，隨機位姿下實驗與模型仿真的誤差都能控制在10%以內(nèi)，說明所述方法可準確預(yù)測機器人銑削系統(tǒng)的刀尖頻響函數(shù)，減小了對大量實驗結(jié)果的依賴，并在保證足夠預(yù)測精度的基礎(chǔ)上簡化了模型建立過程，提高了計算效率。