陳修龍 崔夢強(qiáng)

(山東科技大學(xué)機(jī)械電子工程學(xué)院， 青島 266590)

0 引言

并聯(lián)機(jī)構(gòu)具有工作精度高、末端執(zhí)行器靈敏性好、負(fù)載能力強(qiáng)等優(yōu)點(diǎn)，故而廣泛應(yīng)用于精密制造業(yè)和空間機(jī)械手等領(lǐng)域[1-5]。隨著機(jī)器人技術(shù)的不斷進(jìn)步，對其設(shè)備創(chuàng)新性提出了更高的要求。具備結(jié)構(gòu)冗余的空間并聯(lián)機(jī)構(gòu)可以較好地提高其機(jī)構(gòu)剛度[6]，但在其高速運(yùn)動時，并聯(lián)機(jī)構(gòu)中支鏈的彈性變形會產(chǎn)生系統(tǒng)振動、內(nèi)應(yīng)力增大和運(yùn)動精度降低等不良影響[7-9]。

近年來，眾多學(xué)者采用不同建模方法對并聯(lián)機(jī)構(gòu)的動力學(xué)性能進(jìn)行了研究，柔性桿件的彈性變形問題對機(jī)構(gòu)動力學(xué)的影響已成為機(jī)構(gòu)學(xué)的研究熱點(diǎn)[10-14]。王宗平等[15]提出具有六自由度的6PSS并聯(lián)機(jī)構(gòu)，基于拉格朗日方程并利用差商求導(dǎo)方式建立了機(jī)構(gòu)的剛體動力學(xué)模型，并分析了運(yùn)動學(xué)和動力學(xué)中各參數(shù)的變化。陳修龍等[16]基于牛頓歐拉法建立了4-UPS-UPU空間機(jī)構(gòu)的動力學(xué)模型，分析了相關(guān)支鏈的受力情況，最后通過虛擬樣機(jī)驗(yàn)證了Matlab理論結(jié)果的正確性。ZHANG等[17]基于有限元描述了柔性桿件梁單元模型，并通過拉格朗日法建立了考慮柔性桿件的動力學(xué)方程，采用一次迭代的方式求解其方程，并研究了彈性振動對機(jī)構(gòu)的影響。WANG等[18]根據(jù)并聯(lián)機(jī)構(gòu)剛?cè)狁詈系奶攸c(diǎn)，基于有限元法通過離散柔性構(gòu)件得到梁單元模型，然后利用KED法對3PRR平面并聯(lián)機(jī)構(gòu)的動力學(xué)方程進(jìn)行組合裝配，研究了其末端構(gòu)件的運(yùn)動特性。YU等[19]利用ANSYS Workbench軟件對一種八連桿機(jī)構(gòu)進(jìn)行了剛?cè)狁詈蟿恿W(xué)分析，研究了柔性構(gòu)件對機(jī)構(gòu)傳動特性曲線的影響，對軸孔尺寸進(jìn)行優(yōu)化發(fā)現(xiàn)，主連桿處的最大等效應(yīng)力明顯降低。LIU等[20]以空間三連桿焊接機(jī)器人為研究對象，基于拉格朗日法建立了剛?cè)狁詈蟿恿W(xué)模型，分析了柔性關(guān)節(jié)處不同剛度系數(shù)對焊接機(jī)器人運(yùn)動精度的影響，結(jié)果表明，系統(tǒng)累積誤差與運(yùn)動軌跡有關(guān)。YAN等[21]提出一種用于臂架大范圍回轉(zhuǎn)運(yùn)動和彈性振動的混合坐標(biāo)模型，運(yùn)用初應(yīng)力法建立了鋼絲繩的剛度矩陣，對超長柔性臂的動力學(xué)研究表明，隨著載荷的增大，其偏航角的周期和幅值也相應(yīng)增大。CHEN等[22]結(jié)合ANSYS和ADAMS對2PRR平面并聯(lián)機(jī)構(gòu)進(jìn)行剛?cè)狁詈蟿恿W(xué)仿真，結(jié)果驗(yàn)證了模型運(yùn)動學(xué)和動力學(xué)正逆解的正確性，最后還研究了柔性桿件對從、動滑塊運(yùn)動特性的影響。姜振海等[23]采用絕對節(jié)點(diǎn)坐標(biāo)方法(ANCF)推導(dǎo)了平面3RRR并聯(lián)機(jī)構(gòu)考慮柔性桿件的動力學(xué)方程，利用直接積分方法(BDF)對方程進(jìn)行求解，發(fā)現(xiàn)柔性構(gòu)件的彈性變形對動平臺影響很大，同時還發(fā)現(xiàn)主動桿件比從動桿變形嚴(yán)重。XING等[24]研究了6PSS并聯(lián)機(jī)構(gòu)剛?cè)狁詈系膭恿W(xué)特性，結(jié)合ADAMS和ANSYS對機(jī)構(gòu)進(jìn)行運(yùn)動仿真求解，將得到的響應(yīng)曲線和機(jī)構(gòu)動力學(xué)模型結(jié)果進(jìn)行相互驗(yàn)證，證明了模型的正確性。

上述彈性動力學(xué)的研究對象大多數(shù)為串聯(lián)機(jī)構(gòu)和平面并聯(lián)機(jī)構(gòu)，而對于具有結(jié)構(gòu)冗余的空間并聯(lián)機(jī)構(gòu)的建模分析研究較少。本文以具有結(jié)構(gòu)冗余3-RRPaR空間并聯(lián)機(jī)構(gòu)為研究對象，基于絕對節(jié)點(diǎn)坐標(biāo)法建立柔性梁單元模型，利用拉格朗日乘子方法推導(dǎo)機(jī)構(gòu)的剛?cè)狁詈蟿恿W(xué)方程，并通過廣義α方法求解其動力學(xué)方程，以分析構(gòu)件材料的彈性模量對機(jī)構(gòu)動力學(xué)的影響。

1 基于ANCF方法的梁單元模型

1.1 ANCF梁單元空間模型

采用三維二節(jié)點(diǎn)的梁單元模型對柔性桿件進(jìn)行建模，其中梁單元模型如圖1所示，圖中OXYZ為系統(tǒng)坐標(biāo)系；梁單元的兩端節(jié)點(diǎn)分別為點(diǎn)i和點(diǎn)j；ri和rj分別為單元兩節(jié)點(diǎn)在OXYZ下的位置矢量。

該模型通過插值法可得出單元上隨機(jī)一點(diǎn)的位置坐標(biāo)(位移場)為

(1)

(2)

其中

(α=i,j;β=x,y,z)S=[s1s2…s8]?I3

(3)

s1=1-3ζ2+2ζ3s2=l(ζ-2ζ2+ζ3)s3=l(ψ-ζψ)s4=l(ξ-ζξ)s5=3ζ2-2ζ3s6=l(-ζ2+ζ3)s7=lζψs8=lζξζ=x/lψ=y/lξ=z/l

式中mi、ni、ci——插值系數(shù)，i=0,1,…,7

S——梁單元形函數(shù)

e——梁單元節(jié)點(diǎn)坐標(biāo)

?——克羅內(nèi)克積運(yùn)算符

I3——三階單位矩陣

l——梁單元變形前長度

從式(2)可看出，梁單元的廣義坐標(biāo)個數(shù)共有24個，分別由節(jié)點(diǎn)處的位置坐標(biāo)和描述3個方向的坐標(biāo)矢量組成。

1.2 單元質(zhì)量矩陣

梁單元質(zhì)量矩陣可通過其動能求出，在計算梁單元動能前需求出其單元任一點(diǎn)的速度，計算式為

(4)

其動能計算式為

(5)

(6)

式中Ve——單元體積

ρ——桿件材料密度

Me——單元質(zhì)量矩陣

1.3 單元廣義彈性力

如圖2所示，當(dāng)點(diǎn)P0通過一段距離u后到點(diǎn)P，根據(jù)連續(xù)介質(zhì)力學(xué)[25]相關(guān)知識可得到變形梯度J表達(dá)式為

(7)

其中

r0=Se0

式中r0——變形體運(yùn)動前點(diǎn)P0的位置矢量

e0——梁單元在初始狀態(tài)下的位置坐標(biāo)向量

J0——單位矩陣

v——單元坐標(biāo)系下的位置矢量

圖3所示單元坐標(biāo)系與系統(tǒng)坐標(biāo)系重合時，可得到用右柯西-格林表達(dá)的變形張量為

C=JTJ

(8)

可推導(dǎo)得拉格朗日應(yīng)變張量為

(9)

其中

式中Sα——形函數(shù)矩陣的第α行元素，α=1,2,3

為了便于后續(xù)積分計算，將式(9)應(yīng)變張量轉(zhuǎn)換成向量表達(dá)形式為

εk=[ε11ε22ε332ε122ε132ε23]T

(10)

通過式(10)推導(dǎo)出梁單元的彈性變形能為

(11)

式中H——構(gòu)件材料的彈性常數(shù)矩陣

對單元變形能求導(dǎo)可得到單元彈性力為

(12)

其中

式中K(e)——單元彈性力剛度矩陣，非常數(shù)矩陣，具有時變性

1.4 單元廣義外力

外力F作用在單元上任意一點(diǎn)時，可通過虛功原理求得該點(diǎn)的力F對單元所做的虛功為

(13)

通過式(13)可得系統(tǒng)坐標(biāo)系OXYZ下沿X軸正方向的重力的一般表達(dá)形式為

(14)

式中me——梁單元質(zhì)量

2 剛?cè)狁詈蟿恿W(xué)模型

2.1 自由度分析

如圖4所示，該機(jī)構(gòu)由定平臺、3條型式相同的支鏈和動平臺組成，其中每條支鏈包括主動臂和具有平行四邊形構(gòu)型的從動臂，可以看出該機(jī)構(gòu)具有結(jié)構(gòu)冗余，有利于增大機(jī)構(gòu)運(yùn)作過程中構(gòu)件的剛度，提高穩(wěn)定性。選擇修正的Grubler-Kutzbach公式[26]計算自由度P，計算式為

(15)

式中n——機(jī)構(gòu)構(gòu)件總數(shù)，取17

m——機(jī)構(gòu)運(yùn)動副總數(shù)，取21

fi——第i個運(yùn)動副的自由度

s——冗余約束個數(shù)，取12

將各個參數(shù)數(shù)值代入式(15)中，可求得本機(jī)構(gòu)自由度為3，且均是移動方向。

2.2 約束方程

多體系統(tǒng)的約束方程需要根據(jù)運(yùn)動副種類的不同建立相關(guān)的位置及方向約束，每個運(yùn)動副相鄰的兩構(gòu)件需要在其自身建立局部坐標(biāo)系以便描述其在空間的位置及姿態(tài)。由于本機(jī)構(gòu)運(yùn)動副均為轉(zhuǎn)動副，故剛性構(gòu)件用自然坐標(biāo)法(NCF)描述，柔性構(gòu)件采用ANCF描述，其機(jī)構(gòu)簡圖如圖5所示，以此確定每個構(gòu)件在系統(tǒng)坐標(biāo)系下的位姿。

該并聯(lián)機(jī)構(gòu)由主動臂作為驅(qū)動構(gòu)件，為了便于約束方程的建立及后續(xù)的計算仿真，主動臂采用質(zhì)心坐標(biāo)法描述，其中一條由定平臺上一點(diǎn)到動平臺上一點(diǎn)組成的支鏈的坐標(biāo)命名，如圖6所示，可知機(jī)構(gòu)廣義坐標(biāo)為

q=(qN,qe)

(16)

式中qN=(Oi1,Oi2,Oi3,Oi4,Oi5,O′i) (i=1,2,3)，主動臂坐標(biāo)包含6個廣義坐標(biāo)，其他均為3個廣義坐標(biāo)；qe=(Oi6,Oi7,Oi8,Oi9)(i=1,2,3)，每個節(jié)點(diǎn)的廣義坐標(biāo)個數(shù)均為12，故該機(jī)構(gòu)一共有207個廣義坐標(biāo)，但為了方便描述廣義坐標(biāo)與自由度的關(guān)系，取節(jié)點(diǎn)的前3個位置坐標(biāo)作為約束廣義坐標(biāo)，可以看出該并聯(lián)機(jī)構(gòu)一共有99個約束廣義坐標(biāo)。

本機(jī)構(gòu)共有3個移動自由度，因此建立線性無關(guān)的約束方程96個，以及3個驅(qū)動約束。Oi(ai,bi,ci)(i=1,2,3)是定平臺上的3個固定點(diǎn)，主動臂處的轉(zhuǎn)動副約束為

(17)

(18)

(19)

其中

式中Ri——第i個支鏈主動臂變換矩陣

li1——第i個支鏈主動臂桿長

O——零矩陣

對于剛性構(gòu)件與柔性構(gòu)件的約束，可采用對應(yīng)節(jié)點(diǎn)重合的形式進(jìn)行約束，即

(20)

Φ5i=[Oi3-Oi2]T(Ai[0,0,1]T)=O1×1

(21)

(22)

(23)

式中Ai——柔性桿件的變換矩陣[27]

根據(jù)動平臺的幾何特性利用NCF法建立三角板單元，其約束方程為

(24)

(25)

該并聯(lián)機(jī)構(gòu)通過運(yùn)動學(xué)反解可求出每條支鏈大臂的運(yùn)動函數(shù)，可建立其驅(qū)動方程的約束為

(26)

式中θi——關(guān)于時間的驅(qū)動函數(shù)，i=1,2,3

由式(17)～(26)可組成3-RRPaR并聯(lián)機(jī)構(gòu)剛?cè)狁詈系募s束方程為

(27)

2.3 質(zhì)量矩陣

質(zhì)量矩陣是反映構(gòu)件在空間中質(zhì)量分布的最佳方法，是構(gòu)建機(jī)構(gòu)動力學(xué)方程不可或缺的一部分。質(zhì)量矩陣為常數(shù)矩陣時可提高仿真計算效率，故該并聯(lián)機(jī)構(gòu)采用NCF和ANCF方法對構(gòu)件進(jìn)行描述，并且柔性構(gòu)件的常數(shù)質(zhì)量矩陣沒有離心力和科氏力。

對于主動臂采用常見的質(zhì)心坐標(biāo)方法進(jìn)行描述，故其質(zhì)量矩陣為

Mk=diag(mk,mk,mk,Ikx,Iky,Ikz) (k=1,6,11)

(28)

對于NCF方法描述的剛性桿件，采用兩點(diǎn)零矢形式通過慣性力的虛功率可求得剛性桿單元質(zhì)量矩陣為

(29)

剛性三角板單元質(zhì)量矩陣為

(30)

對于ANCF方法描述的柔性梁單元，由式(6)可求得其對應(yīng)的質(zhì)量矩陣為

(31)

由式(28)～(31)可組建機(jī)構(gòu)的質(zhì)量矩陣為

MΣ=diag(M1,M2,…,M16)

(32)

2.4 動力學(xué)方程

選擇多體系統(tǒng)動力學(xué)中指標(biāo)-3的DAEs描述動力學(xué)方程，基于拉格朗日乘子法構(gòu)成的3-RRPaR空間并聯(lián)機(jī)構(gòu)剛?cè)狁詈蟿恿W(xué)方程為

(33)

式中Φq——約束方程Φ(q,t)關(guān)于廣義坐標(biāo)的導(dǎo)數(shù)

λ——拉格朗日乘子

QΣ——系統(tǒng)的廣義外力

FΣ——式(12)組裝的廣義彈性力

由于式(33)DAEs不能直接進(jìn)行求解，需要經(jīng)過廣義α方法進(jìn)行差分直接離散成代數(shù)方程的形式求解，具體方法見文獻(xiàn)[28]，將式(33)構(gòu)建成矩陣的形式為

(34)

對式(34)進(jìn)行求導(dǎo)可得其對應(yīng)的Jacobian矩陣為

(35)

3 數(shù)值仿真與分析

3.1 仿真參數(shù)

利用三維軟件建立3-RRPaR空間并聯(lián)機(jī)構(gòu)模型，經(jīng)測量可得構(gòu)件材料相關(guān)參數(shù)如表1～4所示。

表1 3-RRPaR空間并聯(lián)機(jī)構(gòu)主動臂的轉(zhuǎn)動慣量Tab.1 Moment of inertia of active arm of 3-RRPaR spatial parallel mechanism kg·m2

表2 構(gòu)件結(jié)構(gòu)參數(shù)Tab.2 Structural parameters of components

表3 定、動平臺坐標(biāo)參數(shù)Tab.3 Coordinate parameters of fixed and moving platforms m

表4 柔性桿件材料參數(shù)Tab.4 Material parameters of flexible rod

3.2 仿真結(jié)果與分析

并聯(lián)機(jī)構(gòu)擁有3個平移方向的自由度，為了便于分析，規(guī)劃動平臺為一個圓形軌跡，其軌跡函數(shù)表達(dá)式為

(36)

利用Matlab對其進(jìn)行多體系統(tǒng)剛?cè)狁詈蟿恿W(xué)建模，并通過廣義α方法進(jìn)行數(shù)值求解，其中求解步長為0.002 s，譜半徑是0.7，為了避免泊松閉鎖其泊松比取值0[28]。

不同彈性模量E下動平臺的位移響應(yīng)曲線如圖8所示。由8a、8c可知，3種不同彈性模量與剛體機(jī)構(gòu)的位移曲線基本重合。從圖8b、8d可以發(fā)現(xiàn)，彈性模量為116 GPa時與理想情況下的位移誤差曲線最為接近，由圖8d可得，最大位移誤差僅為-2.2×10-7m，在彈性模量為10 GPa時，最大誤差為-6.615×10-6m，彈性模量為1 GPa時最大誤差為-9.579×10-6m。綜上分析發(fā)現(xiàn)彈性模量與理想時的位移誤差基本在數(shù)量級10-5m以內(nèi)，其彈性模量對動平臺位移的影響很小，體現(xiàn)出該具有結(jié)構(gòu)冗余的并聯(lián)機(jī)構(gòu)有著較高的精度。

圖9為3種彈性模量下的速度響應(yīng)曲線。由圖9a、9c可以看出，速度曲線一開始有微小波動，但后面的曲線基本重合，其波動原因可能是柔性桿件突然變形導(dǎo)致的數(shù)值解不收斂，并且Z軸方向與Y軸方向相比波動較大。從圖9b、9d發(fā)現(xiàn)，隨著彈性模量的減小速度誤差變大；由圖9d可知，在彈性模量為1 GPa時幅值的跨距最大，且最大誤差達(dá)到-6.098×10-4m/s，彈性模量為10 GPa時最大誤差為3.92×10-4m/s，彈性模量為116 GPa時最大誤差為2.158×10-4m/s，但是3種誤差基本在數(shù)量級10-3m/s以內(nèi)，可得彈性模量對動平臺速度的影響大于位移的影響。

圖10為3種彈性模量下的加速度響應(yīng)曲線。由圖10a、10c可得，其曲線一開始有著較大波動，Y軸和Z軸兩方向的最大波動分別達(dá)到了-34.78 m/s2和-248.3 m/s2，但是在經(jīng)歷0.1 s后曲線基本重合；由圖10b、10d可知，隨著彈性模量的減小其加速度誤差變大；由圖10d可得，彈性模量為1 GPa時誤差最大可達(dá)到0.109 1 m/s2，彈性模量為10 GPa時最大誤差為0.035 67 m/s2，彈性模量為116 GPa時最大誤差為0.014 72 m/s2。綜合對動平臺的位移、速度、加速度曲線分析可發(fā)現(xiàn)，考慮柔性構(gòu)件的運(yùn)動特性曲線與剛體機(jī)構(gòu)的動力學(xué)特性曲線基本一致，通過誤差發(fā)現(xiàn)隨著彈性模量的減小其位移、速度和加速度曲線偏離理想狀態(tài)的程度越大，且彈性模量對加速度的影響最大，速度次之，位移最小。

采用第四強(qiáng)度理論的等效應(yīng)力進(jìn)行分析，其中在不同彈性模量的影響下該機(jī)構(gòu)各支鏈處桿件的等效應(yīng)力如圖11所示。第1支鏈開始時，E=116 GPa的等效應(yīng)力達(dá)到8.926 MPa，但之后隨著彈性模量減少其等效應(yīng)力曲線波動范圍增大，并且不同曲線之間的差距隨著時間的增大逐漸變小。由圖11b可知，第2支鏈在E=1 GPa時等效應(yīng)力最大達(dá)到了2.446 MPa，并且彈性模量越小，對應(yīng)桿件的等效應(yīng)力曲線波動范圍也就越大。由圖11c可知，第3支鏈在E=116 GPa時等效應(yīng)力最大達(dá)到了15.01 MPa，其3種彈性模量對應(yīng)的曲線規(guī)律與圖11a類似，發(fā)現(xiàn)由于第1支鏈與第3支鏈的對稱性，其兩支鏈的等效應(yīng)力曲線具有明顯的相似性。綜上分析發(fā)現(xiàn)其柔性桿件的彈性模量越小，對應(yīng)等效應(yīng)力曲線的波動范圍也就越大，同時在曲線穩(wěn)定時支鏈處3種彈性模量對應(yīng)的等效應(yīng)力基本在數(shù)量級1 MPa以內(nèi)。

4 結(jié)論

(1)基于絕對節(jié)點(diǎn)坐標(biāo)法描述了三維二節(jié)點(diǎn)柔性梁單元空間模型，并根據(jù)功能關(guān)系推導(dǎo)了梁單元的質(zhì)量矩陣，根據(jù)連續(xù)介質(zhì)力學(xué)推導(dǎo)了彈性力矩陣，最后利用虛功原理推導(dǎo)得出重力方向上單元廣義外力的表達(dá)式。

(2)分析了3-RRPaR空間并聯(lián)機(jī)構(gòu)的結(jié)構(gòu)特征，采用自然坐標(biāo)法和絕對節(jié)點(diǎn)坐標(biāo)法建立了機(jī)構(gòu)約束方程和具有常數(shù)的質(zhì)量矩陣，最后基于拉格朗日乘子法推導(dǎo)得出機(jī)構(gòu)的剛?cè)狁詈蟿恿W(xué)方程，并在Matlab中對方程進(jìn)行數(shù)值求解。

(3)得到不同彈性模量下的動力學(xué)響應(yīng)曲線，3種模量下的運(yùn)動曲線與理想軌跡基本重合，驗(yàn)證了所建模型的正確性，結(jié)果表明，柔性構(gòu)件的彈性模量越小，對機(jī)構(gòu)運(yùn)動精度的影響越大，且彈性模量對機(jī)構(gòu)動平臺加速度的影響最大，對位移的影響最小，體現(xiàn)出具有結(jié)構(gòu)冗余的空間并聯(lián)機(jī)構(gòu)精度高、穩(wěn)定性好的特點(diǎn)。給出了機(jī)構(gòu)在高速運(yùn)動過程中等效應(yīng)力的變化情況，并對比分析了3種不同彈性模量對等效應(yīng)力的影響。