劉 涼 趙新華 周海波 王嘉斌

(天津理工大學(xué)天津市先進(jìn)機(jī)電系統(tǒng)設(shè)計(jì)與智能控制重點(diǎn)實(shí)驗(yàn)室， 天津 300384)

0 引言

柔性多體系統(tǒng)由于建模復(fù)雜和求解困難一直是動(dòng)力學(xué)領(lǐng)域研究的熱點(diǎn)[1-5]。在建模方法上，運(yùn)動(dòng)彈性動(dòng)力學(xué)法的模型簡單，但忽略了大范圍剛體運(yùn)動(dòng)與彈性變形之間的耦合特性，僅將二者進(jìn)行線性疊加，故無法對(duì)彈性體大變形進(jìn)行精確描述[6-8]。浮動(dòng)坐標(biāo)系法將構(gòu)件的彈性變形與大范圍剛體運(yùn)動(dòng)構(gòu)建為零次耦合模型，適用于描述小變形的場合，由于剛?cè)狁詈享?xiàng)的存在，使復(fù)雜柔性機(jī)構(gòu)的求解變得十分困難[9-11]。SHABANA等[12-14]提出了絕對(duì)節(jié)點(diǎn)坐標(biāo)法，它使用絕對(duì)位置坐標(biāo)及其變形梯度作為柔性體單元的節(jié)點(diǎn)坐標(biāo)，可直接描述柔性體大范圍剛體運(yùn)動(dòng)和大彈性變形以及二者的非線性耦合特性，但復(fù)雜的應(yīng)變能形式降低了動(dòng)力系統(tǒng)的求解效率。雖可通過引入稀疏不變矩陣加以改善[15]，但求解時(shí)應(yīng)規(guī)避泊松閉鎖問題。

通常隱式積分法在求解多體系統(tǒng)的剛性微分方程時(shí)更為有效[16]。HUSSEIN等[17]將隱式HHT-I3法與顯式ADAMS法作了比較；由于Generalized-α數(shù)值迭代法可對(duì)高頻模態(tài)實(shí)現(xiàn)可控的人工耗散[18-19]，求解效率較高。目前，非因果解仍然是評(píng)判柔性系統(tǒng)逆動(dòng)力學(xué)解特性的重要標(biāo)準(zhǔn)，因?yàn)樗芮蟪鑫ㄒ挥薪缃鈁20]。但該方法不適于具有高度非線性彈性力項(xiàng)的動(dòng)力系統(tǒng)，應(yīng)構(gòu)建合理的求解方法來尋求多體系統(tǒng)高效穩(wěn)定的數(shù)值因果解，同時(shí)還應(yīng)考慮約束條件下動(dòng)力系統(tǒng)的相容性問題。本文研究3-RRRU空間剛?cè)狁詈喜⒙?lián)機(jī)器人動(dòng)力學(xué)模型穩(wěn)定因果解的求解策略。

1 剛性構(gòu)件與柔性構(gòu)件的質(zhì)量矩陣

由于空間3-RRRU并聯(lián)機(jī)構(gòu)中含有剛性梁單元、剛性三角動(dòng)平臺(tái)和柔性梁單元，所以應(yīng)先分別求出其質(zhì)量矩陣以便構(gòu)建系統(tǒng)的動(dòng)力學(xué)方程。

1.1 空間剛性梁單元的質(zhì)量矩陣

圖1基于自然坐標(biāo)法描述了一種各向同性、質(zhì)量分布均勻且橫截面一致的空間剛性梁單元。梁的中性線上任意一點(diǎn)的絕對(duì)位置矢量可表示為

ri(x)=rji+x(rj+1i-rji)/|rj+1i-rji|=Srqi

(1)

其中

(2)

(3)

式中Sr——?jiǎng)傂粤旱男魏瘮?shù)，Sr∈R3×6

qi——?jiǎng)傂粤旱淖匀蛔鴺?biāo)矢量，qi∈R6

?——克羅內(nèi)科張量積

圖1 NCF法描述的空間剛性梁單元模型Fig.1 Spatial rigid beam described in NCF

根據(jù)慣性力的虛功率計(jì)算公式[20]，剛性梁單元的質(zhì)量矩陣可表示為

(4)

式中ρi、mi——?jiǎng)傂粤好芏?、質(zhì)量

顯然，式(4)描述的質(zhì)量矩陣為對(duì)稱的常數(shù)矩陣。

1.2 剛性三角形平臺(tái)的質(zhì)量矩陣

一種基于自然坐標(biāo)法描述的各向同性、質(zhì)量分布均勻且厚度一致的空間正三角形平臺(tái)如圖2所示，其外接圓的半徑為r。為求其質(zhì)量矩陣，引入兩個(gè)空間基點(diǎn)A41與A42以及在局部坐標(biāo)系oxyz下描述的兩個(gè)單位矢量u和v。這兩個(gè)基點(diǎn)亦可在局部系下分別用矢量ri和rj來描述。三角形平臺(tái)上任意一點(diǎn)在全局坐標(biāo)系下的度量為[20]

r=Cqa

(5)

式中C——三角平臺(tái)形函數(shù)，C∈R3×12

qa——三角平臺(tái)自然坐標(biāo)矢量，qa∈R12

由慣性力的虛功率可得出該平臺(tái)的質(zhì)量矩陣

(6)

式中ρ、mP——三角形平臺(tái)密度、質(zhì)量

該質(zhì)量矩陣為對(duì)稱的常數(shù)矩陣。對(duì)應(yīng)的自然坐標(biāo)矢量為

(7)

式中xi、yi、zi(i=3,6,9)——平臺(tái)頂點(diǎn)的自然坐標(biāo)

圖2 NCF法描述的空間剛性三角形平臺(tái)模型Fig.2 Spatial rigid triangular platform described in NCF

1.3 空間柔性梁單元的質(zhì)量矩陣

利用絕對(duì)節(jié)點(diǎn)坐標(biāo)法[12]，可求出各向同性、質(zhì)量分布均勻且橫截面一致的空間柔性梁單元的質(zhì)量矩陣

Mij=?VijρijSijTSijdVij

(8)

式中ρij、Vij——柔性梁密度、體積

Sij——柔性梁單元形函數(shù)

柔性梁單元的質(zhì)量矩陣為對(duì)稱的常數(shù)矩陣，該單元的節(jié)點(diǎn)坐標(biāo)矢量eij為

(9)

式中m、n——梁單元的首端點(diǎn)、末端點(diǎn)

rijm、rijn——首末端點(diǎn)處絕對(duì)位置矢量

2 剛?cè)狁詈喜⒙?lián)機(jī)器人的動(dòng)力學(xué)方程

2.1 3-RRRU剛?cè)狁詈喜⒙?lián)機(jī)器人

3-RRRU并聯(lián)機(jī)器人如圖3所示，它由靜平臺(tái)、動(dòng)平臺(tái)及3條結(jié)構(gòu)對(duì)稱的支鏈組成，其中靜平臺(tái)和動(dòng)平臺(tái)為等邊三角形，各支鏈包含3個(gè)運(yùn)動(dòng)桿件(前2個(gè)桿件為剛性桿，最后1個(gè)桿件為柔性桿)、3個(gè)轉(zhuǎn)動(dòng)副和1個(gè)虎克鉸。所有運(yùn)動(dòng)關(guān)節(jié)都是剛性關(guān)節(jié)，與靜平臺(tái)相連的轉(zhuǎn)動(dòng)副為驅(qū)動(dòng)副。圖4給出了利用自然坐標(biāo)法和絕對(duì)節(jié)點(diǎn)坐標(biāo)法構(gòu)建動(dòng)力學(xué)方程時(shí)各個(gè)桿件廣義坐標(biāo)的定義方式。

圖3 3-RRRU剛?cè)狁詈喜⒙?lián)機(jī)器人Fig.3 3-RRRU rigid-flexible parallel manipulator1.靜平臺(tái) 2.剛性桿 3.柔性桿 4.動(dòng)平臺(tái) 5.虎克鉸 6.轉(zhuǎn)動(dòng)關(guān)節(jié) 7.驅(qū)動(dòng)關(guān)節(jié)

圖4 并聯(lián)機(jī)器人的廣義坐標(biāo)Fig.4 Generalized coordinates of parallel mechanism

2.2 并聯(lián)機(jī)器人正動(dòng)力學(xué)模型

根據(jù)能量原理和前述的剛性、柔性構(gòu)件的質(zhì)量矩陣可求出系統(tǒng)的拉格朗日動(dòng)力學(xué)方程，由于此時(shí)驅(qū)動(dòng)桿件的動(dòng)力學(xué)參數(shù)是已知的，因此該機(jī)構(gòu)在約束條件下的正動(dòng)力學(xué)模型為

(10)

其中，系統(tǒng)的廣義坐標(biāo)矢量可表示為

(11)

式中M——系統(tǒng)的常質(zhì)量矩陣

Qg——廣義重力矢量

Q′——廣義驅(qū)動(dòng)力矩矢量

Qe——廣義彈性力矢量

Φ——系統(tǒng)約束方程

Φq——系統(tǒng)約束方程的雅可比矩陣

λ——拉格朗日乘子矢量

q——系統(tǒng)的廣義坐標(biāo)矢量

為了便于求解，廣義坐標(biāo)矢量中包含了A11、A12和A13處已知的自然坐標(biāo)xi和zi(i=1，4，7)。

由于正動(dòng)力學(xué)模型在系統(tǒng)約束條件下進(jìn)行求解，因此，首先列出第1支鏈的幾何約束方程

(12)

式(12)中，前2個(gè)約束方程描述了2個(gè)剛性桿的長度約束；第3個(gè)約束方程表明第2個(gè)剛性桿只做平面運(yùn)動(dòng)；其他方程則是對(duì)剛性轉(zhuǎn)動(dòng)關(guān)節(jié)和剛性虎克鉸進(jìn)行約束。同理，另2條支鏈亦可列寫12個(gè)約束方程。最后，剛性動(dòng)平臺(tái)的約束方程為

(13)

因此，正動(dòng)力學(xué)模型的全部約束方程為21個(gè)。

2.3 并聯(lián)機(jī)器人逆動(dòng)力學(xué)模型

并聯(lián)機(jī)器人逆動(dòng)力學(xué)模型與式(10)具有相同的形式，但是由于各支鏈的驅(qū)動(dòng)力矩τ1、τ2和τ3均為未知量，故需要添加3個(gè)新的約束方程才能進(jìn)行求解。約束條件下的逆動(dòng)力學(xué)模型為

(14)

式中ΦT——新添加的3個(gè)約束方程

(15)

式中xp、yp、zp——?jiǎng)悠脚_(tái)中心點(diǎn)P處運(yùn)動(dòng)軌跡

由于其他約束方程與正動(dòng)力學(xué)模型中的完全相同，故逆動(dòng)力學(xué)模型的全部約束方程為24個(gè)。

3 剛?cè)狁詈蟿?dòng)力學(xué)系統(tǒng)的求解策略

3.1 動(dòng)力學(xué)方程的數(shù)值積分法與相容性問題

由于描述系統(tǒng)動(dòng)力學(xué)模型的微分方程具有剛性特性，所以利用Generalized-α法對(duì)高頻分量實(shí)現(xiàn)可控的數(shù)值耗散能準(zhǔn)確地求出系統(tǒng)方程的數(shù)值解。作為一種隱式數(shù)值積分法，在迭代求解時(shí)應(yīng)將系統(tǒng)正、逆動(dòng)力學(xué)方程分別改寫為

(16)

(17)

為保證動(dòng)力學(xué)模型的數(shù)值解滿足系統(tǒng)的求解精度要求，應(yīng)對(duì)式(10)和式(14)中的動(dòng)力學(xué)方程和約束方程的誤差函數(shù)進(jìn)行定義。其中，正動(dòng)力學(xué)方程和約束方程的誤差函數(shù)分別定義為

(18)

(19)

對(duì)于逆動(dòng)力學(xué)系統(tǒng)，其約束方程的誤差函數(shù)還需在式(19)中添加由ΦT描述的3個(gè)約束方程。由于誤差函數(shù)采用了平方和的根式形式，所以對(duì)一個(gè)給定的誤差預(yù)設(shè)值，能夠保證任意一個(gè)迭代方程的計(jì)算誤差都小于該預(yù)設(shè)值。為了綜合考慮該動(dòng)力系統(tǒng)的迭代誤差，將動(dòng)力學(xué)方程和約束方程的誤差函數(shù)之和定義為該系統(tǒng)的誤差函數(shù)。

利用上述方法迭代求解時(shí)，初始階段動(dòng)力學(xué)方程的收斂速度比較快，系統(tǒng)誤差明顯下降，但由于系統(tǒng)方程具有剛性特性，該誤差的下降幅度會(huì)逐漸減小。隨著迭代過程的不斷深入，系統(tǒng)誤差不會(huì)進(jìn)一步降低反而會(huì)不斷升高，產(chǎn)生這一現(xiàn)象的原因是系統(tǒng)約束方程的迭代誤差正在不斷增加，這說明動(dòng)力學(xué)方程在求解過程中與幾何約束方程所保持的相容性原則被破壞。而出現(xiàn)該問題的根源在于系統(tǒng)中剛性構(gòu)件與柔性構(gòu)件的運(yùn)動(dòng)狀態(tài)是在不同時(shí)間尺度下變化的，但二者卻混合在一起進(jìn)行積分，造成了描述剛性構(gòu)件運(yùn)動(dòng)的自然坐標(biāo)與描述柔性構(gòu)件運(yùn)動(dòng)的絕對(duì)節(jié)點(diǎn)坐標(biāo)在迭代時(shí)的更新速率不同，從而破壞了系統(tǒng)動(dòng)力學(xué)方程與約束方程的相容性。

為解決該問題，提出一種瞬態(tài)剛體校正法，它的基本思想是：在迭代求解過程中，將當(dāng)前已產(chǎn)生變形的并聯(lián)機(jī)構(gòu)視為一個(gè)剛性機(jī)構(gòu)，倘若此時(shí)相容性原則被破壞，則必須在該瞬時(shí)狀態(tài)對(duì)該剛性機(jī)構(gòu)的運(yùn)動(dòng)學(xué)模型進(jìn)行重構(gòu)求解，由于此刻的運(yùn)動(dòng)學(xué)模型與剛?cè)狁詈蟿?dòng)力系統(tǒng)的約束方程相容，同時(shí)能保持柔性桿的變形狀態(tài)與彈性力不變，所以該校正法提供的輔助搜索路徑或迭代路徑能夠保持動(dòng)力學(xué)方程與幾何約束方程的相容性。此時(shí)可利用自然坐標(biāo)法描述各個(gè)構(gòu)件的瞬時(shí)運(yùn)動(dòng)狀態(tài)，即先將圖4中A2i(i=1，2，3)點(diǎn)處的絕對(duì)節(jié)點(diǎn)坐標(biāo)依次改為自然坐標(biāo)(x2,y2,z2)、(x5,y5,z5)和(x8,y8,z8)；然后再將A4i(i=1,2,3)點(diǎn)處的絕對(duì)節(jié)點(diǎn)坐標(biāo)依次改為自然坐標(biāo)(x3,y3,z3)、(x6,y6,z6)和(x9,y9,z9)。

3.2 正動(dòng)力學(xué)瞬態(tài)剛體校正法

由于已知驅(qū)動(dòng)桿件的運(yùn)動(dòng)狀態(tài)，故只需對(duì)其他剛性構(gòu)件和被視為剛性構(gòu)件的柔性桿的運(yùn)動(dòng)狀態(tài)進(jìn)行校正，第1支鏈的瞬時(shí)幾何約束方程為

(20)

其中，前2個(gè)方程約束了第2桿和第3桿的端點(diǎn)長度；第3個(gè)方程是柔性桿上A41點(diǎn)到一個(gè)給定平面的距離約束(距離為d1)，該平面通過A21點(diǎn)并與第2轉(zhuǎn)動(dòng)副的軸線相垂直；最后一個(gè)方程約束了虎克鉸2個(gè)轉(zhuǎn)動(dòng)副的軸線，其中u3x、u3y和u3z為虎克鉸第1轉(zhuǎn)動(dòng)副軸線的全局坐標(biāo)。同理，可列出其他2個(gè)支鏈的8個(gè)瞬時(shí)幾何約束方程。最后，三角形剛性動(dòng)平臺(tái)的瞬時(shí)幾何約束方程為

(21)

正動(dòng)力學(xué)瞬態(tài)剛體校正法是綜合利用這15個(gè)約束方程來重構(gòu)該瞬時(shí)剛性機(jī)構(gòu)的運(yùn)動(dòng)學(xué)模型。

3.3 逆動(dòng)力學(xué)瞬態(tài)剛體校正法

此時(shí)已知末端執(zhí)行器的運(yùn)動(dòng)狀態(tài)，需要校正各支鏈構(gòu)件及動(dòng)平臺(tái)的運(yùn)動(dòng)狀態(tài)，第1支鏈的5個(gè)瞬時(shí)幾何約束方程包含式(20)和以下約束方程

(22)

以上方程約束了驅(qū)動(dòng)桿的長度；同理，可列出其他2支鏈的10個(gè)瞬時(shí)幾何約束方程。最后三角形剛性動(dòng)平臺(tái)的6個(gè)瞬時(shí)幾何約束方程包含式(21)和以下約束方程

(23)

逆動(dòng)力學(xué)瞬態(tài)剛體校正法是綜合利用這21個(gè)約束方程來重構(gòu)該瞬時(shí)剛性機(jī)構(gòu)的運(yùn)動(dòng)學(xué)模型。

3.4 剛?cè)狁詈喜⒙?lián)機(jī)機(jī)器人的動(dòng)力學(xué)求解策略

為了獲得剛?cè)狁詈舷到y(tǒng)穩(wěn)定的數(shù)值因果解，除了要滿足相容性條件外，還應(yīng)該從系統(tǒng)能量、求解穩(wěn)定性、初值點(diǎn)的選取等角度來評(píng)價(jià)和篩選動(dòng)力系統(tǒng)的數(shù)值解，求解原則如下：

(1)從理論上講，約束條件下動(dòng)力學(xué)方程的解有無窮多組，故在求解過程中應(yīng)尋求一組能量最優(yōu)解。若其解析形式難以求出，則求解方法必須有能力對(duì)解集進(jìn)行篩選，并最大限度地減少系統(tǒng)耗能，來降低殘余能量在非采樣點(diǎn)處引起的振動(dòng)。

(2) 求解方法應(yīng)保證各個(gè)采樣點(diǎn)處動(dòng)力學(xué)參數(shù)的解具有良好的平滑性與一致性，在降低系統(tǒng)耗能的前提下，根據(jù)具體運(yùn)動(dòng)情況將其變化率控制在合理的范圍內(nèi)，來保證系統(tǒng)能量輸出的可持續(xù)性與穩(wěn)定性，避免出現(xiàn)大幅度的能量波動(dòng)現(xiàn)象。

(3) 由于動(dòng)力學(xué)模型存在剛性特性，為了提高求解效率和計(jì)算精度，應(yīng)對(duì)運(yùn)動(dòng)起始點(diǎn)的初始狀態(tài)(初值點(diǎn))單獨(dú)進(jìn)行迭代求解。若初值點(diǎn)為靜止?fàn)顟B(tài)，則應(yīng)根據(jù)系統(tǒng)的動(dòng)力學(xué)模型推導(dǎo)出相應(yīng)的靜力學(xué)方程，并對(duì)其進(jìn)行迭代求解。

(4) 在求解過程中，針對(duì)動(dòng)力學(xué)方程與約束方程之間存在的相容性問題，應(yīng)合理高效地尋求與之相容的迭代路徑搜索方法，不斷降低系統(tǒng)的迭代誤差，避免求解失敗。

(5) 大變形柔性構(gòu)件的建模缺陷會(huì)導(dǎo)致在求解過程中出現(xiàn)泊松閉鎖問題[21]，可采用縮減積分法解決該問題。另一種方法是建模時(shí)使用新型柔性梁單元模型[22]，來徹底規(guī)避該閉鎖問題的發(fā)生。

4 求解實(shí)例與實(shí)驗(yàn)分析

剛?cè)狁詈?-RRRU并聯(lián)機(jī)器人的仿真物理參數(shù)如表1所示。其靜平臺(tái)和動(dòng)平臺(tái)的外接圓半徑分別為0.175 m和0.06 m。各柔性桿的楊氏彈性模量E為6.9×108Pa，泊松比ν為0.3，截面為圓形，截面半徑為0.01 m。各桿件的質(zhì)量分布均勻。

表1 3-RRRU并聯(lián)機(jī)器人的仿真物理參數(shù)Tab.1 Physical parameters of 3-RRRU parallel robot

4.1 正動(dòng)力學(xué)模型求解實(shí)例

正動(dòng)力學(xué)模型在求解時(shí)需要已知驅(qū)動(dòng)桿的運(yùn)動(dòng)狀態(tài)，為此先按照一個(gè)給定軌跡求出理想的全剛性3-RRRU并聯(lián)機(jī)器人驅(qū)動(dòng)桿的運(yùn)動(dòng)狀態(tài)，然后將其作為當(dāng)前剛?cè)狁詈喜⒙?lián)機(jī)器人正動(dòng)力學(xué)模型的已知量，來求解末端執(zhí)行器的運(yùn)動(dòng)狀態(tài)。為驗(yàn)證求解策略的可行性，將末端執(zhí)行器運(yùn)動(dòng)軌跡設(shè)為平面ZG=-0.8 m上半徑為0.1 m的一個(gè)標(biāo)準(zhǔn)圓周，其額定運(yùn)動(dòng)速度設(shè)為0.5 m/s。為使運(yùn)動(dòng)軌跡更平滑，這里引入S-型加減速機(jī)制。在求解過程中，將動(dòng)力學(xué)方程和約束方程的誤差函數(shù)閾值分別設(shè)為10-6和10-8。圖5為已知的剛性并聯(lián)機(jī)構(gòu)的驅(qū)動(dòng)力矩，此時(shí)求出的末端執(zhí)行器實(shí)際運(yùn)動(dòng)軌跡及其與理想圓周的位置誤差如圖6所示；動(dòng)平臺(tái)法向量與全局坐標(biāo)系ZG軸的夾角如圖7所示。

圖5 正動(dòng)力學(xué)模型的輸入驅(qū)動(dòng)力矩Fig.5 Actuated torques of each limb for forward dynamics

圖6 末端執(zhí)行器的空間運(yùn)動(dòng)軌跡和位置誤差Fig.6 Spatial trajectory and errors of end-effector

圖7 動(dòng)平臺(tái)法向量與ZG軸的夾角Fig.7 Angle between platform normal vector and ZG-axis

從上述結(jié)果可知，末端執(zhí)行器的運(yùn)動(dòng)軌跡不是標(biāo)準(zhǔn)圓周，在x軸方向的軌跡偏差最大，其最大值為8.7×10-3m，而其余2個(gè)方向的偏差均低于2×10-3m，這與坐標(biāo)系的擺放位置和機(jī)構(gòu)的空間位姿有關(guān)。由于柔性桿的變形使得動(dòng)平臺(tái)無法保持平動(dòng)狀態(tài)。所以根據(jù)剛性模型的逆動(dòng)力學(xué)解對(duì)剛?cè)狁詈蠙C(jī)構(gòu)實(shí)施控制無法得到令人滿意的效果。

各支鏈柔性桿中心點(diǎn)處的橫向變形[12]與剪切角如圖8所示。二者均呈現(xiàn)出非線性的變化趨勢，其中前者的變化過程不具有周期對(duì)稱性；而后者的幅值比較小，其原因是柔性桿的縱向尺寸遠(yuǎn)大于其橫向尺寸，這符合Euler-Bernoulli梁的假設(shè)。圖9記錄了約束方程和動(dòng)力學(xué)方程的誤差函數(shù)值，它們均滿足系統(tǒng)求解前預(yù)設(shè)的誤差閾值。這驗(yàn)證了正動(dòng)力學(xué)瞬態(tài)剛體校正法的可行性和有效性。

圖10 末端執(zhí)行器的運(yùn)動(dòng)狀態(tài)Fig.10 States of end-effector

圖8 柔性桿中心點(diǎn)的橫向變形和剪切角Fig.8 Transverse deformation and shear angles of midpoints of flexible links

圖9 迭代求解過程中約束方程誤差和動(dòng)力學(xué)方程誤差Fig.9 Iterative errors of constraint and dynamic equations

4.2 逆動(dòng)力學(xué)模型求解實(shí)例

依然采用正動(dòng)力學(xué)中末端執(zhí)行器的規(guī)劃軌跡和規(guī)劃速度，根據(jù)求解策略求出驅(qū)動(dòng)構(gòu)件和柔性桿的運(yùn)動(dòng)狀態(tài)，通過引入S型加減速機(jī)制來適當(dāng)縮短系統(tǒng)的加減速過程。圖10a、10b給出了末端執(zhí)行器運(yùn)動(dòng)狀態(tài)的迭代求解值與理想值的誤差，其中運(yùn)動(dòng)軌跡和速度與理想值都十分接近，二者的誤差分別低于2×10-11m和3×10-8m/s，這滿足大多數(shù)理論分析與工程應(yīng)用的要求。其次，圖10c中給出的動(dòng)平臺(tái)法向量與ZG軸的夾角呈現(xiàn)出周期性的變化規(guī)律，其最大值為5.3°，出現(xiàn)在起止位置處；勻速運(yùn)動(dòng)段中該夾角低于0.5°，而且變化過程平穩(wěn)，無大幅度振動(dòng)的現(xiàn)象發(fā)生，這說明所求的逆動(dòng)力學(xué)解能通過控制驅(qū)動(dòng)桿的運(yùn)動(dòng)狀態(tài)實(shí)現(xiàn)剛?cè)狁詈蠙C(jī)器人平穩(wěn)的軌跡跟蹤控制。

此時(shí)驅(qū)動(dòng)力矩和驅(qū)動(dòng)角位移如圖11a、11c所示，由于動(dòng)平臺(tái)不再保持平動(dòng)狀態(tài)，所以各驅(qū)動(dòng)力矩之間以及角位移之間的對(duì)稱性已被完全破壞。其中，驅(qū)動(dòng)力矩最大值出現(xiàn)在第2支鏈為-11.1 N·m，最小值出現(xiàn)在第1支鏈為-8.3 N·m；驅(qū)動(dòng)角位移的正向最大值出現(xiàn)在第1支鏈為9.8°。圖11b、11d將驅(qū)動(dòng)力矩和驅(qū)動(dòng)角位移與理想剛性模型的對(duì)應(yīng)值進(jìn)行了比較，其偏差呈現(xiàn)出非線性的變化特征，其中第2支鏈驅(qū)動(dòng)力矩的偏差值最大為-0.39 N·m；而第1支鏈的角度偏差最大為0.84°。各柔性桿中心點(diǎn)處橫向變形與剪切角如圖11e、11f所示，二者均呈現(xiàn)出非常規(guī)范的周期性變化規(guī)律，其中第1支鏈柔性桿的橫向變形量最大為7.73 mm；而剪切角均小于0.000 4°，符合細(xì)長梁的假設(shè)。從以上分析可知，基于逆動(dòng)力學(xué)解的控制方法是通過控制驅(qū)動(dòng)桿的關(guān)節(jié)變量來合理地控制彈性桿的變形狀態(tài)，從而保證末端執(zhí)行器準(zhǔn)確平穩(wěn)地按照理想圓周軌跡進(jìn)行運(yùn)動(dòng)。

圖11 驅(qū)動(dòng)桿與柔性桿的運(yùn)動(dòng)狀態(tài)Fig.11 States of actuated links and flexible links

約束方程和動(dòng)力學(xué)方程的誤差函數(shù)值如圖12所示，二者均低于預(yù)設(shè)值10-8和10-6，這也充分驗(yàn)證了逆動(dòng)力學(xué)瞬態(tài)剛體校正法的有效性。

圖12 約束方程與動(dòng)力學(xué)方程的迭代誤差Fig.12 Iterative errors of constraint and dynamic equations

4.3 運(yùn)動(dòng)控制實(shí)驗(yàn)

為驗(yàn)證求解策略的有效性，所搭建的實(shí)驗(yàn)平臺(tái)包含剛?cè)狁詈?-RRRU并聯(lián)機(jī)器人，其柔性桿的彈性模量為6.9×1010Pa，其他參數(shù)與計(jì)算實(shí)例中相同；驅(qū)動(dòng)關(guān)節(jié)通過高精密減速器由交流伺服電機(jī)驅(qū)動(dòng)，并采用絕對(duì)位置編碼器提供關(guān)節(jié)角度，分辨率為0.005 5°；多軸運(yùn)動(dòng)控制器采用Power PMAC，插補(bǔ)周期221 μs；測量設(shè)備采用FARO型激光跟蹤儀，重復(fù)定位精度為0.02 mm。實(shí)驗(yàn)前首先要對(duì)系統(tǒng)中的機(jī)械設(shè)備、電控設(shè)備和測量設(shè)備進(jìn)行調(diào)試和標(biāo)定。然后，將逆動(dòng)力學(xué)模型的穩(wěn)定因果解輸入給控制器進(jìn)行軌跡跟蹤控制。利用激光跟蹤儀對(duì)測量靶標(biāo)的圓心點(diǎn)進(jìn)行定位，實(shí)時(shí)記錄末端執(zhí)行器的空間運(yùn)動(dòng)軌跡，如圖13所示。由于定位靶標(biāo)的接觸面為球面，其圓心到末端執(zhí)行器的垂直距離為25 mm，因此，實(shí)際測量的空間軌跡與末端執(zhí)行器的運(yùn)動(dòng)軌跡會(huì)存在偏差。

圖13 激光跟蹤儀測量系統(tǒng)Fig.13 Laser tracker measuring system

實(shí)驗(yàn)包含兩部分：令動(dòng)平臺(tái)承載1.6 kg負(fù)載，首先，根據(jù)逆動(dòng)力學(xué)求解策略研究求解實(shí)例中圓周軌跡的跟蹤性能；其次，測量柔性桿上特征點(diǎn)的軸向主應(yīng)變，并與理論值作比較。為便于對(duì)比分析，實(shí)驗(yàn)過程中將分別基于理想剛性模型和剛?cè)狁詈夏Ｐ偷哪鎰?dòng)力學(xué)解來控制機(jī)器人，通過測量末端執(zhí)行器的空間位置、軌跡精度以及柔性桿上特征點(diǎn)的應(yīng)變值，來驗(yàn)證求解策略的有效性。

在研究軌跡跟蹤性能時(shí)，利用激光跟蹤儀分別測量兩種模型下的圓周軌跡，如圖14所示。圖中圓弧上引出的直線段表示運(yùn)動(dòng)軌跡與擬合圓之間的誤差，其長度代表擬合誤差的幅值。由表2的測量結(jié)果可知，二者在擬合直徑上的差距為0.153 mm，但后者的圓度誤差降低了1.46 mm，軌跡最大誤差降低了0.372 mm，最小誤差降低了0.029 mm。從以上對(duì)比可知，基于剛?cè)峄旌蟿?dòng)力學(xué)模型的控制效果優(yōu)于基于理想剛性模型的控制效果。這進(jìn)一步驗(yàn)證了所述控制方法的有效性。

圖14 帶載情況下的軌跡跟蹤實(shí)驗(yàn)Fig.14 Trajectory tracking experiments with payload

表2 圓周軌跡跟蹤測量結(jié)果Tab.2 Circular trajectory tracking results mm

由表2可知，基于剛?cè)狁詈夏Ｐ偷目刂品椒ㄔ谲壽E跟蹤精度上優(yōu)于基于剛性模型時(shí)的精度。這說明利用穩(wěn)定因果解實(shí)施控制能提高剛?cè)狁詈喜⒙?lián)機(jī)器人的空間軌跡精度，從而驗(yàn)證了該方法的有效性。必須指出，實(shí)驗(yàn)過程中定位靶標(biāo)的實(shí)際測量軌跡與末端執(zhí)行器的真實(shí)運(yùn)動(dòng)軌跡存在偏差。由于單臺(tái)激光跟蹤儀無法測量動(dòng)平臺(tái)的旋轉(zhuǎn)運(yùn)動(dòng)，所以該偏差很難通過定量計(jì)算加以分析和補(bǔ)償。為解決該問題，未來可通過修改動(dòng)平臺(tái)的結(jié)構(gòu)來提高系統(tǒng)的軌跡測量精度。

為進(jìn)一步驗(yàn)證動(dòng)力學(xué)模型及求解方法的有效性，在第1支鏈柔性桿外壁上選擇3個(gè)特征點(diǎn)，它們與上關(guān)節(jié)的距離依次為200、407.5、600 mm，在各個(gè)點(diǎn)處粘貼應(yīng)變片并利用DH5927N型動(dòng)態(tài)應(yīng)變儀對(duì)其軸向的主應(yīng)變進(jìn)行測量，如圖15所示。通過計(jì)算與測量發(fā)現(xiàn)，各特征點(diǎn)的主應(yīng)變理論值以及實(shí)際測量值的變化趨勢基本相同，因此以特征點(diǎn)1為例，對(duì)其主應(yīng)變的理論計(jì)算值和測量結(jié)果進(jìn)行分析和說明，如圖16所示。

圖15 柔性桿上的特征點(diǎn)Fig.15 Typical points on flexible link1.特征點(diǎn)1 2.特征點(diǎn)2 3.特征點(diǎn)3

通過對(duì)比圖16a和16b中特征點(diǎn)1的主應(yīng)變理論計(jì)算值與測量值可知，二者的變化趨勢基本相同，而且都是在同一個(gè)數(shù)量級(jí)上發(fā)生變化，這進(jìn)一步驗(yàn)證了理論模型與求解方法的有效性。但二者的變化范圍具有一定的差距，其原因之一是測量應(yīng)變前會(huì)進(jìn)行調(diào)零平衡操作，它會(huì)消除測點(diǎn)的初始應(yīng)變；另一個(gè)原因是建模過程中沒有考慮各運(yùn)動(dòng)關(guān)節(jié)的裝配間隙，運(yùn)動(dòng)過程中該間隙的變化會(huì)影響特征點(diǎn)的軸向主應(yīng)變，使其在幅值和方向上產(chǎn)生偏差。因此，可對(duì)該間隙進(jìn)行建模和補(bǔ)償來進(jìn)一步提高系統(tǒng)理論模型與軌跡跟蹤的精度。表3是各個(gè)特征點(diǎn)處測量的正向主應(yīng)變峰值和負(fù)向主應(yīng)變峰值，通過變化幅度的對(duì)比可知，柔性桿中點(diǎn)處特征點(diǎn)的主應(yīng)變變化幅度更小。

圖16 特征點(diǎn)1處的主應(yīng)變Fig.16 Axial principal stresses of typical point 1

帶載實(shí)驗(yàn)下的測量點(diǎn)正向峰值負(fù)向峰值特征點(diǎn)19.655×10-6-1.2739×10-5特征點(diǎn)24.692×10-6-5.4250×10-6特征點(diǎn)37.612×10-6-1.0965×10-5

5 結(jié)論

(1)針對(duì)空間剛?cè)狁詈喜⒙?lián)機(jī)構(gòu)，提出了基于瞬態(tài)剛體校正法的動(dòng)力學(xué)模型求解方法，該方法有效地解決空間閉鏈機(jī)構(gòu)多體動(dòng)力學(xué)模型在求解過程中出現(xiàn)的相容性問題，改善了系統(tǒng)綜合收斂性能。

(2)從系統(tǒng)能量、初值點(diǎn)的選取、求解的穩(wěn)定性與相容性等角度總結(jié)出剛?cè)狁詈蟿?dòng)力系統(tǒng)的求解策略，適用于求解具有閉鏈結(jié)構(gòu)的空間多體系統(tǒng)的穩(wěn)定因果解，由此構(gòu)建的控制策略可提高并聯(lián)機(jī)器人的軌跡跟蹤精度。

(3)針對(duì)3-RRRU剛?cè)狁詈喜⒙?lián)機(jī)器人進(jìn)行了仿真分析與實(shí)驗(yàn)對(duì)比驗(yàn)證，與基于剛性模型的控制方法相比，該方法將空間軌跡最大誤差降低了0.372 mm，圓度誤差降低了1.46 mm。通過對(duì)比柔性桿上特征點(diǎn)的主應(yīng)變理論值與測量值，進(jìn)一步驗(yàn)證了該理論模型與求解策略的有效性。

1 LONG P, KHALIL W, MARTINET P. Dynamic modeling of parallel robots with flexible platforms[J]. Mechanism and Machine Theory, 2014, 81：21-35.

2 鄭恩來, 張航, 朱躍, 等. 含間隙超精密壓力機(jī)柔性多連桿機(jī)構(gòu)動(dòng)力學(xué)建模與仿真[J/OL]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2017, 48(1)：375-385. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20170150 & journal_id= jcsam. DOI：10.6041/j.issn.1000-1298.2017.01.050.

ZHENG Enlai, ZHANG Hang, ZHU Yue, et al. Dynamic modeling and simulation of flexible multi-link mechanism including joints with clearance for ultra-precision press[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(1)：375-385. (in Chinese)

3 趙磊, 范夢然, 趙新華, 等. 柔性并聯(lián)機(jī)器人非線性摩擦動(dòng)力學(xué)建模與速度規(guī)劃[J/OL]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2017, 48(5)：390-396. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20170550&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2017.05.050.

ZHAO Lei, FAN Mengran, ZHAO Xinhua, et al. Nonlinear friction dynamic modeling and velocity planning of flexible parallel robot[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(5)：390-396. (in Chinese)

4 邱麗芳, 王棟, 印思琪, 等. Deform-X柔性鉸鏈設(shè)計(jì)與分析[J/OL]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2017, 48(4)：370-376. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20170449&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2017.04.049.

QIU Lifang, WANG Dong, YIN Siqi, et al. Design and analysis of Deform-X flexure hinge[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(4)：370-376. (in Chinese)

5 BAUCHAU O A, BETSCH P, CARDONA A, et al. Validation of flexible multibody dynamics beam formulations using benchmark problems[J]. Multibody System Dynamics, 2016, 37(1)：29-48.

6 MARTINI A, TRONCOSSI M, CARRICATO M, et al. Elastodynamic behavior of balanced closed-loop mechanisms: numerical analysis of a four-bar linkage[J]. Meccanica, 2014, 49(3)：601-614.

7 YU Y Q, DU Z C, YANG J X, et al. An experimental study on the dynamics of a 3-RRR flexible parallel robot[J]. IEEE Transactions on Robotics, 2011, 27(5)：992-997.

8 WU G L, BAI S P, HJ?RNET P. Architecture optimization of a parallel Sch?nflies-motion robot for pick-and-place applications in a predefined workspace[J]. Mechanism and Machine Theory, 2016, 106：148-165.

9 LIANG D, SONG Y M, SUN T. Nonlinear dynamic modeling and performance analysis of a redundantly actuated parallel manipulator with multiple actuation modes based on FMD theory[J]. Nonlinear Dynamics, 2017, 89(1)：391-428.

10 ORZECHOWSKI G, MATIKAINEN M K, MIKKOLA A M. Inertia forces and shape integrals in the floating frame of reference formulation[J]. Nonlinear Dynamics, 2017, 88(3)：1953-1968.

11 IRSCHIK H, KROMMER M, NADER M, et al. The equations of Lagrange for a continuous deformable body with rigid body degrees of freedom, written in a momentum based formulation[J]. Journal of Sound and Vibration, 2015, 335：269-285.

12 SHABANA A A. Dynamics of multibody systems[M]. New York：Cambridge University Press, 2013.

13 WANG Z, TIAN Q, HU H Y. Dynamics of spatial rigid-flexible multibody systems with uncertain interval parameters[J]. Nonlinear Dynamics, 2016, 84(2)：527-548.

14 CHEN T, WEN H, HU H Y, et al. Quasi-time-optimal controller design for a rigid-flexible multibody system via absolute coordinate-based formulation[J]. Nonlinear Dynamics, 2017, 88(1)：623-633.

16 WANG J L. Application of Radau IIA algorithms to flexible multibody system with holonomic constraints[J]. Nonlinear Dynamics, 2017, 88(4)：2391-2401.

17 HUSSEIN B, NEGRUT D, SHABANA A A. Implicit and explicit integration in the solution of the absolute nodal coordinate differential/algebraic equations[J]. Nonlinear Dynamics, 2008, 54(4)：283-296.

18 ERLICHER S, BONAVENTURA L, BURSI O S. The analysis of the Generalized-α method for non-linear dynamic problems[J]. Computational Mechanics, 2002, 28(2)：83-104.

19 ARNOLD M, BRüLS O. Convergence of the generalized-α scheme for constrained mechanical systems[J]. Multibody System Dynamics, 2007, 18(2)：185-202.

21 EBEL H, MATIKAINEN M K, HURSKAINEN V V, et al. Higher-order beam elements based on the absolute nodal coordinate formulation for three-dimensional elasticity[J]. Nonlinear Dynamics, 2017, 88(2)：1075-1091.

22 NACHBAGAUER K, GRUBER P, GERSTMAYR J. Structural and continuum mechanics approaches for a 3D shear deformable ANCF beam finite element：application to static and linearized dynamic examples[J]. Journal of Computational and Nonlinear Dynamics, 2012, 8(2)：021004-021004-7.