軟體連續(xù)機(jī)械臂動(dòng)力學(xué)建模與仿真1)

白爭(zhēng)鋒 孔清峰 趙 起

(哈爾濱工業(yè)大學(xué)(威海)機(jī)械工程系,山東威海 264209)

引言

軟體機(jī)器人以其無(wú)限自由度、被動(dòng)變形適應(yīng)環(huán)境和柔順接觸等優(yōu)勢(shì)在醫(yī)學(xué)、航天、農(nóng)業(yè)、汽車(chē)、制造、核能、醫(yī)療保健和康復(fù)、環(huán)境保護(hù)和海洋勘探、反恐偵察等方面有著廣闊的應(yīng)用前景[1-3].在過(guò)去十多年中,軟體機(jī)器人成為眾多學(xué)者研究的對(duì)象.仿象鼻、章魚(yú)觸手等結(jié)構(gòu)的機(jī)械臂式連續(xù)體機(jī)器人,以其固有的柔順性,安全性及適應(yīng)非結(jié)構(gòu)化環(huán)境的優(yōu)良特性使其顯示出巨大的應(yīng)用潛力,成為目前機(jī)器人領(lǐng)域的熱門(mén)研究方向之一[4-8].

目前,研究人員針對(duì)軟體機(jī)械臂進(jìn)行了一系列的研究,并取得了一些研究成果.Walker 等[9-13]對(duì)軟體機(jī)械臂建模做了大量研究,Renda 等[14]設(shè)計(jì)的仿章魚(yú)軟體機(jī)械臂樣機(jī),基于Cosserat 梁理論將軟體機(jī)械臂軸線表示為空間中的一維曲線,并考慮了纜索驅(qū)動(dòng)機(jī)器人的相關(guān)彈性力、驅(qū)動(dòng)力和重力,建立了肌腱驅(qū)動(dòng)的軟體機(jī)械臂運(yùn)動(dòng)學(xué)及靜力學(xué)模型;基于Kelvin-Voight 線性黏彈性本構(gòu)方程建立Cosserat 應(yīng)變模型[15],建立軟體章魚(yú)臂機(jī)器人在水下環(huán)境中的動(dòng)力學(xué)模型;其后提出一種基于離散化Cosserat 理論的分段恒應(yīng)變[16]方法,并建立了軟體機(jī)械臂的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)模型.然而該方法推導(dǎo)出的動(dòng)力學(xué)方程是一組耦合偏微分方程,在數(shù)值求解上十分困難.Kang 等[17-19]利用章魚(yú)生理解剖結(jié)構(gòu),基于分布式的彈簧-阻尼模型模擬章魚(yú)臂4 個(gè)縱向和4 個(gè)徑向肌肉,建立了章魚(yú)機(jī)械臂的三維動(dòng)力學(xué)模型,但是該模型沒(méi)有從能量的角度考慮集中質(zhì)量模型,并且使用歐拉遞推法建立動(dòng)力學(xué)方程,方程推導(dǎo)十分繁瑣.Rone 和Ben-tzvi[20]使用一組有限的運(yùn)動(dòng)學(xué)變量來(lái)捕獲機(jī)械臂的曲率變化,基于虛功原理為肌腱驅(qū)動(dòng)的軟體機(jī)械臂建立了動(dòng)力學(xué)方程;在推導(dǎo)動(dòng)力學(xué)方程中考慮了慣性、致動(dòng)、摩擦、彈性和重力等效應(yīng).然而該方法推導(dǎo)連續(xù)臂的動(dòng)力學(xué)方程會(huì)導(dǎo)致復(fù)雜的非線性積分,除了固有的數(shù)值不穩(wěn)定性之外,這些積分在計(jì)算上是不可行的,因此該方法有一定的局限性.

Tatlicioglu 等[21]基于恒定曲率和分布式質(zhì)量建立了可變長(zhǎng)度連續(xù)體機(jī)械手平面動(dòng)力學(xué)模型,并考慮了重力勢(shì)能項(xiàng)和彈性勢(shì)能項(xiàng),使模型與機(jī)械臂的真實(shí)行為相一致.Falkenhahn 等[22-23]提出了一類稱為仿生操作助手的連續(xù)機(jī)械手的動(dòng)力學(xué)模型,該模型采用歐拉-拉格朗日方法,并基于分段常曲率假設(shè),將每個(gè)段的質(zhì)量看作集中質(zhì)量放置在每段的末端位置.該模型雖然在數(shù)值計(jì)算上可行,但是將每段質(zhì)量看作位于頭部的集中質(zhì)量并不能從能量上與實(shí)際的分布式質(zhì)量相匹配.Giri 和Walker[24]建立了三段式平面軟體機(jī)械臂動(dòng)力學(xué)模型,然而該模型將每段的質(zhì)量看作位于末端的集中質(zhì)量.針對(duì)軟體機(jī)器人變形后呈圓弧形狀彎曲,Mazzolaii 等[25]提出分段常曲率(PCC)理論模型.對(duì)于變曲率彎曲問(wèn)題,可將軟體機(jī)器人本體結(jié)構(gòu)分為多段,假設(shè)每段曲率恒定,每段形狀可通過(guò)曲率半徑、兩截面所在平面間夾角以及兩截面扭轉(zhuǎn)角來(lái)描述[26-27],然后建立末端截面與坐標(biāo)系的關(guān)系.

軟體機(jī)械臂動(dòng)力學(xué)模型數(shù)值求解精度與計(jì)算效率對(duì)于軟體機(jī)械臂的動(dòng)態(tài)控制至關(guān)重要.本文基于模態(tài)方法推導(dǎo)軟體機(jī)械臂運(yùn)動(dòng)方程,明確地從能量的角度考慮軟體機(jī)械臂動(dòng)力學(xué)特性,提出一種基于質(zhì)心集中質(zhì)量的軟體機(jī)械臂動(dòng)力學(xué)模型,該模型將軟體機(jī)械臂的分布質(zhì)量等效為位于質(zhì)心的集中質(zhì)量,通過(guò)動(dòng)能等效系數(shù)實(shí)現(xiàn)二者的動(dòng)能匹配.與將集中質(zhì)量放置在任意位置(例如,軟體機(jī)械臂的末端)的集中質(zhì)量模型相比,該模型兼顧了采用連續(xù)分布質(zhì)量進(jìn)行拉格朗日方法分析時(shí)的準(zhǔn)確性,以及采用離散化的集中質(zhì)量模型分析時(shí)的計(jì)算高效性.仿真結(jié)果表明,本文提出的將集中質(zhì)量放置在與機(jī)械臂的能量行為相匹配的位置所建立的軟體機(jī)械臂集中質(zhì)量動(dòng)力學(xué)模型,能夠準(zhǔn)確、高效地捕獲軟體機(jī)械臂的動(dòng)力學(xué)特性,并且數(shù)值計(jì)算穩(wěn)定.

1 運(yùn)動(dòng)學(xué)建模

1.1 軟體機(jī)械臂模型抽象

將軟體機(jī)械臂抽象為如圖1 所示的結(jié)構(gòu),軟體機(jī)械臂由三個(gè)結(jié)構(gòu)完全相同的驅(qū)動(dòng)器組成,且均布在間距為120°以r為半徑的輕質(zhì)圓形剛性框架上,剛性框架實(shí)現(xiàn)了三個(gè)驅(qū)動(dòng)器在空間上隔離與固定.結(jié)構(gòu)中也等間隔分布圓形剛性框架,以保持驅(qū)動(dòng)器軸線方向始終平行于中性軸,假設(shè)剛性框架相對(duì)于驅(qū)動(dòng)器的質(zhì)量可以忽略不計(jì).三個(gè)驅(qū)動(dòng)器的初始長(zhǎng)度均為 ?,各驅(qū)動(dòng)器的總長(zhǎng)度隨時(shí)間變化規(guī)律為?+?k(t),k∈{1,2,3},k表示軟體機(jī)械臂上三個(gè)驅(qū)動(dòng)器的編號(hào),其中 ?k(t) 表示各驅(qū)動(dòng)器長(zhǎng)度的變化量,并且驅(qū)動(dòng)器長(zhǎng)度的變化量在一定的范圍內(nèi),即 ?min,?max均為定值,其物理意義是驅(qū)動(dòng)器的最大壓縮量或最大伸長(zhǎng)量.如圖1(a)描述的是軟體機(jī)械臂在空間中彎曲后的形狀,軟體機(jī)械臂在空間中彎曲成一段圓弧,C點(diǎn)表示圓弧的曲率中心,全局參考系 {R} 的原點(diǎn)O與軟體機(jī)械臂下底面的中心重合,將驅(qū)動(dòng)器1 的底端落在X軸正方向上的A點(diǎn),驅(qū)動(dòng)器2 和3 按逆時(shí)針?lè)较虿贾?分別對(duì)應(yīng)于B點(diǎn)和D點(diǎn).

圖1 空間軟體機(jī)械臂示意圖Fig.1 Schematic diagram of space soft manipulator

為了簡(jiǎn)化分析,在軟體機(jī)械臂運(yùn)動(dòng)學(xué)建模時(shí)提出以下基本假設(shè):

(1) 軟體機(jī)械臂變形過(guò)程中,各驅(qū)動(dòng)器假定為彎曲曲率相等的光滑連續(xù)曲線;

(2) 驅(qū)動(dòng)器變形過(guò)程中不考慮扭轉(zhuǎn)變形;

(3) 驅(qū)動(dòng)器僅沿長(zhǎng)度方向變化,徑向尺寸不發(fā)生改變;

(4) 軟體機(jī)械臂沿長(zhǎng)度方向具有均勻面密度.

1.2 驅(qū)動(dòng)空間到配置空間的映射

為描述軟體機(jī)械臂的運(yùn)動(dòng),引入如下坐標(biāo)系:全局參考系 {R}固連在機(jī)械臂的底端,局部參考系{R1}固連在上端,并且z1軸與中性軸相切.根據(jù)假設(shè)(1),軟體機(jī)械臂變形后在空間中的配置可以由一組曲線參數(shù)完全定義,曲率半徑λ ∈(0,+∞),瞬時(shí)曲率中心C;彎曲張角?,以及彎曲平面相對(duì)于X軸正方向的轉(zhuǎn)角 φ∈[-π,π] .為了獲得末端執(zhí)行器的位姿與驅(qū)動(dòng)器的長(zhǎng)度變化量的關(guān)系,定義三個(gè)空間與兩個(gè)映射.如圖2 所示,三個(gè)空間分別為:代表驅(qū)動(dòng)空間關(guān)節(jié)變量;代表配置空間關(guān)節(jié)變量;H代表局部參考系 {R1}到全局參考系{R}的齊次變換矩陣,x代表任務(wù)空間,即軟體機(jī)械臂末端的笛卡爾空間.兩個(gè)映射分別為:fspecific表示從驅(qū)動(dòng)空間到配置空間的映射,即 κ=κ(q),fgeneral表示從配置空間到任務(wù)空間的映射,即H=H(κ) .根據(jù)文獻(xiàn)[28]可以建立驅(qū)動(dòng)空間到配置空間的映射關(guān)系為

圖2 不同空間的運(yùn)動(dòng)學(xué)映射Fig.2 Kinematic mapping of different spaces

1.3 配置空間到任務(wù)空間的映射

考慮配置空間κ到任務(wù)空間x的映射關(guān)系,即求解局部參考系{R1}到全局參考系 {R} 的齊次變換矩陣H.齊次變換矩陣需要經(jīng)過(guò)如圖3 所示的5 個(gè)步驟得到,分別為:首先,繞全局坐標(biāo)系 {R} 的Z軸旋轉(zhuǎn)角度φ;然后沿著當(dāng)前參考系X軸正向移動(dòng)距離λ;接著繞當(dāng)前參考系的Y軸旋轉(zhuǎn)角度?;再沿著當(dāng)前參考系的X軸負(fù)方向移動(dòng)距離λ;最后繞著當(dāng)前參考系Z軸旋轉(zhuǎn)角度-φ,最終得到配置空間到任務(wù)空間的齊次變換矩陣為

圖3 齊次變換矩陣的生成步驟Fig.3 Steps of homogeneous transformation matrix

式中,Θ(κ)∈R3×3表示局部參考系 {R1} 到全局參考系{R}的旋轉(zhuǎn)矩陣,ψ(κ)∈R3×1表示局部參考系{R1} 到全局參考系 {R} 的位置平移向量.

Ψ(κ,η)∈R3×1表示動(dòng)參考系到全局參考系的位置平移向量,該向量中的3 個(gè)元素表達(dá)式如下

1.4 配置參數(shù)的局限性

式(3)給出的齊次變換矩陣包含關(guān)節(jié)空間中連續(xù)體中性軸上任意點(diǎn)的位置和方向,但該矩陣中的每個(gè)元素都是由三角函數(shù)構(gòu)成.由于三角函數(shù)屬于復(fù)雜的非線性函數(shù),并且該式中使用的配置參數(shù)為λ和φ,這兩項(xiàng)中的分母分別為和l2+l3-2l1,這將產(chǎn)生一個(gè)不可避免的問(wèn)題,即當(dāng)l1=l2=l3時(shí),這兩項(xiàng)的分母都將變?yōu)榱?隨之導(dǎo)致數(shù)值計(jì)算的奇異性與不穩(wěn)定性.如果采取一種處理方式,將齊次變換矩陣中各元素近似為數(shù)值穩(wěn)定的多元多項(xiàng)式,則可避免奇異性問(wèn)題.不同類型的級(jí)數(shù)展開(kāi)法,如冪級(jí)數(shù)、傅里葉級(jí)數(shù)以及泰勒級(jí)數(shù)均可用于簡(jiǎn)化復(fù)雜的非線性函數(shù)[29].考慮到軟體機(jī)械臂連續(xù)變形的特性,多元泰勒級(jí)數(shù)展開(kāi)不含三角項(xiàng)的多項(xiàng)式與具有驅(qū)動(dòng)空間關(guān)節(jié)變量的分母,因此能夠有效消除數(shù)值計(jì)算的奇異性與不穩(wěn)定性.

本文在描述軟體機(jī)械臂變形后的空間位置描時(shí),使用模態(tài)方法代替配置參數(shù)導(dǎo)出的運(yùn)動(dòng)學(xué)方程,從而避免奇異性.模態(tài)方法是利用簡(jiǎn)單的多項(xiàng)式函數(shù)逼近復(fù)雜的非線性函數(shù),同時(shí),利用模態(tài)坐標(biāo)將會(huì)極大提高動(dòng)力學(xué)方程求解效率.鑒于此,將式(3)展開(kāi)為三元泰勒級(jí)數(shù),齊次變換矩陣H用模態(tài)齊次變換矩陣T近似表達(dá),本文給出截?cái)嗾`差具有11 階精度的模態(tài)變換矩陣的表達(dá)式,為表示方便,令

則模態(tài)變換矩陣T中模態(tài)旋轉(zhuǎn)矩陣R的表達(dá)式如下

模態(tài)變換矩陣T中模態(tài)位置平移向量p的表達(dá)式如下

則軟體連續(xù)臂中性軸上任意點(diǎn)到全局參考系的模態(tài)齊次變換矩陣為

基于泰勒級(jí)數(shù)的模態(tài)方法提供了一個(gè)直觀和直接的選擇,使齊次變換矩陣H中的每個(gè)元素產(chǎn)生唯一且精確的多項(xiàng)式函數(shù),因此,它避免了模型轉(zhuǎn)換和模型奇異性.此外,軟體機(jī)器人正運(yùn)動(dòng)學(xué)和逆運(yùn)動(dòng)學(xué)都可以在空間關(guān)節(jié)中直接計(jì)算,而不需要額外的奇異性解決方法,并且所提出的方法可以簡(jiǎn)便地?cái)U(kuò)展到多節(jié)連續(xù)機(jī)械臂的動(dòng)力學(xué)建模中.

2 連續(xù)機(jī)械臂動(dòng)力學(xué)建模

通常認(rèn)為機(jī)械臂的質(zhì)量是沿臂的長(zhǎng)度均勻連續(xù)分布的,如圖4(a)所示.本節(jié)首先考慮機(jī)械臂的質(zhì)量是沿臂的長(zhǎng)度均勻連續(xù)分布的,基于拉格朗日方法建立軟體機(jī)械臂動(dòng)力學(xué)模型;然后提出基于質(zhì)心集中質(zhì)量的動(dòng)力學(xué)模型,機(jī)械臂的質(zhì)量位于質(zhì)心處,如圖4(b)所示.

圖4 兩種動(dòng)力學(xué)模型Fig.4 Two kinds of dynamics model

2.1 分布質(zhì)量動(dòng)力學(xué)建模

軟體機(jī)械臂的拉格朗日函數(shù)可以寫(xiě)為 L=K-P,其中 K 表示軟體機(jī)械臂的總動(dòng)能,P 表示軟體機(jī)械臂的總勢(shì)能.總動(dòng)能 K=Kυ+Kω,Kυ表示軟體機(jī)械臂的平動(dòng)動(dòng)能,Kω表示軟體機(jī)械臂的轉(zhuǎn)動(dòng)動(dòng)能.總勢(shì)能 P=Pg+Pe,Pg表示軟體機(jī)械臂的重力勢(shì)能,Pe表示軟體機(jī)械臂的彈性勢(shì)能.為獲得以上各能量項(xiàng)的表達(dá)式,如圖5 所示基于微分的思想將軟體機(jī)械臂沿中性軸看作無(wú)窮多個(gè)切片組成,中性軸上η處切片的示意圖如圖6 所示,將切片看作質(zhì)量為mdη,厚度為λ?dη,半徑為r的均質(zhì)圓盤(pán).位于η處切片的動(dòng)能與勢(shì)能得到后,通過(guò)從底部到頂部對(duì)η從0 到1 積分確定軟體機(jī)械臂的總能量.

圖5 切片示意圖Fig.5 Schematic of slice

圖6 切片放大示意圖Fig.6 Slice enlargement

2.1.1 動(dòng)能計(jì)算

為確定軟體機(jī)械臂的總動(dòng)能,采用式(9)推導(dǎo)的模態(tài)變換矩陣獲得切片的瞬時(shí)速度,包括線速度ν與角速度ω.將式(8)對(duì)時(shí)間求導(dǎo)可得中性軸上任意點(diǎn)η處切片相對(duì)于全局參考系 {R} 的線速度為

式中,線速度雅可比矩陣Jn表達(dá)式為

中性軸上點(diǎn)η處切片在體參考系 {R1} 中的角速度為

式中,操作符 ∨ 表示對(duì)反對(duì)稱矩陣提取非零元素,且角速度的三個(gè)分量為

由于假定軟體機(jī)械臂的質(zhì)量恒定且具有均勻的面密度,因此切片質(zhì)量mdη具有的線慣性矩陣為Mυ=mdηI3,I3為三階單位矩陣.又由于假定切片具有圓形截面,因此其具有的慣性張量為Mω=Ixxdiag(1,1,2),其中表示對(duì)角矩陣.將質(zhì)量相關(guān)的切片能量對(duì) η 進(jìn)行積分,則軟體機(jī)械臂的平動(dòng)動(dòng)能為

同理,軟體機(jī)械臂的轉(zhuǎn)動(dòng)動(dòng)能為

則軟體機(jī)械臂總動(dòng)能為

2.1.2 勢(shì)能計(jì)算

軟體機(jī)械臂的勢(shì)能主要包括兩部分,分別為重力勢(shì)能和彈性勢(shì)能,依次進(jìn)行計(jì)算.中性軸上任意點(diǎn)η處切片具有的質(zhì)量為mdη,則其重力勢(shì)能為

故整個(gè)軟體機(jī)械臂的重力勢(shì)能為

式中,g=[0,0,g]T,g表示重力加速度常數(shù).

彈性勢(shì)能是由于三個(gè)驅(qū)動(dòng)器的伸縮效應(yīng)引起的.驅(qū)動(dòng)器的彈性勢(shì)能需要根據(jù)驅(qū)動(dòng)器所使用的不同的材料建立不同的函數(shù)關(guān)系[30].一般而言,彈性勢(shì)能是變形量的函數(shù),即 Pe=f(q) .本文忽略驅(qū)動(dòng)器材料的剪切變形,假設(shè)驅(qū)動(dòng)器材料是線彈性的,即驅(qū)動(dòng)器的彈性勢(shì)能正比于伸長(zhǎng)量的平方.因此,軟體機(jī)械臂的剛度矩陣可以寫(xiě)為

從而軟體機(jī)械臂的彈性勢(shì)能為

獲得以上各能量項(xiàng)后,則軟體機(jī)械臂的拉格朗日方程為

進(jìn)而,該軟體機(jī)械臂驅(qū)動(dòng)關(guān)節(jié)空間中的拉格朗日動(dòng)力學(xué)方程以矩陣形式表示為

式中,M∈R3×3為機(jī)器人的質(zhì)量矩陣;C∈R3×3為離心力與科里奧利矩陣;G∈R3×3為保守力矢量;Q∈R3×3為外力矢量.

質(zhì)量矩陣為

保守力矩陣為

2.2 能量分析與方程簡(jiǎn)化

由于軟體機(jī)械臂的軸向尺寸遠(yuǎn)大于其徑向尺寸,本節(jié)對(duì)軟體機(jī)械臂的平動(dòng)動(dòng)能和轉(zhuǎn)動(dòng)動(dòng)能進(jìn)行定量評(píng)估.平動(dòng)動(dòng)能與轉(zhuǎn)動(dòng)動(dòng)能是驅(qū)動(dòng)空間變量q和驅(qū)動(dòng)空間速度的函數(shù),但驅(qū)動(dòng)空間變量受驅(qū)動(dòng)力的影響,隨時(shí)間變化的規(guī)律多樣,導(dǎo)致是變化多樣的,盡管q與變化規(guī)律不定,但受到材料物理極限的限制,兩者的變化量必位于某一范圍內(nèi).假定軟體機(jī)械臂驅(qū)動(dòng)器的長(zhǎng)度變化范圍為 [0,30 mm],并且其速度變化范圍為 [-50,50 mm/s],在此條件下,通過(guò)使用全局搜索計(jì)算得到 0≤Kω/K ≤3% .

進(jìn)而,采用一具體的數(shù)值算例對(duì)動(dòng)能進(jìn)行分析,以便直觀的揭示轉(zhuǎn)動(dòng)動(dòng)能 Kω占總動(dòng)能 K 的比值.假定軟體機(jī)械臂的三個(gè)驅(qū)動(dòng)器長(zhǎng)度變化量隨時(shí)間變化規(guī)律為,0≤t≤2π,根據(jù)驅(qū)動(dòng)器的變化規(guī)律,首先計(jì)算軟體機(jī)械臂末端的線速度與角速度,分別如圖7 和圖8 所示.

圖8 角速度變化曲線Fig.8 Angular velocity

由圖7 可以看出線速度的大小為每秒幾十毫米,角速度的大小為1 rad/s 左右,這顯然符合正常的速度值范圍,因此無(wú)法通過(guò)線速度與角速度的大小關(guān)系來(lái)判斷平動(dòng)動(dòng)能和轉(zhuǎn)動(dòng)動(dòng)能的大小關(guān)系.值得注意的是,在推導(dǎo)二者表達(dá)式時(shí)得到軟體機(jī)械臂的平移質(zhì)量矩陣為 Mυ=mI3,轉(zhuǎn)動(dòng)質(zhì)量矩陣為Mω=,兩者的數(shù)量級(jí)是不同的,進(jìn)而可以定性的判定出平動(dòng)動(dòng)能要遠(yuǎn)大于轉(zhuǎn)動(dòng)動(dòng)能.

圖7 線速度變化曲線Fig.7 Linear speed

進(jìn)一步經(jīng)定量計(jì)算得到二者隨時(shí)間的變化規(guī)律如圖9 所示,為 Kω和 Kυ隨時(shí)間變化曲線,從圖9 中可以直觀地看出轉(zhuǎn)動(dòng)動(dòng)能遠(yuǎn)小于平動(dòng)動(dòng)能.如圖10所示為 Kω占總動(dòng)能 K 的比值隨時(shí)間的變化圖,可以發(fā)現(xiàn)兩者比值最大值不足 1.2%,由此可見(jiàn),轉(zhuǎn)動(dòng)動(dòng)能Kω占總動(dòng)能 K 的比值很小.另外,根據(jù)式(11)和式(12)可以看出,角速度雅可比Jω的計(jì)算量要比線速度雅可比Jυ的計(jì)算復(fù)雜的多,然而其對(duì)整個(gè)動(dòng)能的占比又非常小,因此為了簡(jiǎn)化計(jì)算,在建模過(guò)程中轉(zhuǎn)動(dòng)動(dòng)能可以忽略不計(jì).本文下述動(dòng)力學(xué)分析過(guò)程中,均不考慮轉(zhuǎn)動(dòng)動(dòng)能,即 K≈Kυ.

圖9 動(dòng)能變化曲線Fig.9 Kinetic energy

圖10 轉(zhuǎn)動(dòng)動(dòng)能占總動(dòng)能百分比Fig.10 Percentage of rotational kinetic energy

如圖11 所示為軟體機(jī)械臂的彈性勢(shì)能以及重力勢(shì)能隨時(shí)間的變化規(guī)律.初始時(shí),軟體機(jī)械臂處于未致動(dòng)狀態(tài),此時(shí)三個(gè)驅(qū)動(dòng)器的長(zhǎng)度均未發(fā)生變化,因此彈性勢(shì)能為零.但重力勢(shì)能恰好為mg?/2=0.112 5 J,隨著驅(qū)動(dòng)器長(zhǎng)度的變化,重力勢(shì)能的變化范圍并未出現(xiàn)較大的波動(dòng),這是因?yàn)檐涹w機(jī)械臂的質(zhì)心位置在Z 軸方向上的位移沒(méi)有發(fā)生太大變化,該現(xiàn)象從圖12 也能明顯發(fā)現(xiàn).彈性勢(shì)能的變化范圍較大,這主要是驅(qū)動(dòng)器長(zhǎng)度變化范圍引起的.需要特別指出的是,彈性勢(shì)能與重力勢(shì)能的大小在同一數(shù)量級(jí)上,這一點(diǎn)不同于平動(dòng)動(dòng)能和轉(zhuǎn)動(dòng)動(dòng)能,因此,二者均需要在動(dòng)力學(xué)方程中考慮.

圖11 勢(shì)能比較Fig.11 Potential energy comparison

圖12 質(zhì)心軌跡Fig.12 Trajectory of mass center

通過(guò)以上對(duì)軟體機(jī)械臂能量進(jìn)行計(jì)算分析得出如下結(jié)論:第一,平動(dòng)動(dòng)能遠(yuǎn)遠(yuǎn)大于轉(zhuǎn)動(dòng)動(dòng)能,因此在動(dòng)力學(xué)方程中忽略轉(zhuǎn)動(dòng)動(dòng)能;第二,重力勢(shì)能與彈性勢(shì)能處于同一數(shù)量級(jí),二者均要在動(dòng)力學(xué)方程中考慮.

2.3 基于質(zhì)心集中質(zhì)量動(dòng)力學(xué)建模

考慮基于質(zhì)心集中質(zhì)量的動(dòng)力學(xué)模型,質(zhì)心在動(dòng)力學(xué)分析中具有重要的地位,特別是在剛體力學(xué)中,其地位顯得尤為重要.事實(shí)上在處理大變形問(wèn)題時(shí),也可以基于質(zhì)心進(jìn)行分析,在軟體機(jī)械臂建模方面,已經(jīng)提出了許多集中參數(shù)模型.然而,近似軟體機(jī)械臂的平滑彎曲需要將其劃分成很多小段或集中質(zhì)量點(diǎn),這顯著增加了模型的自由度和計(jì)算的復(fù)雜性[31],迄今為止,集中質(zhì)量建模研究的重點(diǎn)主要是降低計(jì)算復(fù)雜性,而沒(méi)有充分考慮各種近似值在能量域中是否具有良好的匹配精度.本小節(jié)將從能量的角度考慮軟體機(jī)器臂動(dòng)力學(xué),將集中質(zhì)量放置在與機(jī)械臂的能量行為相匹配的位置(即質(zhì)心位置).與將集中質(zhì)量放置在任意位置(例如每個(gè)部分的尖端)的集中質(zhì)量模型相比,基于質(zhì)心集中質(zhì)量的動(dòng)力學(xué)模型更接近基于分布質(zhì)量建立的動(dòng)力學(xué)模型,且沒(méi)有積分項(xiàng),計(jì)算效率高,并且數(shù)值穩(wěn)定.

2.3.1 基于質(zhì)心的集中質(zhì)量動(dòng)力學(xué)建模

基于質(zhì)心集中質(zhì)量模型對(duì)軟體機(jī)械臂進(jìn)行動(dòng)力學(xué)分析,隨著驅(qū)動(dòng)器長(zhǎng)度的變化,軟體機(jī)械臂在運(yùn)動(dòng)時(shí),其構(gòu)型時(shí)刻發(fā)生改變,導(dǎo)致其質(zhì)心位置不斷發(fā)生改變,并且軟體機(jī)械臂在運(yùn)動(dòng)過(guò)程中,其質(zhì)心位置不一定位于幾何體上.圖12 顯示了隨著機(jī)械臂的構(gòu)型發(fā)生變化,其質(zhì)心的位置也在不斷改變.

盡管質(zhì)心的位置在不斷變化,但可以確定的是,軟體機(jī)械臂質(zhì)心位置矢量 β∈R3×3必然是驅(qū)動(dòng)空間q的函數(shù),且具有如下表達(dá)式

上式對(duì)時(shí)間求導(dǎo)可得質(zhì)心的速度為

式中,Jβ表示速度雅可比矩陣且具有如下表達(dá)式

因此基于質(zhì)心的動(dòng)能為

同理,基于質(zhì)心的重力勢(shì)能為

對(duì)比式(18)與式(29),不難發(fā)現(xiàn)基于質(zhì)心集中質(zhì)量模型所求軟體機(jī)械臂的重力勢(shì)能和采用連續(xù)分布質(zhì)量所求的動(dòng)力學(xué)模型相等.進(jìn)一步分析基于質(zhì)心所求的軟體機(jī)械臂的動(dòng)能 Kβ與采用連續(xù)分布質(zhì)量積分所求的動(dòng)能 K,顯然,兩種方式計(jì)算出的動(dòng)能都是驅(qū)動(dòng)空間變量q以及其導(dǎo)數(shù)的函數(shù),但驅(qū)動(dòng)空間受驅(qū)動(dòng)力的影響,其隨時(shí)間變化的規(guī)律是不確定的,將影響的變化,盡管其變化規(guī)律不定,但受到物理極限的限制,兩者的變化必定位于某一范圍內(nèi).為了比較 Kβ與 K 的數(shù)量關(guān)系,現(xiàn)假定兩者的變化范圍為[0,30 mm],[-50,50 mm/s],在該范圍內(nèi)隨機(jī)生成10 000 個(gè)數(shù)據(jù)點(diǎn)來(lái)比較兩者的數(shù)量關(guān)系,結(jié)果如圖13 所示.

圖13 隨機(jī)樣本下動(dòng)能大小Fig.13 Kinetic energy under a random sample

顯然兩者的大小并不相等.進(jìn)一步研究?jī)烧叩谋戎店P(guān)系,即對(duì)每一個(gè)隨機(jī)數(shù)據(jù)點(diǎn)計(jì)算,計(jì)算結(jié)果如圖14 所示,由圖發(fā)現(xiàn)比值的范圍落在區(qū)間 [0.5,0.75] 內(nèi),并且大致以0.56 為中心上下波動(dòng).將區(qū)間 [0.5,0.75] 以0.05 為間隔劃分為五個(gè)子區(qū)間,對(duì)隨機(jī)生成的10 000 個(gè)數(shù)據(jù)點(diǎn)落在各區(qū)間的數(shù)量進(jìn)行統(tǒng)計(jì),統(tǒng)計(jì)結(jié)果見(jiàn)表1.

表1 樣本點(diǎn)落入各區(qū)間的個(gè)數(shù)Table 1 The number of sample points falling into each interval

圖14 動(dòng)能等效系數(shù)Fig.14 Kinetic energy equivalent coefficient

計(jì)算落在每個(gè)子區(qū)間的概率,繪制成如圖15 所示的條形圖,由圖15 可知,落在0.5 到0.6 之間的概率約為 93% ,均值記為 ξ,計(jì)算結(jié)果為 ξ=0.56 .

圖15 概率條形圖Fig.15 Probability bar

接著,以積分動(dòng)能 K 為橫坐標(biāo),以質(zhì)心動(dòng)能Kβ為縱坐標(biāo)繪制散點(diǎn)圖,如圖16 所示,顯然散點(diǎn)數(shù)據(jù)大致散落在某直線附近,然后對(duì)數(shù)據(jù)進(jìn)行最小二乘擬合,擬合結(jié)果為一條過(guò)原點(diǎn)的直線,從上面的分析可以看出該直線斜率的物理意義即為 ξ .

圖16 積分動(dòng)能與質(zhì)心動(dòng)能比值關(guān)系擬合曲線Fig.16 Fitting curve of the ratio between integral kinetic energy and center-of-mass kinetic energy

由此可見(jiàn)將軟體機(jī)械臂的質(zhì)量等效為位于質(zhì)心處的集中質(zhì)量所計(jì)算得到的質(zhì)心動(dòng)能 Kβ與實(shí)際的連續(xù)分布質(zhì)量所計(jì)算得到的積分動(dòng)能 K 是不相等的,但是二者的比值近似為一定值 ξ,將該定值稱為動(dòng)能等效系數(shù),為了補(bǔ)償二者的不等,引入基于質(zhì)心的等效質(zhì)量mβ,并且使mβ=m/ξ,如此,便可近似認(rèn)為 Kβ≈K,從而使得基于質(zhì)心所得的動(dòng)能在能量上與分布式積分方法所得動(dòng)能匹配起來(lái).

綜上分析,基于質(zhì)心集中質(zhì)量模型的軟體機(jī)械臂拉格朗日函數(shù)為 Lβ=Kβ-Pβ,其中 Kβ由式(28)得到,.最終,將軟體機(jī)械臂的質(zhì)量等效為位于質(zhì)心的集中質(zhì)量模型的動(dòng)力學(xué)方程以矩陣形式可以寫(xiě)為

3 數(shù)值仿真算例

本小節(jié)對(duì)上述建立的軟體機(jī)械臂動(dòng)力學(xué)模型進(jìn)行數(shù)值仿真分析,動(dòng)力學(xué)仿真參數(shù)參考文獻(xiàn)[32],如表2 所示.

表2 動(dòng)力學(xué)仿真數(shù)值Table 2 Parameters used in the dynamic simulation

3.1 純伸長(zhǎng)模擬

圖17 驅(qū)動(dòng)器長(zhǎng)度變化曲線Fig.17 Actuator length

圖18 配置空間變量變化曲線Fig.18 Configure space variables

3.2 計(jì)算效率比較

為了定量比較動(dòng)力學(xué)建模過(guò)程使用的兩種方法(即采用分布質(zhì)量建立的積分動(dòng)力學(xué)模型以及將質(zhì)量等效成位于軟體機(jī)械臂質(zhì)心的集中質(zhì)量動(dòng)力學(xué)模型)的數(shù)值計(jì)算結(jié)果與計(jì)算效率,對(duì)相同輸入條件下兩種動(dòng)力學(xué)模型的求解結(jié)果以及數(shù)值仿真時(shí)間進(jìn)行比較.圖19 表示的是采用分布質(zhì)量的積分動(dòng)力學(xué)模型求得的結(jié)果,圖20 是采用基于質(zhì)心的集中質(zhì)量動(dòng)力學(xué)模型所求結(jié)果,比較圖19 和圖20,可以發(fā)現(xiàn)二者的求解結(jié)果一致,但使用兩種方法所用的時(shí)間卻大不相同,表3 給出了相同條件下兩種動(dòng)力學(xué)模型的計(jì)算時(shí)間,結(jié)果表明基于質(zhì)心的集中質(zhì)量模型的求解效率比分布質(zhì)量模型提高了大約8 倍.

表3 兩種動(dòng)力學(xué)模型的計(jì)算效率比較Table 3 Comparison of computational efficiency of the two models

圖19 分布質(zhì)量模型響應(yīng)曲線Fig.19 Responses from distributed mass model

圖20 基于質(zhì)心的集中質(zhì)量模型響應(yīng)曲線Fig.20 Responses from concentrated mass model based on mass center

分析其原因在于本文所提出的基于質(zhì)心集中質(zhì)量動(dòng)力學(xué)模型的主要優(yōu)點(diǎn)是它不需要可積項(xiàng)乘積的預(yù)積分.例如,式(14)必須預(yù)先通過(guò)對(duì)中性軸上任意點(diǎn)對(duì)位置求導(dǎo)得到速度雅可比矩陣Jυ,而該矩陣是含有積分變量η的函數(shù)矩陣,這必然會(huì)導(dǎo)致(Jυ)TJυ是關(guān)于積分變量η更加復(fù)雜的矩陣,然后對(duì)η積分后才能得到動(dòng)能.與之相反的是,使用基于質(zhì)心的集中質(zhì)量模型恰好避免了這一復(fù)雜的積分計(jì)算,這是因?yàn)槭?25)是預(yù)先對(duì)積分變量η進(jìn)行積分,顯然這一步的積分運(yùn)算是簡(jiǎn)單的,其后再求速度并計(jì)算動(dòng)能,顯著降低了計(jì)算的復(fù)雜性.

4 結(jié)論

本文以一類具有常曲率彎曲與伸長(zhǎng)能力的軟體機(jī)械臂為研究對(duì)象,對(duì)其空間運(yùn)動(dòng)學(xué)及動(dòng)力學(xué)建模與仿真進(jìn)行了研究,主要研究結(jié)論如下.

(1) 本文采用模態(tài)方法對(duì)軟體機(jī)械臂進(jìn)行運(yùn)動(dòng)學(xué)描述,該方法以多項(xiàng)式函數(shù)描述軟體機(jī)械臂的變形場(chǎng),能夠簡(jiǎn)便地推導(dǎo)出精確、無(wú)奇異且唯一的模態(tài)變換矩陣,該方法在沒(méi)有使用中間變換或曲線參數(shù)的情況下,直接以驅(qū)動(dòng)關(guān)節(jié)變量推導(dǎo)軟體機(jī)械臂運(yùn)動(dòng)學(xué)方程.仿真結(jié)果表明該方法既能有效模擬軟體機(jī)械臂的彎曲行為,又可以模擬軟體機(jī)械臂的純伸長(zhǎng)運(yùn)動(dòng),不需要將純伸長(zhǎng)運(yùn)動(dòng)狀態(tài)單獨(dú)考慮,有效克服了以往模型假定軟體機(jī)械臂變形后曲率不為零(只能模擬彎曲變形,而不能模擬直線運(yùn)動(dòng))的局限性.

(2) 軟體機(jī)械臂動(dòng)力學(xué)建模時(shí),轉(zhuǎn)動(dòng)動(dòng)能計(jì)算復(fù)雜,影響動(dòng)力學(xué)方程的求解效率,但與平動(dòng)動(dòng)能相比,在給定條件下,轉(zhuǎn)動(dòng)動(dòng)能占總動(dòng)能的百分比不超過(guò)3%,其對(duì)動(dòng)力學(xué)結(jié)果影響很小.因此在動(dòng)力學(xué)建模過(guò)程中可以忽略,進(jìn)而提高動(dòng)力學(xué)方程求解效率,為軟體機(jī)械臂的實(shí)時(shí)動(dòng)態(tài)控制奠定了理論基礎(chǔ).

(3) 本文提出的基于質(zhì)心集中質(zhì)量動(dòng)力學(xué)模型,將軟體機(jī)械臂的連續(xù)分布質(zhì)量模型等效為位于質(zhì)心的集中質(zhì)量模型,利用統(tǒng)計(jì)的方法計(jì)算出動(dòng)能等效系數(shù),進(jìn)而通過(guò)動(dòng)能等效系數(shù)實(shí)現(xiàn)兩種模型的動(dòng)能匹配.仿真結(jié)果表明,本文所建立的模型既具有分布質(zhì)量模型的準(zhǔn)確性,又具有集中質(zhì)量模型的計(jì)算高效性,能夠準(zhǔn)確、高效地獲得軟體機(jī)械臂的動(dòng)力學(xué)特性.