3-RRPaR并联机构刚体动力学建模与分析

陈修龙 郭景尧 贾永皓

(山东科技大学机械电子工程学院， 青岛 266590)

0 引言

并联机构相比于串联机构具有刚度大、精度高、承载能力强等优点，但由于并联机构奇异点的影响，使其优势难以发挥[1-4]。具有冗余结构的并联机构可以避免并联机构的奇异位型，提高并联机构的刚度，在一定程度上提高了并联机构的运动性能，因此受到研究者的广泛关注[5-7]。结构冗余并联机构通过构型的优势弥补了普通并联机构的不足，因此，对结构冗余并联机构的研究具有重要意义。

并联机构的刚体动力学建模是实现动力学控制和研究机构性能的前提，是在机构运动学的基础上对并联机构的进一步研究与探索[8-11]，可为并联机构驱动电机的选择、控制和实际机构的研制奠定理论基础，同时为其他结构冗余并联机构的刚体动力学建模提供可行的思路。近年来，国内外学者采用不同的方法对并联机构动力学模型进行了大量研究，但针对具有冗余结构的并联机构研究较少。刘芳华等[12]基于拉格朗日动力学模型，建立了6-UPS并联平台的全动力学模型，并求解了理想驱动力和时间的关系。XIN等[13]以三自由度空间并联机器人为例，利用拉格朗日与虚功原理相结合的方法建立了动力学模型，并通过仿真和实验加以证明。桑董辉等[14]提出了一种新型的二自由度球面并联机构，利用牛顿-欧拉法建立了动力学方程，并在给定外力的作用下计算了关节作用力和作用力矩。王英波等[15]利用凯恩法建立了Stewart平台的刚体动力学模型，并利用动力学建模仿真验证了其动力学模型的正确性。陈修龙等[16]采用牛顿-欧拉法和虚拟仿真相结合的方法对4-UPS-RPS空间并联机构进行了动力学分析。ABDELLATIF等[17]利用拉格朗日公式和广义坐标重新推导了动力学方程，验证了坐标选取形式的重要性。朱伟等[18]以一种新型四自由度高速并联机构为例，运用虚功原理建立了动力学方程，并通过ADAMS虚拟仿真进行了验证。GOSSELIN等[19]介绍了一种具有冗余结构的并联机构，并证明该机构所有奇异点都可通过运动学冗余来避免。ABADI等[20]证明了结构冗余是避免奇异点和增大工作空间的有效途径，并提出了一种新的运动规划策略，以避免奇异点的产生。QU等[21]提出了一种平面结构冗余并联机构，并通过静力学仿真预测了机构中较为容易损坏的位置。

本文以具有冗余结构的3-RRPaR空间并联机构为研究对象，考虑其冗余构型的影响，利用拉格朗日乘子法建立该并联机构的刚体动力学模型。首先，分析3-RRPaR空间并联机构的结构特征，计算机构整体自由度；然后，利用闭环矢量法建立并联机构的运动学模型，求解出驱动杆的角位移、角速度和角加速度；采用质心坐标和自然坐标相结合的方式建立系统约束方程，计算机构整体的质量矩阵，利用拉格朗日乘子法建立并联机构的动力学方程，并根据拉格朗日乘子与约束力矩的关系求解驱动力矩；最后，利用Matlab对动力学方程进行数值求解，将理论模型求解结果与ADAMS虚拟仿真结果进行对比，分析并联机构动力学响应曲线。

1 结构特征

3-RRPaR空间并联机构如图1所示。该并联机构由定平台、动平台和3条运动支链组成，每条运动支链包括主动臂和从动臂，从动臂为由4根杆件组成的平行四边形构型，其中R表示转动副，Pa表示平行四边形构型，该并联机构的各个构件之间均通过转动副连接。

图1 3-RRPaR空间并联机构的三维模型Fig.1 Three dimensional model of 3-RRPaR spatial parallel mechanism1.主动臂 2.转动副 3.从动臂 4.定平台 5.动平台

3-RRPaR为闭环空间并联机构，采用修正的Grubler-Kutzbach公式计算其自由度[22]，计算公式为

(1)

式中P——机构整体自由度

p——机构构件总数

q——机构运动副总数

fk——第k个运动副的自由度

h——冗余约束数

并联机构的3条运动支链分布在2个平面内，增加了机构的刚度和平稳性，通过对3个主动臂的旋转驱动可实现动平台的三维移动。对于运动支链中的平行四边形构型，由于引入了冗余约束，在动力学建模过程中需去掉多余的约束方程。

2 运动学模型

机构运动学是研究物体的位置、速度和加速度与时间的关系，是机构动力学研究的基础。本文利用闭环矢量法建立并联机构运动学模型[23-25]，分别求出机构位置、速度和加速度的逆解。

2.1 位置模型

为方便表示并联机构构件间的几何关系，对该机构作机构简图如图2所示。在定平台下建立全局坐标系Oxyz，动平台视为质点O′，则O′点的位置矢量g=(x,y,z)T可表示为

gi=Mi+l1iui+l2isi-Ni(i=1,2,3)

(2)

式中Mi——铰点Mi在Oxyz中的位置矢量

Ni——铰点Ni在O′x′y′z′中的位置矢量

l1i——机构主动臂杆长

l2i——机构从动臂杆长

ui——主动臂单位方向矢量

si——从动臂单位方向矢量

图2 3-RRPaR空间并联机构简图Fig.2 3-RRPaR spatial parallel mechanism diagram

Mi在定系Oxyz下的位置矢量为

Ni在局部坐标系O′x′y′z′下的位置矢量为

为方便计算，引入单位矢量wi，wi垂直于主动臂的旋转面，使得wi×ui能够表示主动臂的旋转方向。ui和wi的表达式为

将式(2)化简，得gi-Mi-l1iui+Ni=l2isi，等式两边左乘各自的转置，可得

(3)

将式(3)转换为三角函数式

Aisinθi+Bicosθi+Ci=0

其中A1=-2l11x

B1=2l11(y+a-m)

A2=2l12(z-b+n)

B2=2l12y

A3=-2l13x

B3=-2l13(y-a+m)

解得主动臂的角位移为

(4)

因此，可得从动臂单位方向矢量si为

si=(g-Mi-l1iui+Ni)/l2i

(5)

2.2 速度模型

将式(2)两端关于时间求导，得到动平台质心点O′的速度表达式为

(6)

ωi——从动臂角速度矢量

(7)

对式(6)两端同时左叉乘ωi，可求得机构从动臂角速度为

(8)

则机构从动臂质心速度为

(9)

2.3 加速度模型

将式(6)两端关于时间求导，得到动平台质心点O′的加速度表达式为

(10)

(11)

(12)

则机构从动臂质心加速度为

(13)

通过对并联机构进行运动学逆解，得到主动臂角位移、角速度和角加速度。

3 动力学模型

3.1 约束方程

根据空间并联机构各构件之间的几何关系建立其约束方程。质心坐标和欧拉角转换矩阵可以准确描述构件的位置和姿态，易于通过驱动臂的角度变化来实现机构驱动，但对于特殊情况下的旋转变换，欧拉角存在奇异问题。自然坐标是利用全局坐标系下的绝对坐标来描述构件的运动状态，两点零矢的自然坐标不涉及到构件的方向矢量，适用于描述空间构件的相对一维转动，但两点零矢的自然坐标不能通过角度变化来实现机构驱动。

3-RRPaR空间并联机构为结构冗余的并联机构，各个构件之间均通过转动副连接。为了易于驱动，方便地构建约束方程并去掉多余的约束，本文对3个支链的主动臂采用质心坐标和欧拉角来表示机构构件的位置和姿态，对其他构件采用自然坐标来表示机构构件的位置和姿态。建立自然坐标时，对于杆构件，在杆的两端建立自然坐标，对于三角板构件，在三角板的3个端点处建立自然坐标。则机构整体的广义坐标表示为

q=(O1,O2, …,O30)

其中Ok=(xk,yk,zk)T(k=2,3,…,9,11,12,…, 18,20,21,…,30)，Ok=(xk,yk,zk,αk,βk,γk)T(k=1,10,19)。

该并联机构共有99个广义坐标，且机构有三维移动自由度，因此需要添加96个线性无关的约束方程。由于机构复杂且广义坐标数目较多，本文以1号支链为例构建其约束方程，1号支链构件坐标表示如图3所示。

图3 并联机构1号支链坐标示意图Fig.3 Coordinate diagram of No.1 branch chain of parallel mechanism

用质心坐标和欧拉角表示构件位姿时，可以利用运动副的约束关系建立约束方程。空间转动副只有一维转动，限制了构件的三维移动和二维转动，因此需构建5个约束方程。为了约束两构件之间的相对移动，需使两构件在旋转轴中心上的两点重合。为了约束两构件之间两个方向的相对转动，需沿旋转轴在两构件上分别建立两个单位方向向量，使得这两个方向向量共线，由于已知1号支链主动臂绕z轴旋转，因此可直接限制构件另外2个方向的相对旋转。可得转动副约束方程为

(14)

其中

式中Ru1——1号支链主动臂绕z轴的旋转矩阵

O3×1、O2×1——零矩阵

l11——1号支链主动臂杆长

利用自然坐标表示构件位姿时，可直接根据构件的几何关系建立约束方程

(15)

(16)

Φ4=O4-O2=O3×1

(17)

Φ5=O6-O3=O3×1

(18)

Φ6=|O5-O4|-l21=0

(19)

Φ7=O8-O5=O3×1

(20)

Φ8=O9-O7=O3×1

(21)

(22)

(23)

式中l21——1号支链从动臂杆长

lx——平行四边形构型中小杆杆长

根据动平台的尺寸和运动状态建立约束方程

(24)

(25)

由运动学模型可知主动臂角位移变化函数，则对机构动力学模型建立驱动约束方程为

(26)

式中Φd——机构主动臂驱动约束方程

由式(26)组成3-RRPaR空间并联机构理想约束方程

Φ(q)=[Φ1Φ2…Φ32Φd]T=O99×1

(27)

式中Φ1～Φ10——1号支链约束方程

Φ11～Φ20——2号支链约束方程

Φ21～Φ30——3号支链约束方程

3.2 质量矩阵

质量矩阵是质量到广义坐标的推广，是机构动力学方程的重要组成部分，而质心坐标和自然坐标可以方便地表示构件的质量矩阵。

对于质心坐标下构件对应的质量矩阵，由构件局部坐标原点的移动分量对应的质量矩阵和构件局部坐标原点的转动分量对应的惯性张量矩阵组成，可分别表示为

Nj=diag(mj,mj,mj) (j=1,6,11)

(28)

Jj=diag(Ixj,Iyj,Izj) (j=1,6,11)

(29)

则质心坐标下构件j对应的质量矩阵为

Mj=diag(Nj,Jj) (j=1,6,11)

(30)

对于自然坐标下构件对应的质量矩阵，可利用惯性力的虚功率分别求出杆单元和三角板单元的质量矩阵，杆单元和三角板单元(动平台)的质量矩阵分别为

(31)

(32)

式中I3——三阶单位矩阵

则系统整体质量矩阵为

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

(33)

3.3 动力学方程

利用拉格朗日乘子法来建立并联机构刚体动力学方程。将式(27)对时间求一阶导数，得到系统速度约束方程为

(34)

式中Φq——系统约束方程的Jacobian矩阵

将式(27)对时间求二阶导数，得到系统加速度约束方程为

(35)

利用拉格朗日法建立3-RRPaR空间并联机构带有拉格朗日乘子的动力学方程为

(36)

式中δ、η——大于0的修正参数[25]

λ——系统拉格朗日乘子矢量

Q——系统广义外力矢量

求解该动力学方程可得到系统的广义坐标对应的广义加速度和拉格朗日乘子矢量。通过系统的广义加速度，积分得到系统的广义速度和广义坐标，可用来描述间隙对系统动力学的影响，为研究含间隙动力学打下基础。通过拉格朗日乘子矢量，可求得主动臂转动副处的驱动力矩，完成3-RRPaR空间并联机构动力学逆解。则i号支链驱动杆转动副处驱动力矩为

(37)

Ri——i号支链主动臂局部坐标系关于全局坐标系的转换矩阵

ρi——i号支链主动臂转动铰点在其局部坐标系中的位置矢量

4 算例仿真与分析

4.1 仿真参数

根据上述理论模型利用Matlab编写程序，进行3-RRPaR空间并联机构动力学方程的数值求解。利用Matlab自带的ode45算法进行数值求解，使得计算效率更高，结果更精确。

3-RRPaR空间并联机构的构件转动惯量和构件参数如表1～3所示。

表1 3-RRPaR空间并联机构主动臂转动惯量Tab.1 Moment of inertia of active arm of 3-RRPaR spatial parallel mechanism kg·m2

图4 动平台动态响应曲线Fig.4 Dynamic response curves of moving platform

表2 3-RRPaR空间并联机构的结构参数Tab.2 Structural parameters of 3-RRPaR spatial parallel mechanism

表3 定平台和动平台尺寸参数Tab.3 Size parameters of fixed platform and moving platform m

4.2 仿真结果与分析

为了分析并联机构的动力学响应，需要对动平台规划一条运动轨迹。已知并联机构有三维移动自由度，使得动平台中心点在全局坐标系Y-Z平面运动轨迹为圆轨迹，其运动方程描述为

(38)

动平台在全局坐标系下的起始点为(0.54, 0.15, 0.08)，利用Matlab对动力学方程进行数值求解，得到该机构动平台位移、速度和加速度曲线，然后利用拉格朗日乘子和约束力矩的关系求解出驱动杆处的驱动力矩，并将结果与ADAMS仿真结果进行对比，对机构动态响应曲线进行分析。

图4是该机构动平台动态响应曲线，分别求出了机构刚体时动平台的位移、速度和加速度，与规划的运动轨迹相同，验证了该动力学模型的正确性，为分析含间隙时的动力学模型打好基础。

图5 并联机构驱动杆的驱动力矩曲线Fig.5 Driving torque curves of driving rod

图5a为该并联机构空载时驱动杆的驱动力矩曲线。从图5a可以看出，Matlab求解动力学模型所得驱动力矩曲线和ADAMS虚拟仿真所得驱动力矩曲线基本吻合，但有较小误差。其中3号驱动杆的驱动力矩误差最大，为0.027 N·m；1号驱动杆的驱动力矩误差次之，为0.015 3 N·m；2号驱动杆的驱动力矩误差最小，为0.002 3 N·m。空载时，3号支链驱动杆的驱动力矩最大，最大值达到4 N·m，1号支链和2号支链驱动力矩较小且相差不大，为0.8 N·m左右。

图5b为该并联机构动平台加载1 kg时驱动杆的驱动力矩曲线。从图5b可以看出，动平台加载1 kg之后，驱动杆的驱动力矩数值相应增大，Matlab与ADAMS数值结果的误差也随之增大，但幅值和周期基本没有变化。3号驱动杆的驱动力矩最大为6 N·m ，误差为0.027 N·m ；1号驱动杆和2号驱动杆驱动力矩相差不大，为1 N·m，误差也相应较小，1号驱动杆误差为 0.014 N·m，2号驱动杆误差为0.004 N·m。

由以上分析可知，通过刚体动力学模型求解出的并联机构动平台的位移、速度和加速度，结果与所规划轨迹完全吻合，分别对比了机构空载和加载时驱动杆的驱动力矩，Matlab与ADAMS结果基本一致，验证了3-RRPaR空间并联机构刚体动力学模型的正确性。

5 结论

(1)分析了3-RRPaR并联机构的结构特征，利用闭环矢量法建立了该机构的逆运动学模型，得到了主动臂的角位移、角速度和角加速度。

(2)通过质心坐标法和自然坐标法建立了并联机构的约束方程，分别计算出质心坐标和自然坐标下构件的质量矩阵，利用拉格朗日乘子法建立了3-RRPaR空间并联机构的动力学方程，并根据拉格朗日乘子和驱动力矩的关系求解出驱动杆的驱动力矩。

(3)通过Matlab求解3-RRPaR并联机构的刚体动力学方程，得到并联机构动平台的位移、速度和加速度，结果与所规划轨迹完全吻合，通过ADAMS虚拟仿真得到机构空载和加载时驱动杆的驱动力矩曲线，与Matlab数值计算结果进行对比分析，验证了3-RRPaR空间并联机构刚体动力学模型的正确性。