一种具有部分运动解耦和符号式位置正解的空间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法的优点。