基于絕對節點坐標法的空間并聯機構剛柔耦合動力學建模

陳修龍 崔夢強

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

0 引言

并聯機構具有工作精度高、末端執行器靈敏性好、負載能力強等優點，故而廣泛應用于精密制造業和空間機械手等領域[1-5]。隨著機器人技術的不斷進步，對其設備創新性提出了更高的要求。具備結構冗余的空間并聯機構可以較好地提高其機構剛度[6]，但在其高速運動時，并聯機構中支鏈的彈性變形會產生系統振動、內應力增大和運動精度降低等不良影響[7-9]。

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

上述彈性動力學的研究對象大多數為串聯機構和平面并聯機構，而對于具有結構冗余的空間并聯機構的建模分析研究較少。本文以具有結構冗余3-RRPaR空間并聯機構為研究對象，基于絕對節點坐標法建立柔性梁單元模型，利用拉格朗日乘子方法推導機構的剛柔耦合動力學方程，并通過廣義α方法求解其動力學方程，以分析構件材料的彈性模量對機構動力學的影響。

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

1.1 ANCF梁單元空間模型

采用三維二節點的梁單元模型對柔性桿件進行建模，其中梁單元模型如圖1所示，圖中OXYZ為系統坐標系；梁單元的兩端節點分別為點i和點j；ri和rj分別為單元兩節點在OXYZ下的位置矢量。

該模型通過插值法可得出單元上隨機一點的位置坐標(位移場)為

(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——插值系數，i=0,1,…,7

S——梁單元形函數

e——梁單元節點坐標

?——克羅內克積運算符

I3——三階單位矩陣

l——梁單元變形前長度

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

1.2 單元質量矩陣

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

(4)

其動能計算式為

(5)

(6)

式中Ve——單元體積

ρ——桿件材料密度

Me——單元質量矩陣

1.3 單元廣義彈性力

如圖2所示，當點P0通過一段距離u后到點P，根據連續介質力學[25]相關知識可得到變形梯度J表達式為

(7)

其中

r0=Se0

式中r0——變形體運動前點P0的位置矢量

e0——梁單元在初始狀態下的位置坐標向量

J0——單位矩陣

v——單元坐標系下的位置矢量

圖3所示單元坐標系與系統坐標系重合時，可得到用右柯西-格林表達的變形張量為

C=JTJ

(8)

可推導得拉格朗日應變張量為

(9)

其中

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

為了便于后續積分計算，將式(9)應變張量轉換成向量表達形式為

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

(10)

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

(11)

式中H——構件材料的彈性常數矩陣

對單元變形能求導可得到單元彈性力為

(12)

其中

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

1.4 單元廣義外力

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

(13)

通過式(13)可得系統坐標系OXYZ下沿X軸正方向的重力的一般表達形式為

(14)

式中me——梁單元質量

2 剛柔耦合動力學模型

2.1 自由度分析

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

(15)

式中n——機構構件總數，取17

m——機構運動副總數，取21

fi——第i個運動副的自由度

s——冗余約束個數，取12

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

2.2 約束方程

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

該并聯機構由主動臂作為驅動構件，為了便于約束方程的建立及后續的計算仿真，主動臂采用質心坐標法描述，其中一條由定平臺上一點到動平臺上一點組成的支鏈的坐標命名，如圖6所示，可知機構廣義坐標為

q=(qN,qe)

(16)

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

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

(17)

(18)

(19)

其中

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

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

O——零矩陣

對于剛性構件與柔性構件的約束，可采用對應節點重合的形式進行約束，即

(20)

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

(21)

(22)

(23)

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

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

(24)

(25)

該并聯機構通過運動學反解可求出每條支鏈大臂的運動函數，可建立其驅動方程的約束為

(26)

式中θi——關于時間的驅動函數，i=1,2,3

由式(17)～(26)可組成3-RRPaR并聯機構剛柔耦合的約束方程為

(27)

2.3 質量矩陣

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

對于主動臂采用常見的質心坐標方法進行描述，故其質量矩陣為

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

(28)

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

(29)

剛性三角板單元質量矩陣為

(30)

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

(31)

由式(28)～(31)可組建機構的質量矩陣為

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

(32)

2.4 動力學方程

選擇多體系統動力學中指標-3的DAEs描述動力學方程，基于拉格朗日乘子法構成的3-RRPaR空間并聯機構剛柔耦合動力學方程為

(33)

式中Φq——約束方程Φ(q,t)關于廣義坐標的導數

λ——拉格朗日乘子

QΣ——系統的廣義外力

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

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

(34)

對式(34)進行求導可得其對應的Jacobian矩陣為

(35)

3 數值仿真與分析

3.1 仿真參數

利用三維軟件建立3-RRPaR空間并聯機構模型，經測量可得構件材料相關參數如表1～4所示。

表1 3-RRPaR空間并聯機構主動臂的轉動慣量Tab.1 Moment of inertia of active arm of 3-RRPaR spatial parallel mechanism kg·m2

表2 構件結構參數Tab.2 Structural parameters of components

表3 定、動平臺坐標參數Tab.3 Coordinate parameters of fixed and moving platforms m

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

3.2 仿真結果與分析

并聯機構擁有3個平移方向的自由度，為了便于分析，規劃動平臺為一個圓形軌跡，其軌跡函數表達式為

(36)

利用Matlab對其進行多體系統剛柔耦合動力學建模，并通過廣義α方法進行數值求解，其中求解步長為0.002 s，譜半徑是0.7，為了避免泊松閉鎖其泊松比取值0[28]。

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

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

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

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

4 結論

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

(2)分析了3-RRPaR空間并聯機構的結構特征，采用自然坐標法和絕對節點坐標法建立了機構約束方程和具有常數的質量矩陣，最后基于拉格朗日乘子法推導得出機構的剛柔耦合動力學方程，并在Matlab中對方程進行數值求解。

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