一種具有部分運動解耦和符號式位置正解的空間2T1R并聯機構拓撲設計與動力學建模

黃凱偉 沈惠平 李 菊 朱忠頎 楊廷力

常州大學現代機構學研究中心，常州，213016

0 引言

少自由度并聯機構具有驅動件少、控制簡單、制造容易等優勢，已成為機構學領域的研究熱點之一，其中，具有三自由度的兩平移一轉動(2T1R)并聯機構可用于空間抓放定位、調姿等場合，應用范圍廣泛[1]，該機構主要分為平面型、空間型兩種：平面型2T1R機構主要有典型的3-RRR并聯機構及其衍生機構[2-4](如3-PRR、3-RPR等)，其轉動軸線方向平行于動平臺的法線；空間型2T1R并聯機構[5-10]的轉動軸線平行于動平臺的平面。

并聯機構的動力學建模對機構的動態性能分析和實時控制策略至關重要[11]，常用的動力學分析方法有Newton-Euler法[12-13]、動力學普遍方程[14-15]、Lagrange方程[16-17]、Hamilton正則方程[18]、Kane方法[19-20]等。

SHIAU等[12]采用New-Euler法對3-RPS混聯機構進行了動力學建模；郝秀清等[13]運用Newton-Euler法推導出了3-PTT并聯機構的動力學方程；KALANI等[14]基于虛功原理提出了一種能縮短計算時間、提高精度的改進型動力學普遍方程，并對Gough-Stewart機構進行了正逆動力學分析；賈曉輝等[15]依據虛功原理構建了3-RRPR柔性機構的動力學逆解模型，并利用其中的質量矩陣和剛度矩陣確定了系統固有頻率的求解表達式；THANH等[16]通過Lagrange方程解決了一種冗余并聯機構的動力學建模問題；徐奕柳等[17]基于Lagrange方程建立了PURU+RR+S球面踝關節機構學模型；尤晶晶等[18]利用Hamilton正則方程對并聯式六維加速度傳感器進行了動力學研究；姚建濤等[19]采用Kane方法分別建立了基座固定和基座運動兩種情況下并聯調整機構的動力學方程；ZHAO等[20]建立了基于Kane方程的四足機器人動力學模型，得到了每個滑塊上的驅動力。

上述方法中，有些建模步驟較為繁瑣，需要對所有構件進行力學分析(如Newton-Euler法)；有些無法求出系統運動副的支反力(如Lagrange方程)；而有些則需要引入稍復雜的偏速度、廣義速率等概念(如Kane方法)。

楊廷力[21]提出了一種基于虛功原理的序單開鏈法，該方法以子運動鏈(sub-kinematic chain，SKC)為基本單元，建模思路清晰，能得到若干重要運動副(即連接各SKC的運動副)處的支反力，這對機械結構的強度設計是至關重要的。但該方法自推出以來主要用于平面機構，還未曾用于空間機構的逆向動力學建模。

本文首先根據基于方位特征(position and orientation characteristic，POC)方程的并聯機構拓撲設計理論及方法[22]，設計并提出了一種含有一條冗余驅動支鏈、零耦合度且運動部分解耦的空間兩平移一轉動(2T1R)并聯機構。該機構的優勢是：①零耦合度使機構具有符號式位置正解；②含三個最小SKC，且分別含有驅動副，使機構具有部分運動解耦特性；③全由轉動副R構成，使機構易于制造；④依據冗余驅動消除奇異位置原理[23-24]，冗余驅動支鏈能避開機構的奇異位置，且能使機構整體的剛度特性提高20%左右[7]。該機構可設計成具有較大剛度的管材彎曲加工并聯裝備，可應用于需要制造小批量、多品種曲線形狀管材的汽車、五金、電子產品等領域。然后根據基于拓撲特征的運動學建模原理，給出了該機構的符號式位置正解，同時，基于雅可比矩陣求出了機構各構件的速度/加速度值，并利用基于虛功原理的序單開鏈法[21]建立了該空間機構的逆向動力學模型，求解分析了機構的驅動力矩及部分運動副的支反力。

1 2T1R并聯機構的設計及其拓撲分析

1.1 機構設計

根據基于POC方程的并聯機構拓撲設計方法以及受文獻[7]的啟發，本文設計了一種全由轉動副組成、零耦合度且部分運動解耦的2T1R并聯機構，如圖1所示。

圖1 全轉動副空間型2T1R并聯機構的設計Fig.1 Design of spatial 2T1R parallel mechanism with full rotation pairs

動平臺1由三條支鏈(Ⅰ、Ⅱ、Ⅲ)連接于靜平臺0，其中，混合支鏈Ⅰ為一平面五桿機構(R11—R12—R13—R22—R21)，也可視為由兩條相同的RRR支鏈(即支鏈A、B)串聯而成；簡單支鏈Ⅱ由三個相互平行的R副(R31—R32—R33)串聯而成；冗余支鏈Ⅲ也由三個相互平行的R副(R41—R42—R43)串聯而成。動平臺1上的轉動副R13、R23處為復合鉸鏈，該機構中各個轉動副的軸線均平行于靜平臺0所在平面的y軸。

1.2 拓撲結構分析

1.2.1機構的POC計算

并聯機構的POC方程為

Mbk=∪MJi

(1)

Mpa=∩Mbk

(2)

式中，MJi為第i個運動副的POC集；Mbk為第k條支鏈末端的POC集，k取Ⅰ、Ⅱ、Ⅲ；Mpa為機構動平臺的POC集。

首先，計算各條支鏈的POC集。第Ⅰ條混合支鏈由兩條相同的RRR支鏈(即支鏈A、B)串聯而成，由式(1)可得

其中，t2表示二維移動；r1表示一維轉動；t2(⊥R13)表示在垂直于R13軸線方向的平面內存在二維移動；r1(∥R13)表示平行于R13軸線方向存在一維轉動。

第Ⅱ、Ⅲ條支鏈分別由三個相互平行的R副串聯而成，其末端的POC集分別為

然后，計算由支鏈Ⅰ、Ⅱ及轉動副R23構成的子并聯機構的POC集，可表示為

(3)

其中，MJR23為轉動副R23的POC集。

最后，計算并聯機構動平臺的POC集。由式(2)可得

(4)

由式(3)和式(4)可知，支鏈Ⅰ、Ⅱ及轉動副R23組成子并聯機構后，即可實現兩平移一轉動的輸出運動，支鏈Ⅲ并不會影響機構動平臺的POC。

1.2.2機構的自由度分析

首先，確定第1個獨立回路的自由度F1。混合支鏈Ⅰ為第1個獨立回路(即第1個子并聯機構)，也是第1條單開鏈(single open chain，SOC)。混合支鏈Ⅰ是一個平面五桿機構，因此第1個獨立回路的獨立位移方程數ξL1=3，則第1個子并聯機構的自由度為

式中，fi為第i個運動副的自由度；c為運動副數；ξLj為第j個SOC的獨立位移方程數。

然后，確定第2個獨立回路的自由度F2。混合支鏈Ⅰ、轉動副R23及支鏈Ⅱ組成第2個獨立回路(即第2個子并聯機構)，其獨立位移方程數為

其中，dim{·}表示POC集的維數。則第2個子并聯機構的自由度為

最后，確定并聯機構整體的自由度Fpa。第2個子并聯機構與支鏈Ⅲ組成第3個獨立回路，其獨立位移方程數為

則機構整體的自由度為

因此，該機構整體的自由度為3，當取R11、R21、R31作為驅動副時，動平臺1可實現oxz平面內的兩平移和繞y軸旋轉(2T1R)的輸出運動。

由于支鏈Ⅲ既不影響機構整體的自由度，又不影響機構動平臺的POC，因此支鏈Ⅲ為冗余支鏈。此外，支鏈Ⅲ上可以存在驅動，只有當機構處于奇異位置時，該驅動電機才會工作，并帶動機構避開奇異位置，因此，支鏈Ⅲ可以作為冗余驅動支鏈。

1.2.3機構的耦合度計算

由基于單開鏈單元的機構組成原理[22]可知，任意機構可分解為若干個子運動鏈(SKC)，每一個SKC可分解為約束度為正、零、負的單開鏈(SOC)，將第j個SOC的約束度定義為

(5)

式中，cj為第j個SOC的運動副數；Ij為第j個SOC的驅動副數。

因此，機構的耦合度κ定義為

(6)

1.2.2節已求得ξL1=ξL2=ξL3=3，由式(5)得到各個回路的約束度分別為

上述三個回路分別構成對應的SKC1、SKC2和SKC3，因此，該機構包含3個SKC，它們的耦合度可由式(6)求得，即

因此，機構的耦合度為零。三個SKC連接處的運動副R23、R43為本機構受力的薄弱處，本文將R23、R43定義為重要運動副。重要運動副的受力求出后，各SKC內其他運動副的受力則易求出。

至此，該機構的三個主要拓撲特征(方位特征POC、自由度F、耦合度κ)已求出[25]，為后續的運動學建模和動力學建模奠定了基礎。

2 機構的位置分析

2.1 位置正解

2.1.1基于拓撲特征的運動學建模原理

由基于單開鏈的機構組成原理[22]可知，機構可分解為若干個SKC，而每個SKC又可分解出約束度為正值、零、負值3種形式的SOC，因此，機構位置正解的求解可轉換為3種不同形式SOC的位置求解。對于本機構，3個SKC的耦合度均為零(即所有SOC的約束度為零)，因此，可直接求出各個SKC的符號式位置正解，從而求出整個機構的符號式位置正解[26]。

2T1R機構的運動學建模如圖2所示(為便于后文表達，圖2中各點表示圖1中對應位置上運動副的中心點，如點A1為運動副R11的中心點，其他類同)，靜平臺0是一個以Ad(d=1，2，3，4)為頂點構成的正方形，其邊長為l1；動平臺1是一個以Cd為頂點構成的正三角形，其邊長為l3；其余各桿件的長度均為l2，設桿件AdBd與x軸負方向的夾角分別為θd。

圖2 2T1R機構的運動學建模Fig.2 Kinematics modeling of 2T1R mechanism

靜坐標系{o}(oxyz)的原點在靜平臺0的質心o處，且x軸與A1A2連線平行，y軸與A1A2連線垂直；動坐標系{p}(puvw)在動平臺1的質心p處建立，且u軸與C3C4連線平行，v軸指向點C1，z軸、w軸的方向由右手定則確定。

設動坐標系{p}的原點p相對于靜坐標系{o}的原點o的位置坐標為p=(xp，yp，zp)，動平臺繞y軸旋轉的角度(即姿態角)為α。

2.1.2SKC1(A1—B1—C1—B2—A2)的位置求解

由拓撲分析可知，SKC1由運動副R11、R12、R13、R22、R21組成，可表示為R11—R12—R13—R22—R21，其中各個運動副的中心點依次為A1、B1、C1、B2、A2，則SKC1可簡寫為A1—B1—C1—B2—A2，其他類同。

已知θ1、θ2及各桿件的尺寸參數，可得

A1=(l1/2，l1/2，0)TA2=(-l1/2，l1/2，0)T

A3=(-l1/2，-l1/2，0)TA4=(l1/2，-l1/2，0)T

B1=(l1/2-l2cosθ1，l1/2，l2sinθ1)T

B2=(-l1/2-l2cosθ2，l1/2，l2sinθ2)T

由1.2節可知，動平臺會產生xoz平面內的兩平移和繞y軸旋轉的2T1R輸出運動，則

由桿長約束B1C1=B2C2=l2，建立如下位置方程：

解得

(7)

m=1+a2n=2(ab-azB1-xB1)

由式(7)可以看出，動平臺1的兩維平移運動量xp、zp僅由驅動副R11、R21的輸入角θ1、θ2確定，因此，該機構具有輸入-輸出部分運動解耦性。

2.1.3SKC2(C2—C3—B3—A3)的位置求解

由桿長約束條件B3C3=l2，可得

(xp+l1/2+l2cosθ3)cosα=0

整理可得

Asinα+Bcosα+C=0

A=zp-l2sinθ3B=-xp-l1/2-l2cosθ3

S=zp-l2sinθ3R=xp+l1/2+l2cosθ3

解得

(8)

從而求得該機構動平臺的姿態角α。

2.1.4SKC3(A4—B4—C4)的位置求解

設冗余驅動支鏈Ⅲ的隨動角為θ4，則點B4、C4的坐標分別為

其中，oC4、pC4分別為點C4在靜坐標系{o}和動坐標系{p}下的坐標。

由桿長約束B4C4=l2，可得

Dsinθ4+Ecosθ4+F=0

D=2l2(l3sinα/2-zp)

E=-2l2(l1/2-xp-l3cosα/2)

F=(l1/2-xp-l3cosα/2)2+(l3sinα/2-zp)2

解得

(9)

至此，機構位置正解已全部完成。

2.2 位置逆解

已知動平臺1上p點坐標(xp，yp，zp)和姿態角α，求輸入角θ1、θ2、θ3。

由2.1節已求得B1、B2、B3、C1、C2、C3點坐標，求解由B1C1、B2C2、B3C3的桿長約束建立的位置方程，可得

(10)

t1=t2=2l2zpt3=2l2zp+l2l3sinα

M3=2l2(-xp+l3cosα/2-l1/2)-

(-xp+l3cosα/2-l1/2)2-(zp+l3sinα/2)2

N3=2l2(-xp+l3cosα/2-l1/2)+

(-xp+l3cosα/2-l1/2)2+(zp+l3sinα/2)2

由此可求得輸入角θ1、θ2、θ3，且由式(9)可得隨動角θ4的值。

綜上所述，當給定動平臺1上p點的坐標(xp，yp，zp)和姿態角α時，輸入角θ1、θ2、θ3以及隨動角θ4各有兩組解，故逆解數為24=16，即機構有16種構型。

2.3 正逆解驗算

取機構尺寸參數分別如下：l1=400 mm，l2=200 mm，l3=462 mm。取輸入角θ1、θ2、θ3分別為10°、20°、30°。由MATLAB軟件計算2.1節的位置正解公式可得機構在理論上存在4組位置正解，見表1。通過對比三維CAD模型可知第4組數據是實際位置正解。

表1 機構的正解數值

取表1中第4組數據代入位置逆解公式(式(10))，可求得輸入角θ1、θ2、θ3的8組理論逆解數值，見表2。

表2 機構的逆解數值

由表2可知，第4*組逆解數值與正解求解時設定的3個輸入角一致，因此，正逆解公式推導正確。

3 機構的速度與加速度分析

3.1 動平臺的速度與加速度

Jpν=Jqω

(11)

f11=-2(xB1-xp)f12=-2(zB1-zp)f13=0

f21=-2(xB2-xp)f22=-2(zB2-zp)f23=0

f31=-2(xB3-xC3)f32=-2(zB3-zC3)

f33=l3sinα(xC3-xB3)+l3cosα(zC3-zB3)

g11=2l2sinθ1(xB1-xp)+2l2cosθ1(zB1-zp)

g22=2l2sinθ2(xB2-xp)+2l2cosθ2(zB2-zp)

g33=2l2sinθ3(xB3-xC3)+2l2cosθ3(zB3-zC3)

當機構不存在奇異位置時Jp可逆，可得

(12)

依據式(12)即可求得動平臺原點的輸出速度。

則動平臺的移動速度矩陣、轉動速度矩陣與三個輸入角之間的關系可表示為

(13)

(14)

將式(12)對時間t求導，得到動平臺p點加速度與輸入加速度之間的映射關系為

(15)

其中，k1、k2、k3為矩陣Jp和Jq中各項元素對時間t的導數(k1、k2、k3的具體值可掃描本文首頁OSID二維碼獲得)。

3.2 桿件的速度與加速度

3.2.1桿件A1B1的速度與加速度

B1點的速度為

vB1=vA1+ω11×(l2c11)

(16)

式中，vA1為驅動副R11的線速度，因驅動副在機架上，故vA1=0；ω11為驅動桿A1B1的角速度；c11為桿件A1B1的單位向量。

對式(16)求導，可得B1點的加速度為

aB1=aA1+l2ε11×c11+l2ω11×(ω11×c11)

(17)

式中，aA1為驅動副R11的線加速度；ε11為驅動桿A1B1的角加速度。

于是，桿件A1B1質心的速度、加速度分別為

vmid，11=vB1/2

(18)

amid，11=aB1/2

(19)

3.2.2桿件B1C1的速度與加速度

vC1=vB1+ω12×(l2c12)

(20)

在式(20)等號兩邊同時叉乘c12，可得桿件B1C1的角速度為

(21)

在式(20)等號兩邊同時對時間t求導，可得

(22)

在式(22)等號兩邊同時叉乘c12，可得桿件B1C1角加速度為

(23)

于是，桿件B1C1質心的速度、加速度分別為

3.2.3其余構件的速度與加速度

其余構件的速度和加速度的求法與前文類似，故直接給出結果。桿件AdBd(d=2，3，4)質心的速度、加速度分別為

(24)

(25)

vBd=vAd+ωd1×(l2cd1)

aBd=aAd+l2εd1×cd1+l2ωd1×(ωd1×cd1)

其中，cd1、ωd1、εd1分別為相應構件的單位向量、角速度、角加速度。

桿件BdCd(d=2，3，4)質心的速度、加速度分別為

(26)

(27)

其中，vCd、aCd分別為點Cd的已知速度和加速度。則桿件BdCd的角速度、角加速度分別為

(28)

(29)

3.3 算例與仿真

給定3個驅動副的運動規律分別為：θ1=πt/18，θ2=πt/15，θ3=πt/10。

依據式(12)，通過MATLAB計算得到動平臺p點的理論速度值，并與ADAMS仿真得到的速度值進行對比，如圖3所示；依據式(15)，通過MATLAB計算得到動平臺p點的理論加速度值，并與ADAMS仿真得到的加速度值進行對比，如圖4所示。

圖3 動平臺的理論速度與仿真速度曲線Fig.3 The theoretical and simulation speed curve of the moving platform

圖4 動平臺的理論加速度與仿真加速度曲線Fig.4 The theoretical and simulation acceleration curve of the moving platform

由圖3和圖4可以發現：該機構速度、加速度的理論值與對應仿真值基本一致，由此驗證了運動學模型的正確性。

4 機構的動力學分析

將主動臂在驅動力矩作用下產生的轉角定義為廣義坐標，即q=(θ1,θ2,θ3)T，所對應的廣義虛位移為Δq=(Δθ1,Δθ2,Δθ3)T。依據式(16)～式(29)，各構件的虛位移與機構的廣義虛位移之間的關系可表示為：ΔX=J1Δq，Δα=J2Δq，ΔXde=JdevΔq，Δθde=JdeωΔq(d=1，2，3，4；e=1，2，其中，e=1表示桿件AB，e=2表示桿件BC)。其中，桿件A1B1的移動和轉動虛位移分別為ΔX11、Δθ11，移動和轉動速度雅可比矩陣分別為J11v、J11ω；桿件B1C1的移動和轉動虛位移分別為ΔX12、Δθ12，移動和轉動速度雅可比矩陣分別為J12v、J12ω；桿件A2B2的移動和轉動虛位移分別為ΔX21、Δθ21，移動和轉動速度雅可比矩陣分別為J21v、J21ω；其余桿件類同。本文采用基于虛功原理的序單開鏈法[21]來構建2T1R機構的動力學模型。

4.1 基于虛功原理的序單開鏈法基本原理

4.2 受力分析

作用于構件質心上的力有重力和慣性力，而作用于構件質心上的力矩則僅為慣性力矩。對于動平臺，作用在質心的力和力矩分別為

(30)

對于各支鏈，假設重力是唯一的外力，則作用在各構件上的力和力矩分別為

Fde=mdeg-mdeamid，de

(32)

Mde=-oIdeεde-ωde×(oIdeωde)

(33)

式中，mde為對應構件的質量；amid，de為對應構件質心的加速度；oIde為靜坐標系{o}下各構件質心處的慣量矩陣；ωde、εde分別對應構件的角速度和角加速度；d=1,2,3,4;e=1,2。

4.3 動力學方程

4.3.1SKC3的動力學求解

對于SKC3，解除運動副R43處的約束，將點C4處的運動副支反力FC4轉化為未知外力，由虛功原理可得

(34)

其中，ΔXC4為點C4處的移動虛位移。

在靜坐標系{o}下，以桿B4C4為研究對象，并將桿B4C4上所有的力對點B4取矩，可得

oI42ε42+ω42×(oI42ω42)

(35)

其中，c42為桿件B4C4的單位向量。

聯立式(34)、式(35)，可計算出點C4處的運動副支反力為

q1=xC4-xB4q2=yC4-yB4q3=zC4-zB4

其中，E1、E2、E3分別為單位矩陣E3×3的第1、2、3行向量；vC4x、vC4z分別為點C4在x、z方向的速度分量。

4.3.2SKC2的動力學求解

對于SKC2，解除運動副R23處的約束，將點C2處的運動副支反力FC2轉化為未知外力，由虛功原理可得

(36)

其中，ΔXC2為點C2處的移動虛位移，M3為驅動副R31的驅動力矩。

在靜坐標系{o}下，以動平臺為研究對象，并將動平臺上所有的力對點C3取矩，可得

(37)

其中，c43、c23、cp3分別為點C4、C2、p到點C3的單位向量。

4.3.3SKC1的動力學求解

SKC1上無需解除運動副，可直接求解，有

(38)

4.4 數值仿真算例

本文機構各構件的尺寸參數如表3所示。本文僅考慮動平臺的自身重力(即f=0、τ=0)，采用與3.3節相同的運動規律，首先通過MATLAB編程計算得到驅動力矩的理論值，然后將虛擬樣機導入ADAMS中，設定機構運動副的約束類型，施加豎直向下的重力，并選取仿真步長為0.01 s、仿真時間為5 s，對虛擬樣機進行動力學仿真研究，得到驅動力矩的仿真值，如圖5所示。同理可分別得到運動副R43、R23處支反力的理論值和仿真值(支反力的曲線圖可掃描本文首頁OSID二維碼獲得)。

將MATLAB理論計算值與ADAMS仿真值進行對比，結果表明理論值與仿真值基本吻合，從而驗證了動力學模型的正確性。存在差異的原因如下：關節運動副存在的間隙使得理論計算的尺寸參數與模型仿真時存在少許差異；ADAMS軟件仿真采用的是Lagrange方程，與本文采用的基于虛功原理的序單開鏈法在計算上存在少許的舍入和累計誤差。

4.5 與傳統虛功原理建模方法的比較

依據傳統虛功原理，可得

化簡可得

(39)

式(39)對任意Δq都成立，因此，可得機構的逆向動力學方程為

(40)

本文研究結果表明，通過式(40)得到的結果與圖5中驅動力矩理論值一致。傳統的虛功原理建模方法采用的是整體建模思路，不區分建模的先后順序，方程中也沒有直接體現求解支反力這一要素，而本文采用的是基于虛功原理的序單開鏈法，該方法按照機構拓撲結構分解的順序，有序地進行動力學分析，建模思路清晰，使得機構的結構學、運動學及動力學具有統一性，并且可以直接求解出連接不同SKC處的運動副支反力，求出這些支反力后，各SKC內其他運動副的受力也便于求出，這有助于機械結構設計。

5 結論

(1)根據基于方位特征(POC)方程的并聯機構拓撲設計理論及方法，設計并提出了一種轉動軸線方向垂直于動平臺法線的空間布置兩平移一轉動(2T1R)并聯機構，該機構具有符號式位置正解、部分運動解耦特性、被動冗余支鏈能避免奇異位置、剛度性能好等特點。該機構可設計成具有較大剛度的管材彎曲加工并聯裝備，可應用于需要制造小批量、多品種曲線形狀管材的汽車、五金、電子產品等領域。

(2)根據基于拓撲特征的運動學建模原理，給出了該機構的符號式位置正解，采用基于虛功原理的序單開鏈法，建立了該機構的逆向動力學模型，并求解出了機構的驅動力矩以及重要運動副處的受力。

(3)與傳統虛功原理建模方法相比，基于虛功原理的序單開鏈法的優勢之處在于：①建模思路清晰，體現了機構結構學、運動學及動力學的統一；②既能求解驅動力矩，又能求解出重要運動副處的支反力，即兼有Newton-Euler法和Lagrange法的優點。