劉 涼 趙新華 周海波 王嘉斌

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

0 引言

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

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

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

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

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

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

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

(1)

其中

(2)

(3)

式中Sr——剛性梁的形函數(shù)，Sr∈R3×6

qi——剛性梁的自然坐標(biāo)矢量，qi∈R6

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

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

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

(4)

式中ρi、mi——剛性梁密度、質(zhì)量

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

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

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

r=Cqa

(5)

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

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

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

(6)

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

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

(7)

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

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

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

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

Mij=?VijρijSijTSijdVij

(8)

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

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

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

(9)

式中m、n——梁單元的首端點、末端點

rijm、rijn——首末端點處絕對位置矢量

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

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

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

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

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

2.2 并聯(lián)機器人正動力學(xué)模型

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

(10)

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

(11)

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

Qg——廣義重力矢量

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

Qe——廣義彈性力矢量

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

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

λ——拉格朗日乘子矢量

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

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

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

(12)

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

(13)

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

2.3 并聯(lián)機器人逆動力學(xué)模型

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

(14)

式中ΦT——新添加的3個約束方程

(15)

式中xp、yp、zp——動平臺中心點P處運動軌跡

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

3 剛?cè)狁詈蟿恿W(xué)系統(tǒng)的求解策略

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

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

(16)

(17)

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

(18)

(19)

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

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

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

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

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

(20)

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

(21)

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

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

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

(22)

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

(23)

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

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

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

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

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

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

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

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

4 求解實例與實驗分析

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

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

4.1 正動力學(xué)模型求解實例

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

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

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

圖7 動平臺法向量與ZG軸的夾角Fig.7 Angle between platform normal vector and ZG-axis

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

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

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

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

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

4.2 逆動力學(xué)模型求解實例

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

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

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

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

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

4.3 運動控制實驗

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

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

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

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

圖14 帶載情況下的軌跡跟蹤實驗Fig.14 Trajectory tracking experiments with payload

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

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

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

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

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

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

帶載實驗下的測量點正向峰值負向峰值特征點19.655×10-6-1.2739×10-5特征點24.692×10-6-5.4250×10-6特征點37.612×10-6-1.0965×10-5

5 結(jié)論

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

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

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

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 鄭恩來, 張航, 朱躍, 等. 含間隙超精密壓力機柔性多連桿機構(gòu)動力學(xué)建模與仿真[J/OL]. 農(nóng)業(yè)機械學(xué)報, 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)機器人非線性摩擦動力學(xué)建模與速度規(guī)劃[J/OL]. 農(nóng)業(yè)機械學(xué)報, 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/OL]. 農(nóng)業(yè)機械學(xué)報, 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.