基于绝对节点坐标法的空间并联机构刚柔耦合动力学建模

陈修龙 崔梦强

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

0 引言

并联机构具有工作精度高、末端执行器灵敏性好、负载能力强等优点，故而广泛应用于精密制造业和空间机械手等领域[1-5]。随着机器人技术的不断进步，对其设备创新性提出了更高的要求。具备结构冗余的空间并联机构可以较好地提高其机构刚度[6]，但在其高速运动时，并联机构中支链的弹性变形会产生系统振动、内应力增大和运动精度降低等不良影响[7-9]。

近年来，众多学者采用不同建模方法对并联机构的动力学性能进行了研究，柔性杆件的弹性变形问题对机构动力学的影响已成为机构学的研究热点[10-14]。王宗平等[15]提出具有六自由度的6PSS并联机构，基于拉格朗日方程并利用差商求导方式建立了机构的刚体动力学模型，并分析了运动学和动力学中各参数的变化。陈修龙等[16]基于牛顿欧拉法建立了4-UPS-UPU空间机构的动力学模型，分析了相关支链的受力情况，最后通过虚拟样机验证了Matlab理论结果的正确性。ZHANG等[17]基于有限元描述了柔性杆件梁单元模型，并通过拉格朗日法建立了考虑柔性杆件的动力学方程，采用一次迭代的方式求解其方程，并研究了弹性振动对机构的影响。WANG等[18]根据并联机构刚柔耦合的特点，基于有限元法通过离散柔性构件得到梁单元模型，然后利用KED法对3PRR平面并联机构的动力学方程进行组合装配，研究了其末端构件的运动特性。YU等[19]利用ANSYS Workbench软件对一种八连杆机构进行了刚柔耦合动力学分析，研究了柔性构件对机构传动特性曲线的影响，对轴孔尺寸进行优化发现，主连杆处的最大等效应力明显降低。LIU等[20]以空间三连杆焊接机器人为研究对象，基于拉格朗日法建立了刚柔耦合动力学模型，分析了柔性关节处不同刚度系数对焊接机器人运动精度的影响，结果表明，系统累积误差与运动轨迹有关。YAN等[21]提出一种用于臂架大范围回转运动和弹性振动的混合坐标模型，运用初应力法建立了钢丝绳的刚度矩阵，对超长柔性臂的动力学研究表明，随着载荷的增大，其偏航角的周期和幅值也相应增大。CHEN等[22]结合ANSYS和ADAMS对2PRR平面并联机构进行刚柔耦合动力学仿真，结果验证了模型运动学和动力学正逆解的正确性，最后还研究了柔性杆件对从、动滑块运动特性的影响。姜振海等[23]采用绝对节点坐标方法(ANCF)推导了平面3RRR并联机构考虑柔性杆件的动力学方程，利用直接积分方法(BDF)对方程进行求解，发现柔性构件的弹性变形对动平台影响很大，同时还发现主动杆件比从动杆变形严重。XING等[24]研究了6PSS并联机构刚柔耦合的动力学特性，结合ADAMS和ANSYS对机构进行运动仿真求解，将得到的响应曲线和机构动力学模型结果进行相互验证，证明了模型的正确性。

上述弹性动力学的研究对象大多数为串联机构和平面并联机构，而对于具有结构冗余的空间并联机构的建模分析研究较少。本文以具有结构冗余3-RRPaR空间并联机构为研究对象，基于绝对节点坐标法建立柔性梁单元模型，利用拉格朗日乘子方法推导机构的刚柔耦合动力学方程，并通过广义α方法求解其动力学方程，以分析构件材料的弹性模量对机构动力学的影响。

1 基于ANCF方法的梁单元模型

1.1 ANCF梁单元空间模型

采用三维二节点的梁单元模型对柔性杆件进行建模，其中梁单元模型如图1所示，图中OXYZ为系统坐标系；梁单元的两端节点分别为点i和点j；ri和rj分别为单元两节点在OXYZ下的位置矢量。

该模型通过插值法可得出单元上随机一点的位置坐标(位移场)为

(1)

(2)

其中

(α=i,j;β=x,y,z)S=[s1s2…s8]⊗I3

(3)

s1=1-3ζ2+2ζ3s2=l(ζ-2ζ2+ζ3)s3=l(ψ-ζψ)s4=l(ξ-ζξ)s5=3ζ2-2ζ3s6=l(-ζ2+ζ3)s7=lζψs8=lζξζ=x/lψ=y/lξ=z/l

式中mi、ni、ci——插值系数，i=0,1,…,7

S——梁单元形函数

e——梁单元节点坐标

⊗——克罗内克积运算符

I3——三阶单位矩阵

l——梁单元变形前长度

从式(2)可看出，梁单元的广义坐标个数共有24个，分别由节点处的位置坐标和描述3个方向的坐标矢量组成。

1.2 单元质量矩阵

梁单元质量矩阵可通过其动能求出，在计算梁单元动能前需求出其单元任一点的速度，计算式为

(4)

其动能计算式为

(5)

(6)

式中Ve——单元体积

ρ——杆件材料密度

Me——单元质量矩阵

1.3 单元广义弹性力

如图2所示，当点P0通过一段距离u后到点P，根据连续介质力学[25]相关知识可得到变形梯度J表达式为

(7)

其中

r0=Se0

式中r0——变形体运动前点P0的位置矢量

e0——梁单元在初始状态下的位置坐标向量

J0——单位矩阵

v——单元坐标系下的位置矢量

图3所示单元坐标系与系统坐标系重合时，可得到用右柯西-格林表达的变形张量为

C=JTJ

(8)

可推导得拉格朗日应变张量为

(9)

其中

式中Sα——形函数矩阵的第α行元素，α=1,2,3

为了便于后续积分计算，将式(9)应变张量转换成向量表达形式为

εk=[ε11ε22ε332ε122ε132ε23]T

(10)

通过式(10)推导出梁单元的弹性变形能为

(11)

式中H——构件材料的弹性常数矩阵

对单元变形能求导可得到单元弹性力为

(12)

其中

式中K(e)——单元弹性力刚度矩阵，非常数矩阵，具有时变性

1.4 单元广义外力

外力F作用在单元上任意一点时，可通过虚功原理求得该点的力F对单元所做的虚功为

(13)

通过式(13)可得系统坐标系OXYZ下沿X轴正方向的重力的一般表达形式为

(14)

式中me——梁单元质量

2 刚柔耦合动力学模型

2.1 自由度分析

如图4所示，该机构由定平台、3条型式相同的支链和动平台组成，其中每条支链包括主动臂和具有平行四边形构型的从动臂，可以看出该机构具有结构冗余，有利于增大机构运作过程中构件的刚度，提高稳定性。选择修正的Grubler-Kutzbach公式[26]计算自由度P，计算式为

(15)

式中n——机构构件总数，取17

m——机构运动副总数，取21

fi——第i个运动副的自由度

s——冗余约束个数，取12

将各个参数数值代入式(15)中，可求得本机构自由度为3，且均是移动方向。

2.2 约束方程

多体系统的约束方程需要根据运动副种类的不同建立相关的位置及方向约束，每个运动副相邻的两构件需要在其自身建立局部坐标系以便描述其在空间的位置及姿态。由于本机构运动副均为转动副，故刚性构件用自然坐标法(NCF)描述，柔性构件采用ANCF描述，其机构简图如图5所示，以此确定每个构件在系统坐标系下的位姿。

该并联机构由主动臂作为驱动构件，为了便于约束方程的建立及后续的计算仿真，主动臂采用质心坐标法描述，其中一条由定平台上一点到动平台上一点组成的支链的坐标命名，如图6所示，可知机构广义坐标为

q=(qN,qe)

(16)

式中qN=(Oi1,Oi2,Oi3,Oi4,Oi5,O′i) (i=1,2,3)，主动臂坐标包含6个广义坐标，其他均为3个广义坐标；qe=(Oi6,Oi7,Oi8,Oi9)(i=1,2,3)，每个节点的广义坐标个数均为12，故该机构一共有207个广义坐标，但为了方便描述广义坐标与自由度的关系，取节点的前3个位置坐标作为约束广义坐标，可以看出该并联机构一共有99个约束广义坐标。

本机构共有3个移动自由度，因此建立线性无关的约束方程96个，以及3个驱动约束。Oi(ai,bi,ci)(i=1,2,3)是定平台上的3个固定点，主动臂处的转动副约束为

(17)

(18)

(19)

其中

式中Ri——第i个支链主动臂变换矩阵

li1——第i个支链主动臂杆长

O——零矩阵

对于刚性构件与柔性构件的约束，可采用对应节点重合的形式进行约束，即

(20)

Φ5i=[Oi3-Oi2]T(Ai[0,0,1]T)=O1×1

(21)

(22)

(23)

式中Ai——柔性杆件的变换矩阵[27]

根据动平台的几何特性利用NCF法建立三角板单元，其约束方程为

(24)

(25)

该并联机构通过运动学反解可求出每条支链大臂的运动函数，可建立其驱动方程的约束为

(26)

式中θi——关于时间的驱动函数，i=1,2,3

由式(17)～(26)可组成3-RRPaR并联机构刚柔耦合的约束方程为

(27)

2.3 质量矩阵

质量矩阵是反映构件在空间中质量分布的最佳方法，是构建机构动力学方程不可或缺的一部分。质量矩阵为常数矩阵时可提高仿真计算效率，故该并联机构采用NCF和ANCF方法对构件进行描述，并且柔性构件的常数质量矩阵没有离心力和科氏力。

对于主动臂采用常见的质心坐标方法进行描述，故其质量矩阵为

Mk=diag(mk,mk,mk,Ikx,Iky,Ikz) (k=1,6,11)

(28)

对于NCF方法描述的刚性杆件，采用两点零矢形式通过惯性力的虚功率可求得刚性杆单元质量矩阵为

(29)

刚性三角板单元质量矩阵为

(30)

对于ANCF方法描述的柔性梁单元，由式(6)可求得其对应的质量矩阵为

(31)

由式(28)～(31)可组建机构的质量矩阵为

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

(32)

2.4 动力学方程

选择多体系统动力学中指标-3的DAEs描述动力学方程，基于拉格朗日乘子法构成的3-RRPaR空间并联机构刚柔耦合动力学方程为

(33)

式中Φq——约束方程Φ(q,t)关于广义坐标的导数

λ——拉格朗日乘子

QΣ——系统的广义外力

FΣ——式(12)组装的广义弹性力

由于式(33)DAEs不能直接进行求解，需要经过广义α方法进行差分直接离散成代数方程的形式求解，具体方法见文献[28]，将式(33)构建成矩阵的形式为

(34)

对式(34)进行求导可得其对应的Jacobian矩阵为

(35)

3 数值仿真与分析

3.1 仿真参数

利用三维软件建立3-RRPaR空间并联机构模型，经测量可得构件材料相关参数如表1～4所示。

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

表2 构件结构参数Tab.2 Structural parameters of components

表3 定、动平台坐标参数Tab.3 Coordinate parameters of fixed and moving platforms m

表4 柔性杆件材料参数Tab.4 Material parameters of flexible rod

3.2 仿真结果与分析

并联机构拥有3个平移方向的自由度，为了便于分析，规划动平台为一个圆形轨迹，其轨迹函数表达式为

(36)

利用Matlab对其进行多体系统刚柔耦合动力学建模，并通过广义α方法进行数值求解，其中求解步长为0.002 s，谱半径是0.7，为了避免泊松闭锁其泊松比取值0[28]。

不同弹性模量E下动平台的位移响应曲线如图8所示。由8a、8c可知，3种不同弹性模量与刚体机构的位移曲线基本重合。从图8b、8d可以发现，弹性模量为116 GPa时与理想情况下的位移误差曲线最为接近，由图8d可得，最大位移误差仅为-2.2×10-7m，在弹性模量为10 GPa时，最大误差为-6.615×10-6m，弹性模量为1 GPa时最大误差为-9.579×10-6m。综上分析发现弹性模量与理想时的位移误差基本在数量级10-5m以内，其弹性模量对动平台位移的影响很小，体现出该具有结构冗余的并联机构有着较高的精度。

图9为3种弹性模量下的速度响应曲线。由图9a、9c可以看出，速度曲线一开始有微小波动，但后面的曲线基本重合，其波动原因可能是柔性杆件突然变形导致的数值解不收敛，并且Z轴方向与Y轴方向相比波动较大。从图9b、9d发现，随着弹性模量的减小速度误差变大；由图9d可知，在弹性模量为1 GPa时幅值的跨距最大，且最大误差达到-6.098×10-4m/s，弹性模量为10 GPa时最大误差为3.92×10-4m/s，弹性模量为116 GPa时最大误差为2.158×10-4m/s，但是3种误差基本在数量级10-3m/s以内，可得弹性模量对动平台速度的影响大于位移的影响。

图10为3种弹性模量下的加速度响应曲线。由图10a、10c可得，其曲线一开始有着较大波动，Y轴和Z轴两方向的最大波动分别达到了-34.78 m/s2和-248.3 m/s2，但是在经历0.1 s后曲线基本重合；由图10b、10d可知，随着弹性模量的减小其加速度误差变大；由图10d可得，弹性模量为1 GPa时误差最大可达到0.109 1 m/s2，弹性模量为10 GPa时最大误差为0.035 67 m/s2，弹性模量为116 GPa时最大误差为0.014 72 m/s2。综合对动平台的位移、速度、加速度曲线分析可发现，考虑柔性构件的运动特性曲线与刚体机构的动力学特性曲线基本一致，通过误差发现随着弹性模量的减小其位移、速度和加速度曲线偏离理想状态的程度越大，且弹性模量对加速度的影响最大，速度次之，位移最小。

采用第四强度理论的等效应力进行分析，其中在不同弹性模量的影响下该机构各支链处杆件的等效应力如图11所示。第1支链开始时，E=116 GPa的等效应力达到8.926 MPa，但之后随着弹性模量减少其等效应力曲线波动范围增大，并且不同曲线之间的差距随着时间的增大逐渐变小。由图11b可知，第2支链在E=1 GPa时等效应力最大达到了2.446 MPa，并且弹性模量越小，对应杆件的等效应力曲线波动范围也就越大。由图11c可知，第3支链在E=116 GPa时等效应力最大达到了15.01 MPa，其3种弹性模量对应的曲线规律与图11a类似，发现由于第1支链与第3支链的对称性，其两支链的等效应力曲线具有明显的相似性。综上分析发现其柔性杆件的弹性模量越小，对应等效应力曲线的波动范围也就越大，同时在曲线稳定时支链处3种弹性模量对应的等效应力基本在数量级1 MPa以内。

4 结论

(1)基于绝对节点坐标法描述了三维二节点柔性梁单元空间模型，并根据功能关系推导了梁单元的质量矩阵，根据连续介质力学推导了弹性力矩阵，最后利用虚功原理推导得出重力方向上单元广义外力的表达式。

(2)分析了3-RRPaR空间并联机构的结构特征，采用自然坐标法和绝对节点坐标法建立了机构约束方程和具有常数的质量矩阵，最后基于拉格朗日乘子法推导得出机构的刚柔耦合动力学方程，并在Matlab中对方程进行数值求解。

(3)得到不同弹性模量下的动力学响应曲线，3种模量下的运动曲线与理想轨迹基本重合，验证了所建模型的正确性，结果表明，柔性构件的弹性模量越小，对机构运动精度的影响越大，且弹性模量对机构动平台加速度的影响最大，对位移的影响最小，体现出具有结构冗余的空间并联机构精度高、稳定性好的特点。给出了机构在高速运动过程中等效应力的变化情况，并对比分析了3种不同弹性模量对等效应力的影响。