陳君杰，李攀磊，韓 威，許楊劍，王效貴

(浙江工業(yè)大學(xué) 機(jī)械工程學(xué)院，浙江 杭州 310014)

0 引 言

1985年，瑞士的Clavel[1]發(fā)明了Delta并聯(lián)機(jī)器人，該型機(jī)器人為三自由度空間平移機(jī)構(gòu)，具有承載能力強(qiáng)、運(yùn)動耦合弱、力控制容易等優(yōu)點。隨著并聯(lián)機(jī)器人的應(yīng)用領(lǐng)域不斷得到拓展，其工作環(huán)境日趨復(fù)雜，并聯(lián)機(jī)器人不斷向高速度、高加速度、高精度、重載荷和輕量化方向發(fā)展[2-3]，導(dǎo)致機(jī)構(gòu)運(yùn)行中彈性振動和運(yùn)動誤差也隨之增加。傳統(tǒng)的剛體動力學(xué)分析方法無法滿足彈性誤差分析的需求，考慮構(gòu)件彈性的動力學(xué)分析成了研究重點。通過運(yùn)動彈性動力學(xué)分析方法(kineto-elasto-dynamic，KED)，將機(jī)構(gòu)位移視作彈性位移與剛體位移(名義位移)的疊加，在給定機(jī)構(gòu)名義運(yùn)動條件規(guī)律的條件下，確定機(jī)構(gòu)的彈性響應(yīng)。

Piras[4]利用有限元理論與彈性動力分析方法(KED)研究了3-PRR平面并聯(lián)機(jī)器人的彈性動力學(xué)問題。劉善增等人[5]建立了剛?cè)狁詈喜⒙?lián)機(jī)構(gòu)系統(tǒng)的整體動力學(xué)方程的步驟與方法，對3-RRS并聯(lián)機(jī)器人的頻率特性進(jìn)行了分析。韓亞峰等人[6]利用有限元理論，采用平面梁單元對Delta機(jī)器人進(jìn)行了彈性動力學(xué)建模。Kuo等人[7]基于D-H方法定義了一組全局變量，在不使用約束方程的情況下，導(dǎo)出了Delta機(jī)器人的彈性動力學(xué)模型。目前，Delta并聯(lián)機(jī)器人的彈性動力學(xué)研究中大都認(rèn)為其四邊形從動臂機(jī)構(gòu)在運(yùn)動過程中兩側(cè)桿保持平行，從而將其四邊形從動臂機(jī)構(gòu)簡化為一根虛擬從動桿進(jìn)行分析，而在實際運(yùn)動過程中，由于驅(qū)動桿末端彈性轉(zhuǎn)角等影響，四邊形機(jī)構(gòu)會產(chǎn)生扭曲，有必要針對機(jī)器人的四邊形機(jī)構(gòu)進(jìn)行動力學(xué)建模研究。

本研究將Delta機(jī)器人四邊形從動臂機(jī)構(gòu)劃分為2剛性短桿與2柔性從動桿，通過分析其運(yùn)動協(xié)調(diào)條件，在有限元理論基礎(chǔ)上建立Delta機(jī)器人的彈性動力學(xué)方程，通過數(shù)值方法求解出機(jī)器人在運(yùn)動軌跡中的彈性誤差，并通過改變桿件截面尺寸，分析桿件截面尺寸對彈性誤差的影響。

1 彈性動力學(xué)方程

Delta并聯(lián)機(jī)器人結(jié)構(gòu)圖如圖1所示。

圖1 Delta并聯(lián)機(jī)器人結(jié)構(gòu)圖

Delta機(jī)器人由定平臺、驅(qū)動器、驅(qū)動桿、四邊形從動臂、動平臺等組成，其中連接動平臺與定平臺的3個支鏈互成120°角對稱分布。每條支鏈包含一個驅(qū)動桿與一個從動臂，驅(qū)動桿一端通過驅(qū)動電機(jī)與定平臺相聯(lián)接，另一端以轉(zhuǎn)動副形式與從動臂相連，從動臂為平行四邊形結(jié)構(gòu)，確保動平臺在工作空間內(nèi)做三維平動。因為構(gòu)件的柔性以及運(yùn)動過程中慣性力和外載荷的影響，動平臺中心點P的實際位置相對于理想位置產(chǎn)生偏差，即彈性位置誤差，通過建立Delta機(jī)器人的彈性動力學(xué)方程可以對其彈性位置誤差進(jìn)行計算。

1.1 單元劃分及位移分析

根據(jù)有限元理論，筆者選擇矩形截面梁單元作為基本單元，用以劃分機(jī)構(gòu)中的柔性桿件，空間梁單元模型如圖2所示。

圖2 空間梁單元模型

單元包含2個節(jié)點。每個節(jié)點有6個彈性位移自由度，表示梁單元彈性位移的廣義坐標(biāo)：

(1)

式中：x1，x2—兩節(jié)點在x軸向的彈性位移；y1，y2—兩節(jié)點沿y向彈性位移；z1，z2—兩節(jié)點沿z向彈性位移；φx1，φx2—兩節(jié)點繞x軸彈性轉(zhuǎn)角；φy1，φy2—兩節(jié)點繞y軸彈性轉(zhuǎn)角；φz1，φz2—兩節(jié)點繞z軸彈性轉(zhuǎn)角。

根據(jù)歐拉-伯努利梁理論，梁單元上任意一點彈性位移可以表示成如下形式：

(2)

(3)

(4)

(5)

1.2 單元彈性動力學(xué)方程

因為單元彈性變形較小，忽略機(jī)構(gòu)剛體運(yùn)動與彈性變形運(yùn)動之間的耦合影響，單元的位移看作是剛體位移與彈性位移的疊加，單元動能為：

(6)

(7)

式中：ρ—單元質(zhì)量密度；L—梁單元長度；A—梁單元截面面積；IP—梁單元橫截面對x軸的極慣性矩。

單元的變形能包括彎矩、軸向力和扭矩作用時所產(chǎn)生的能量，單元總變形能為：

(8)

式中：Ke—單元剛度矩陣。

(9)

式中：E—梁單元材料的楊氏模量；G—梁單元材料的剪切模量；Iy，Iz—梁單元橫截面對y軸和z軸的極慣性矩，形函數(shù)下標(biāo)中的x以及xx分別代表對x的一階偏導(dǎo)和二階偏導(dǎo)，例如：Nux代表Nu對x的一階偏導(dǎo)。

將式(6，8)代入拉格朗日動力學(xué)方程，導(dǎo)出單元彈性動力學(xué)方程：

(10)

Fe=Qe+Pe+Ge

(11)

1.3 運(yùn)動協(xié)調(diào)關(guān)系

1.3.1 支鏈內(nèi)坐標(biāo)系建立

鑒于Delta機(jī)器人支鏈的對稱性，本研究選取一條支鏈進(jìn)行彈性動力學(xué)建模。支鏈有限元模型如圖3所示。

圖3 支鏈有限元模型

圖3中，A為驅(qū)動端，驅(qū)動桿AB被視為空間懸臂梁，B1B2D2D1為支鏈的平行四邊形從動臂結(jié)構(gòu)，P為機(jī)器人末端，因為平臺剛度遠(yuǎn)大于空間梁單元機(jī)構(gòu)，視動平臺與定平臺為剛性體。而在四邊形機(jī)構(gòu)中，上下短桿長度遠(yuǎn)小于兩側(cè)桿，因此本研究將B1B2、D1D2兩桿視作剛性體，并忽略其質(zhì)量影響，將B1D1、B2D2兩桿作為彈性桿件進(jìn)行分析。

本研究在圖3中做驅(qū)動桿AB直線在B1B2D2D1平面上的投影BC，并如圖中定義角度n1、n2、n3。建立支鏈O-XYZ坐標(biāo)系，方向定義為Z軸向上，Y軸平行于驅(qū)動器轉(zhuǎn)動軸線，X軸遵守右手定則指向支鏈方向；動平臺坐標(biāo)系P-x3y3z3方向與支鏈坐標(biāo)系O-XYZ一致；驅(qū)動桿AB單元坐標(biāo)系A(chǔ)-x1y1z1的坐標(biāo)軸方向由支鏈坐標(biāo)系O-XYZ繞Y軸旋轉(zhuǎn)n1得到；從動桿B1D1單元坐標(biāo)系B1-x2y2z2方向由支鏈坐標(biāo)系O-XYZ先繞Y軸旋轉(zhuǎn)n1+n2，后繞Z軸旋轉(zhuǎn)n3得到。

定義支鏈彈性位移廣義坐標(biāo)為：

(12)

其中，前6項元素組成P點在坐標(biāo)系P-x3y3z3下廣義坐標(biāo)UP，描述動平臺因為機(jī)構(gòu)彈性變形影響，P點相對于名義位置的位移；7～12項元素組成廣義坐標(biāo)UB，對應(yīng)AB梁單元坐標(biāo)系A(chǔ)-x1y1z1下B節(jié)點處的彈性位移；φ1-8分別為B1D1、B2D2兩桿端點在各自單元坐標(biāo)系下繞y軸與z軸方向的彈性轉(zhuǎn)角。

1.3.2 支鏈內(nèi)運(yùn)動協(xié)調(diào)關(guān)系

驅(qū)動桿AB為空間懸臂梁，點A處的彈性位移與轉(zhuǎn)角均為零，可以得出AB梁單元與支鏈彈性位移之間的關(guān)系為：

δAB=SABΨ

(13)

其中：

(14)

式中：I—單位矩陣；0—零矩陣；δAB—AB梁單元廣義坐標(biāo)。

記D點在動平臺坐標(biāo)系P-x3y3z3下廣義坐標(biāo)為：

(15)

則D點與P點位移協(xié)調(diào)關(guān)系為[8]：

(16)

根據(jù)圖3中幾何約束關(guān)系可以得到B1D1梁單元在單元坐標(biāo)系下的彈性位移：

x1=(xBc2+zBs2-d(φzBc2-φxBs2))c3+yBs3；

y1=-(xBc2+zBs2-d(φzBc2-φxBs2))s3+yBc3；

z1=zBc2-xBs2+d(φxBc2+φzBs2)；

φx1=0；φy1=φ1；φz1=φ2；

x2=(xDc12+zPs12-d(φzDc12-φxDs12))c3+yDs3；

y2=-(xDc12+zDs12-d(φzDc12-φxDs12))s3+yDc3；

z2=zDc12-xDs12+d(φxDc12+φzDs12)；

φx2=0；φy2=φ3；φz2=φ4

(17)

式中：d=1/2|B1B2|；符號s，c—函數(shù)sin和cos，下標(biāo)數(shù)字對應(yīng)3個角度n1、n2、n3，例如：s12代表sin(n1+n2)，c3代表cosn3。

由式(16，17)可以得出B1D1梁單元與支鏈彈性位移之間的協(xié)調(diào)關(guān)系：

δBD1=SBD1Ψ

(18)

式中：δBD1—B1D1梁單元廣義坐標(biāo)；SBD1—B1D1梁坐標(biāo)協(xié)調(diào)矩陣，同理可以得出B2D2梁單元與支鏈彈性位移之間的協(xié)調(diào)關(guān)系：

δBD2=SBD2Ψ

(19)

動平臺P與支鏈彈性位移協(xié)調(diào)關(guān)系為：

UP=SPΨ

(20)

式中：

1.3.3 系統(tǒng)運(yùn)動協(xié)調(diào)關(guān)系

本研究建立系統(tǒng)廣義坐標(biāo)U∈R48×1，根據(jù)式(13，18，19，20)，以及3條支鏈的對稱性，建立出任意構(gòu)件i與系統(tǒng)廣義坐標(biāo)U之間的協(xié)調(diào)關(guān)系：

δi=SiU

(21)

式中：i—構(gòu)件編號；δi—構(gòu)件i的單元廣義坐標(biāo)；Si—對應(yīng)的坐標(biāo)協(xié)調(diào)矩陣。

1.4 系統(tǒng)彈性動力學(xué)方程

動平臺為剛體，其動力學(xué)方程[9-10]可以表示為：

(22)

其中，動平臺質(zhì)量矩陣：

(23)

將式(19)代入式(10，22)得到各構(gòu)件動力學(xué)方程：

(24)

將各構(gòu)件的單元動力學(xué)方程(24)進(jìn)行總裝得到：

(25)

2 彈性誤差分析

2.1 彈性誤差計算

筆者選用系統(tǒng)參數(shù)：驅(qū)動桿長500 mm，動平臺質(zhì)量0.2 kg，上平臺外接圓半徑100 mm，動平臺外接圓半徑為50 mm，從動臂結(jié)構(gòu)中短桿長度為50 mm，兩側(cè)的從動桿長600 mm，驅(qū)動桿截面與從動桿截面均選用正方形截面，驅(qū)動桿截面邊長尺寸20 mm，從動桿截面邊長尺寸10 mm，材料密度為7 850 kg/m3，彈性模量210 GPa，泊松比0.3，運(yùn)行時間T=3s。給定動平臺運(yùn)動軌跡：

(26)

本研究利用Newmark方法在Matlab中對系統(tǒng)動力學(xué)方程(25)進(jìn)行數(shù)值求解，計算出動平臺末端P點的彈性位置誤差在運(yùn)行時間內(nèi)的變化情況。

圓周軌跡下彈性位置誤差如圖4所示。

圖4 圓周軌跡下彈性位置誤差Δx、Δy、Δz—P點沿x、y、z三軸方向的彈性誤差。

2.2 桿截面尺寸對彈性誤差影響

通過改變驅(qū)動桿與從動桿的截面參數(shù)，可以有效改變機(jī)器人的力學(xué)性能。定義軌跡上誤差均值：

(27)

式中：Δx(t)，Δy(t)，Δz(t)—P點在t時刻沿各軸向的彈性誤差值。

2.2.1 驅(qū)動桿截面尺寸對彈性誤差影響

驅(qū)動桿截面邊長選擇取值20 mm～30 mm的范圍，其他參數(shù)不變。本研究根據(jù)式(27)計算誤差均值隨驅(qū)動桿的截面邊長的變化情況。結(jié)果如圖5所示。

圖5 驅(qū)動桿截面尺寸對機(jī)構(gòu)末端彈性誤差的影響

可以看到：隨著驅(qū)動桿截面邊長增大，機(jī)器人末端的誤差均值明顯減少。究其原因，驅(qū)動桿的受力形式類似于懸臂梁，增加其截面尺寸能夠有效提升其彎曲剛度。

2.2.2 從動桿截面尺寸對彈性誤差影響

考察從動桿截面尺寸對機(jī)構(gòu)彈性誤差的影響，從動桿截面邊長選擇取值5 mm～10 mm的范圍，其他參數(shù)不變。誤差均值隨從動桿截面邊長的變化情況如圖6所示。

圖6 從動桿截面尺寸對機(jī)構(gòu)末端彈性誤差的影響

由圖可以看到：當(dāng)從動桿截面邊長增大時，機(jī)構(gòu)的均值誤差隨之增加。分析原因，從動桿兩端均為球鉸關(guān)節(jié)，其變形方式以拉壓為主，相對于彎曲變形，軸向拉壓變形的尺寸相對較小，其剛度的提升并不能抵消掉由質(zhì)量增加帶來的額外載荷影響。

2.2.3 兩桿截面尺寸對彈性誤差綜合影響

筆者考察兩種桿截面尺寸對機(jī)構(gòu)彈性誤差的綜合影響，從動桿截面邊長選擇取值5 mm～10 mm的范圍，驅(qū)動桿截面邊長選擇取值20 mm～30 mm的范圍，其他參數(shù)不變，機(jī)器人均值誤差隨兩桿截面尺寸變化如圖7所示。

圖7 兩類桿截面尺寸對機(jī)構(gòu)末端彈性誤差的綜合影響

可以看到在區(qū)間內(nèi)均值誤差的變化規(guī)律，誤差均值與驅(qū)動桿截面尺寸呈負(fù)相關(guān)，與從動桿截面尺寸呈正相關(guān)，具有一定的單調(diào)性。

3 結(jié)束語

本研究針對Delta機(jī)器人運(yùn)動過程中的彈性變形誤差進(jìn)行了分析，建立了系統(tǒng)的彈性動力學(xué)控制方程，通過數(shù)值方法進(jìn)行算例分析，求解出了機(jī)器人在運(yùn)動過程中的誤差情況，分析了桿件截面尺寸對彈性誤差的影響。

結(jié)論顯示：通過增加驅(qū)動桿的截面尺寸以及減少從動臂的截面尺寸，能夠有效降低機(jī)器人運(yùn)動過程中的彈性變形。

[1] CLAVEL R, A fast robot with parallel geometry[C]. Proc.Int.symposium. on Industrial Robots, Lausanne: CiNii,1988.

[2] 計時鳴,黃希歡.工業(yè)機(jī)器人技術(shù)的發(fā)展與應(yīng)用綜述[J].機(jī)電工程,2015,32(1):1-13.

[3] 馮李航,張為公,龔宗洋,等.Delta系列并聯(lián)機(jī)器人研究進(jìn)展與現(xiàn)狀[J].機(jī)器人,2014(3):375-384.

[4] PIRAS G, CLEGHORN W L, MILLS J K. Dynamic finite-element analysis of a planar high-speed, high-precision parallel manipulator with flexible links[J].Mechanism&MachineTheory,2005,40(7):849-862.

[5] 劉善增，朱真才,余躍慶,等.空間剛?cè)狁詈喜⒙?lián)機(jī)構(gòu)系統(tǒng)的頻率特性分析[J].機(jī)械工程學(xué)報,2011,47(23):39-48.

[6] 韓亞鋒,馬履中,吳偉光,等.Delta并聯(lián)機(jī)器人彈性動力學(xué)研究[J].農(nóng)業(yè)機(jī)械學(xué)報,2011,42(10):197-202.

[7] KUO Y L. Mathematical modeling and analysis of the Delta robot with flexible links[J].Computers&MathematicswithApplications,2016,71(10):1973-1989.

[8] 黃 真,孔令富,方躍法.并聯(lián)機(jī)器人機(jī)構(gòu)學(xué)理論及控制[M].北京：機(jī)械工業(yè)出版社,1997.

[9] 韓敬虎，俞經(jīng)虎.食品檢測咀嚼機(jī)器人工作空間研究[J].輕工機(jī)械，2016，34(3)：1-4.

[10] 熊艷梅，楊延棟.碼垛機(jī)器人運(yùn)動學(xué)分析與仿真[J].機(jī)械，2015(12)：62-66.

[11] 巴 特.有限元分析中的數(shù)值方法[M].北京:科學(xué)出版社,1985.