机械臂系统的动力学建模与分析

苏华勇，李海艳

（广东工业大学机电工程学院，广州 510006）

0 引言

机械臂具有控制范围广及其操作灵活性高等优点，在工业各领域中得到广泛应用，并且机械臂系统在当下的发展需求材料更轻质、操作范围更广及运行精度更高[1]。因此，需考虑机械臂在运动过程中各节臂发生弹性变形对位置控制的影响，这就是典型的柔性多体动力学问题，其主要特点是机械臂系统在运动的过程中不但经历着大的刚性转动，也伴随着变形运动，并且这两种运动是高度耦合的[2]。机械臂的动力学建模是机械臂运动分析的基础，典型的动力学分析方法主要有拉格朗日法、牛顿-欧拉法和Kane 方法等等，目前较常用的是拉格朗日法，它在分析过程中不必求机械臂的内作用力而只需求其对应的速度[3]。由于弹性变形有无限多自由度，动力学模型的精确解是无法得到的，通常将弹性变形离散成有限自由度来近似[4-6]。典型的离散方法有模态分析法[7]和瑞利-里兹法[8]等等，目前常用的是通过瑞利-里兹法离散弹性变形。为了方便模型的求解，常取弹性变形离散形式前两阶近似，但没有分析不同阶数对模型求解准确度的影响。此外，动力学模型属于高度复杂的非线性方程组，从而忽略了模型中一些小变形项以简化分析的难度，但这同样影响了动力学模型求解的准确度。

本文以三节臂机械臂为研究对象，如图1 所示。通过拉格朗日法和瑞利-里兹法推导出准确的动力学模型，分析弹性变形对机械臂运动的影响及弹性变形离散不同阶数对模型求解准确度的影响。

图1 三节臂机械臂

1 机械臂动力学建模

在拉格朗日法中，只要求出机械臂系统的总动能Ek和总势能Ep，得到它们的差值L＝Ek-Ep，最后建立拉格朗日动力学方程：

式中：Qi为系统的广义坐标；Q˙i为对应的广义速度；Fi为系统的广义力；n为完整约束方程的个数。

1.1 参考坐标系

如图2 所示，分析机械臂系统首先需要建立一系列的坐标系来表示各机械臂的位置，其中xoy 为整体坐标系，x1o1y1、x2o2y2、x3o3y3为局部坐标系，局部坐标系分别以各臂杆产生弹性变形的首末端的连线作为坐标系的横轴，这样建立坐标系的好处是各臂杆自身产生的弹性变形不会影响到其他臂杆的位置描述，从而降低建立动力模型的复杂度。

图2 机械臂坐标系示意图

在机械臂系统中，把各臂杆均考虑为均质杆，驱动电机及机械手夹持装置均考虑为集中质量，并设臂杆i 的质量为mi，长度为li，电机i的质量为mi-1，机械手夹持装置的质量为M3，臂杆i 在整体坐标系中的偏转角为θi，将臂杆的位置向量从局部坐标系过渡到整体坐标系需要借助旋转矩阵Ai：

其中i＝1，2，3。

1.2 机械臂的总动能

假设机械臂系统各臂杆在相对坐标系上一点的位置向量为pi[ xi, 0 ]T，则产生弹性变形wi( xi, t )后的位置为pi′＝[xi,wi(xi,t)]T，末端位置表示为 pli＝[li,0]T，借助旋转矩阵可将在局部坐标系上的位置转换成在整体坐标系上的位置，转换坐标系后的位置向量分别表示如下：

通过对位置向量关于时间t 求导后便可得到对应的速度。因此，臂杆的总动能表示为：

式（6）的前一项表示臂杆i 在运动时产生的动能，后一项表示表示相应臂杆末端负载的动能。

1.3 机械臂的总势能

机械臂在运动过程中，总势能的变化除了各臂杆运动时产生的重力势能Eg1，还有臂杆产生弹性变形产生的重力势能Eg2和变形产生的应变能Es，各势能可以表示为：

式中：EIi为臂杆i的刚度；E为弹性模量；I为惯性矩。

因此，机械臂总势能可以表示为：

1.4 瑞利-里兹法

由于机械臂臂杆的弹性变形是与距离和时间有关的量，通常需要将弹性变形离散成有限自由度来近似，典型的离散方法有模态分析法和瑞利-里兹法。目前常用的是瑞利-里兹法，它也是弹性连续体力学中求近似解的最基本方法。该方法是将弹性变形项离散成里兹基函数与对应时间系数乘积的级数展开式，以往的方法一般都是取级数展开式的前2阶进行近似。为了进行对比验证，此次选取了级数展开式前4阶来近似弹性变形项，各臂杆离散后的弹性变形可表示如下：

1.5 建立动力学方程

将弹性变形的离散形式分别代入到动能和势能的求解式中，最后得到的总动能和总势能表达式都是仅与时间t有关的多项式。

拉格朗日函数可表示为L＝Ek-Ep，方程的广义坐标表示为Q，机械臂系统的广义坐标除了臂杆的偏转角外，还有考虑各臂杆发生弹性变形造成的变化，它也是受时间t影响的量，在弹性变形的离散形式中以各阶基函数对应的时间系数表示。方程的广义力表示为F，在旋转运动中的广义力为驱动转矩，并设电机的转矩为T，臂杆1 受到电机1 的驱动转矩和电机2反作用的驱动转矩，臂杆2受到电机2的驱动转矩和电机3的反作用驱动转矩，臂杆3 受到电机3 的驱动转矩。因此臂杆i受到的驱动转矩τi为：

最后，机械臂的拉格朗日动力学方程表示如下：

式中：Q＝[θ,q]T， θ＝[θ1,θ2,θ3]， q ＝ [q11,q12,q13,q14,q21,，。

从式（13）可以看到，机械臂的动力学模型包含了15个方程，并且该模型属于复杂非线性强耦合的方程组。它既包含了机械臂各臂杆间刚性运动的耦合，还包含了刚性运动与臂杆弹性变形的耦合。前者表示机械臂的大范围刚性运动，而后者会使机械臂系统运动产生扰动。该模型考虑了弹性变形的影响，因此对机械臂系统描述更加全面，但是也很大程度地增加模型的求解难度。

图3 动力学模型求解流程图

1.6 动力学方程求解

以往的方法对方程组的求解过程中，为了方便求解就一定程度地简化了模型，把方程中的小变形项都省略了，但是这无疑对方程的准确度有一定影响。因此，动力学方程的求解不应该舍弃小变形项对模型的贡献。

通过Mathematica对公式推导和整理，拉格朗日方程组最后可以整合成My¨＝z 的形式，其中M15×15是个满秩的对称矩阵，是广义坐标的二阶导数项，将其他剩余的项整合到z15×1里面，且方程组存在唯一解。

由于该方程组属于非线性的常微分方程组，可以借助MATLAB 软件的ode 函数求解，但是方程一般需要满足格式y˙＝f (t,y) 的形式，因此要对动力学方程作进一步处理并求解，流程如图3所示。

2 数值仿真

现在通过具体的参数进行数值仿真，三节臂机械臂的长度分别为 l1＝l2＝l3＝1.5 m，各臂杆的质量分别为 m1＝m2＝3 kg，m3＝2.5 kg，电机的质量分别为M0＝M1＝M2＝5.5 kg，末端夹持装置的质量 M3＝4.5 kg，重力加速度 g＝9.8 m/s2，刚度 EI1＝EI2＝EI3＝406 N·m2，电机驱动力矩T1＝1 150 sin(2πt)3+850(N·m)；T2＝452sin(2πt)3+450(N·m)；T3＝155 sin(2πt)3+150(N·m)。仿真时间为1.5 s，并设初始条件Q ＝Q˙＝ 0 。

为了分析产生变形前后对位置的影响，考虑将机械臂的末端位置进行对比，可以表示为：

如图4 所示，机械臂末端位置在变形前后的曲线变化存在一定的差异，在变形前的曲线是很光滑的，机械臂弹性变形后的曲线在局部放大处更是可以看到很明显的振动，这正是弹性变形产生的抖动现象导致的。另外，机械臂各臂杆的偏转角在变形前变化后角度也产生了一定的偏差。在变形前的曲线也是很光滑的，在机械臂发生弹性变形后在局部放大处也是可以看到很明显的锯齿状。这些现象都可以说明机械臂系统发生的弹性变形对位置控制有一定程度影响。

图4 机械臂变形前后对比

由于弹性变形有无限多自由度，动力学模型的精确解是无法得到的，通常将弹性变形离散成有限自由度来近似。为了分析不同阶数近似的弹性变形对求解的影响，用弹性变形近似的前2阶和前4阶进行对比，其误差分析如图5所示。

可以看出，选取前2阶和前4阶近似弹性变形对模型求解准确度的影响并不大。但是，选取展开式近似的阶数越大，动力学模型规模就会更大更复杂，求解的难度也会陡增，因此选取阶数多少需要根据实际情况进行综合权衡。

图5 弹性变形2阶与4阶近似误差

3 结束语

本文通过以三节臂的机械臂作为研究对象，考虑弹性变形对机械臂系统的影响，建立动力学模型并进行数值仿真求解与分析，可以得到如下的结论。

（1）机械臂系统在运动过程中产生的弹性变形对臂杆的位置控制产生的扰动较大。因此，机械臂系统的控制有必要考虑弹性变形的影响，弹性变形的扰动特性需要进一步分析。

（2）机械臂的弹性变形离散近似取前2阶和前4阶对模型的求解准确度影响不大，并且阶数增加会导致模型更加复杂，阶数的选取也需要根据实际情况综合权衡。