一种基于机械臂柔性变形的工业机器人动力学建模方法，包括以下步骤：首先根据ABAQUS软件建立模型，确定模型网格参数、动力学建模参数和边界条件，然后使用控制方程求解机械臂柔性变形后的单元节点位置，进而根据正运动学计算出旋转后的单元节点位置，对此时的单元逆向运动求解单元节点纯变形，计算出单元节点内力；根据惯性力确定单元节点外力，最后对单元节点位置和计算时间进行更新；对上述步骤进行迭代计算，直至计算时间达到最终时间。本发明将机械臂杆件柔性考虑到动力学模型中，有效模拟出因机械臂柔性产生的末端抖动现象，为工业机器人振动机理的研究及提高定位精度提供了理论基础；并且推导流程简洁，适合编程计算，具有实用意义。
1.一种基于机械臂柔性变形的工业机器人动力学建模方法，特征在于该方法包括以下步骤：步骤1，采用ABAQUS软件使用六面体单元网络建立工业机器人物理模型，模型所在全局T T T坐标系O‑XYZ的坐标轴方向为：ex＝[1 0 0] ,ey＝[0 1 0] ,ez＝[0 0 1] ，获得模型网格参数，确定动力学建模参数和边界条件；
步骤2，根据步骤1的模型网格参数，动力学建模参数以及边界条件，使用结合牛顿运动学和显式中心差分的控制方程求解机械臂柔性变形后的单元节点位置；
式中， 为第n+1步机械臂柔性变形后的单元节点i的位置； 为第n步机械臂柔性变n形后的单元节点i的位置；n为计算步数；ξ为阻尼系数,h为计算步长,m为单元节点质量；Fin为单元节点i的单元节点力， fi 为单元节点i的单元节点内力， 为单元节点i的单元节点外力，i＝a,b,c...,g,h；
步骤3，根据正运动学对单元节点位置 和 以关节为旋转中心旋转角度θ，获得旋转后的机械臂单元节点位置 和步骤4，以第n+1步的单元质心为参考点，将第n+1步的单元分别对应逆向运动至第n步的单元位置上，使得相邻两步的单元充分重合，进而根据经过逆向运动后虚拟单元节点位置与第n步的单元节点位置的差值计算出单元节点i的纯变形qi；
步骤5，以第n步的单元质心为原点建立局部坐标系 坐标轴方向与全局坐标系O‑XYZ相同，将第n+1步单元节点位置和单元节点纯变形转换到局部坐标系中，采用八节点六面体单元形函数建立形函数坐标系O”‑rst，并将局部坐标系内单元节点位置通过雅克比矩阵转换到形函数坐标系中，根据形函数导数计算出单元应变增量，进而计算出单元变形后的应力，利用单元变形描述的虚功和应力应变描述的虚功相等这一原理，计算出虚拟单元节点位置下的单元节点内力，进一步对其进行正向运动，计算出实际位置下的单元节点内力，根据惯性项计算出单元节点外力；
步骤6，在第n步到第n+1步，经过旋转后的机械臂单元节点位置 和 分别作为下一计算步的单元节点位置 和更新计算时间time＝nh；
步骤7，判断计算时间time是否达到最终时间ET；若否，转入到步骤2到步骤6继续迭代计算；若是，计算结束，输出计算结果。
2.如权利要求1所述的一种基于机械臂柔性变形的工业机器人动力学建模方法，其特征在于，所述模型网格参数包括模型单元序号、单元节点序号和单元节点坐标值；所述动力学建模参数包括运动时间ET、材料密度rho、弹性模量E和泊松比pr；所述边界条件为对模型单元节点运动方向添加约束。
3.如权利要求1所述的一种基于机械臂柔性变形的工业机器人动力学建模方法，其特征在于，步骤3中正运动学方法的步骤为：建立关节转动轴线向量RA，确定回转中心点xcet；
式中，θ为旋转角度， RA＝[RAx RAy RAz] ；
4.如权利要求1所述的一种基于机械臂柔性变形的工业机器人动力学建模方法，其特征在于，步骤4中逆向运动步骤为：逆向平移，包括第一次逆向转动和第二次逆向转动；
逆向平移是将第n+1步的单元，以单元质心 为参考点，整体逆向平移至第n步的单元质心 上，经过逆向平移的单元节点位置x′i为：第一次逆向转动首先对经过逆向平移的单元，以单元平面质心为参考点构造参考平面及其法向量vec′1，采用相同的方式对第n步的单元构造参考平面及其法向量vec1,进而计算出两法向量的夹角θA及其垂直向量nvecA:进而计算第一次逆向转动矩阵R′(‑θA):2
R′(‑θA)＝[1‑cos(‑θA)]V′+sin(‑θA)V′式中，
经过第一次逆向转动的单元节点位置x″i为：x″i＝R′(‑θA)x′i
第二次逆向转动除选取的参考平面与第一次逆向转动不同外，其余步骤均相同，经过第二次逆向转动的单元节点位置x″′i为：x″i＝R″(‑θB)x″i
式中，R″(‑θB)为第二次逆向转动矩阵，θB为参考平面法向量夹角；
5.如权利要求1所述的一种基于机械臂柔性变形的工业机器人动力学建模方法，其特征在于，步骤5中将第n+1步单元节点位置转换到局部坐标系公式为：式中， 为局部坐标系下的单元节点i的位置，八节点六面体单元形函数Ni表达式为：
式中，r、s、t分别为形函数表达式参数，ri＝±1，si＝±1，ti＝±1为形函数坐标系内节点i的坐标；
形函数导数计算步骤为：对雅克比矩阵J求逆，建立形函数关于局部坐标系的导数与形函数关于形函数坐标的导数等式为：单元变形后的应力 计算公式为：
式中， 为第n步单元的初始应力值，D为三维应力应变的材料本构矩阵；单元应变增量几何矩阵B＝[B1 B2 B3 B4 B5 B6 B7 B8]，节点纯变形q＝[qa
qb qc qd qe qf qg qh] ，其中qa,qb,qc...qg,qh分别代表单元节点a,b,c...g,h上的单元节点纯变形；
单元变形描述的虚功和应力应变描述的虚功相等原理等式为：式中，Va为六面体单元体积；
节点a,b,c...g,h虚拟位置下的单元节点内力；
式中，fi为实际位置下的单元节点内力；R′(θA)和R″(θB)分别为第一次和第二次正向运动转动矩阵，i＝a,b,c...g,h；
一种基于机械臂柔性变形的工业机器人动力学建模方法
技术领域
[0001] 本发明涉及一种工业机器人动力学建模方法，特别涉及一种基于机械臂柔性变形的工业机器人动力学建模方法，属于动力学建模技术领域。
背景技术
[0002] 随着对工业机器人精度要求的不断增加，工业机器人机械臂柔性变形导致的振动问题已无法忽略，传统消除工业机器人机械臂柔性变形的方法从提升机器人整体刚度的角
度着手，对机械臂材料和结构进行冗余设计，但是这种设计使得工业机器人产生平顺性差、耗能高和响应慢的问题，不得不进行轻量化、柔性悬挂设计，进而又会加剧工业机器人柔性变形的影响。
[0003] 工业机器人动力学建模能够对机械臂柔性变形机理进行分析，对进一步提高工业机器人定位精度和实时性能有着重要意义。传统的工业机器人动力学建模方法主要有
Lagrange、Newton‑Euler以及Kane法，多基于工业机器人机械臂刚性假设，而采用有限单元法为主的柔性变形分析方法，对于工业机器人这种刚体大转动‑多柔体变形耦合的模型，常面临求解效率低，刚度矩阵病态甚至无法求解的问题。
发明内容
[0004] 本发明的目的在于提供一种基于机械臂柔性变形的工业机器人动力学建模方法，使其能够简便、准确的建立工业机器人动力学模型，以解决机械臂柔性变形对末端精度的
影响，为工业机器人振动机理的研究以及提高定位精度，提供理论基础。
[0005] 本发明采用如下技术方案：
[0006] 一种基于机械臂柔性变形的工业机器人动力学建模方法，特征在于该方法包括以下步骤：
[0007] 步骤1，采用ABAQUS软件使用六面体单元网络建立工业机器人物理模型，模型所在T T T
全局坐标系O‑XYZ的坐标轴方向为：ex＝[1 0 0] ,ey＝[0 1 0] ,ez＝[0 0 1] ，获得模型网格参数，确定动力学建模参数和边界条件；
[0008] 步骤2，根据步骤1的模型网格参数，动力学建模参数以及边界条件，使用结合牛顿运动学和显式中心差分的控制方程求解机械臂柔性变形后的单元节点位置；
[0009]
[0010] 式中， 为第n+1步机械臂柔性变形后的单元节点i的位置； 为第n步机械臂柔性变形后的单元节点i的位置；n为计算步数；ξ为阻尼系数,h为计算步长,m为单元节点质
n n
量；Fi 为单元节点i的单元节点力， fi 为单元节点i的单元节点内力， 为单
元节点i的单元节点外力，i＝a,b,c...,g,h；
[0011] 步骤3，根据正运动学对单元节点位置 和 以关节为旋转中心旋转角度θ，获得旋转后的机械臂单元节点位置 和
[0012] 步骤4，以第n+1步的单元质心为参考点，将第n+1步的单元分别对应逆向运动至第n步的单元位置上，使得相邻两步的单元充分重合，进而根据经过逆向运动后虚拟单元节点位置与第n步的单元节点位置的差值计算出单元节点i的纯变形qi；
[0013] 步骤5，以第n步的单元质心为原点建立局部坐标系 坐标轴方向与全局坐标系O‑XYZ相同，将第n+1步单元节点位置和单元节点纯变形转换到局部坐标系中，采用八节点六面体单元形函数建立形函数坐标系O”‑rst，并将局部坐标系内单元节点位置通过雅克比矩阵转换到形函数坐标系中，根据形函数导数计算出单元应变增量，进而计算出单元
变形后的应力，利用单元变形描述的虚功和应力应变描述的虚功相等这一原理，计算出虚
拟单元节点位置下的单元节点内力，进一步对其进行正向运动，计算出实际位置下的单元
节点内力，根据惯性项计算出单元节点外力；
[0014] 步骤6，在第n步到第n+1步，经过旋转后的机械臂单元节点位置 和 分别作为下一计算步的单元节点位置 和
[0015] 更新计算时间time＝nh；
[0016] 步骤7，判断计算时间time是否达到最终时间ET；若否，转入到步骤2到步骤6继续迭代计算；若是，计算结束，输出计算结果。
[0017] 上述技术方案中，所述模型网格参数包括模型单元序号、单元节点序号和单元节点坐标值；所述动力学建模参数包括运动时间ET、材料密度rho、弹性模量E和泊松比pr；所述边界条件为对模型单元节点运动方向添加约束。
[0018] 进一步地，步骤3中所述正运动学方法的步骤为：
[0019] 建立关节转动轴线向量RA，确定回转中心点xcet；
[0020] 建立转动矩阵R：
[0021] R＝(1‑cosθ)·Rot2+sinθ·Rot
[0022] 式中，θ为旋转角度， RA＝[RAx RAy RAz]T；
[0023] 计算旋转后的单元节点位置：
[0024]
[0025] 进一步地，步骤4中逆向运动步骤为：逆向平移包括第一次逆向转动和第二次逆向转动；逆向平移是将第n+1步的单元，以单元质心 为参考点，整体逆向平移至第n步的单
元质心 上，经过逆向平移的单元节点位置x′i为：
[0026]
[0027] 第一次逆向转动首先对经过逆向平移的单元，以单元平面质心为参考点构造参考平面及其法向量vec′1，采用相同的方式对第n步的单元构造参考平面及其法向量vec1,进而计算出两法向量的夹角θA及其垂直向量nvecA:
[0028]
[0029]
[0030] 进而计算第一次逆向转动矩阵R′(‑θA):
[0031] R′(‑θA)＝[1‑cos(‑θA)]V′2+sin(‑θA)V′
[0032] 式中，
[0033] 经过第一次逆向转动的单元节点位置x″i为：
[0034] x″i＝R′(‑θA)x′i
[0035] 第二次逆向转动除选取的参考平面与第一次逆向转动不同外，其余步骤均相同，经过第二次逆向转动的单元节点位置x″′i为：
[0036] x″′i＝R″(‑θB)x″i
[0037] 式中，R″(‑θB)为第二次逆向转动矩阵，θB为参考平面法向量夹角；
[0038] 单元节点i的单元节点纯变形qi计算公式为：
[0039]
[0040] 进一步地，步骤5中将第n+1步单元节点位置转换到局部坐标系公式为：
[0041]
[0042] 式中， 为局部坐标系下的单元节点i的位置，
[0043] 八节点六面体单元形函数Ni表达式为：
[0044]
[0045] 式中，r、s、t分别为形函数表达式参数，ri＝±1，si＝±1，ti＝±1为形函数坐标系内节点i的坐标；
[0046] 雅克比矩阵J为：
[0047]
[0048] 形函数导数计算步骤为：对雅克比矩阵J求逆，建立形函数关于局部坐标系的导数与形函数关于形函数坐标的导数等式为：
[0049]
[0050] 单元变形后的应力 计算公式为：
[0051]
[0052] 式中， 为第n步单元的初始应力值，D为三维应力应变的材料本构矩阵；单元应变增量 几何矩阵B＝[B1 B2 B3 B4 B5 B6 B7 B8]，
[0053] j＝1,2,3...7,8；节点纯变形q＝[qa T
qb qc qd qe qf qg qh] ，其中qa,qb,qc...qg,qh分别代表单元节点a,b,c...g,h上的单元节点纯变形；
[0054] 单元变形描述的虚功和应力应变描述的虚功相等原理等式为：
[0055]
[0056] 式中，Va为六面体单元体积；
[0057] 虚拟位置下单元节点内力
[0058]
[0059] 式中， 其中 分别为单元节点a,b,c...g,h虚拟位置下的单元节点内力；
[0060] 正向运动的方式为：
[0061]
[0062] 式中，fin为实际位置下的单元节点内力；R′(θA)和R″(θB)分别为第一次和第二次正向运动转动矩阵，i＝a,b,c...g,h；
[0063] 单元节点外力pn计算公式为：
[0064]
[0065] 本发明与传统工业机器人动力学建模方法相比，其显著优点为：1)本发明对工业机器人建立的动力学模型，考虑了机械臂柔性变形对末端精度的影响，为工业机器人振动
机理的研究以及通过控制手段提高定位精度，提供了理论基础；2)本发明的工业机器人建
模流程简便快捷，对于不同的模型建模流程都是基于统一的框架，并且使用迭代计算的形
式，适用于计算机编程计算，有助于提高计算效率。
附图说明
[0066] 图1是一种基于机械臂柔性变形的工业机器人动力学建模方法流程图。
[0067] 图2是采用ABAQUS软件建立的工业机器人物理模型图。
[0068] 图3a是单个单元节点空间转动示意图。
[0069] 图3b是单个单元节点旋转平面示意图。
[0070] 图4是六面体单元运动和变形示意图。
[0071] 图5是逆向平移后的单元位置示意图。
[0072] 图6是第一次逆向转动六面体单元参考平面及其法向量示意图。
[0073] 图7是经第一次逆向转动后的单元位置示意图。
[0074] 图8是第二次逆向转动六面体单元参考平面及其法向量示意图。
[0075] 图9是经第二次逆向转动后的单元位置示意图。
[0076] 图10是六面体单元局部坐标系示意图。
[0077] 图11是八节点六面体单元形函数坐标示意图。
具体实施方式
[0078] 下面结合附图对本发明作进一步详细描述。
[0079] 参见图1，一种基于机械臂柔性变形的工业机器人动力学建模方法，其包括以下步骤：
[0080] 步骤1，采用ABAQUS软件使用六面体单元网络建立工业机器人物理模型，模型所在T T T
全局坐标系O‑XYZ的坐标轴方向为：ex＝[1 0 0] ,ey＝[0 1 0] ,ez＝[0 0 1] ，获得模型网格参数，包括模型单元序号、单元节点序号和单元节点坐标值。
[0081] 确定动力学建模参数，包括运动时间ET、材料密度rho、弹性模量E和泊松比pr。
[0082] 确定边界条件，为对模型单元节点运动方向添加约束。
[0083] 步骤2，根据步骤1的模型网格参数，动力学建模参数以及边界条件，使用结合牛顿运动学和显式中心差分的控制方程求解机械臂柔性变形后的单元节点位置。
[0084]
[0085] 式中， 为第n+1步机械臂柔性变形后的单元节点i的位置； 为第n步机械臂柔性变形后的单元节点i的位置；n为计算步数；ξ为阻尼系数,h为计算步长,m为单元节点质
n n
量；Fi 为单元节点i的单元节点力， fi为单元节点i的单元节点内力， 为单
元节点i的单元节点外力，i＝a,b,c...,g,h。
[0086] 步骤3，根据正运动学对第n+1步机械臂柔性变形后的单元节点位置 和第n步机械臂柔性变形前的单元节点位置 以关节为旋转中心旋转角度θ，获得旋转后的机械臂
单元节点位置 和 如图3a所示，以工业机器人模型中六面体单元一个节点A,位置向
量xp为例，进行推导正运动学公式，单元节点A,以关节轴线作为旋转中心，旋转向量θ＝θeθ,eθ为转轴的方向向量，θ为关节转角，转角与转轴方向均以右旋规则为正向。旋转后的位置为B,关节轴线与旋转平面的交点为C,则构成与旋转轴线垂直的旋转平面ABC，如图3b所示，其中A点和B点是在一个以C点为圆心，r为半径的圆弧上。建立圆弧上A点的单位切线向量
et。
[0087]
[0088] 式中，α是eA和eθ之间的夹角，eA为位置向量xp的方向向量。
[0089] 令es为圆弧在A点沿半径方向的向量：
[0090] es＝eθ×et (3)
[0091] 将(eθ,et,es)形成一组基底，空间转动产生的位移向量η可以写成：
[0092]
[0093] 式中半径r＝|eθ×xp|＝xpsinα。
[0094] 令R*＝(1‑cosθ)A2+sinθA，则空间转动产生的位移向量η可以改写为：
[0095] η＝R*xp (4)
[0096] 经过旋转后节点B的位置为：
[0097] xB＝xp+η＝(I+R*)xp (5)
[0098] 根据以上正运动学求解工业机器人模型步骤为：
[0099] 建立关节转动轴线向量RA；确定回转中心点xcet。
[0100] 建立转动矩阵R：
[0101] R＝(1‑cosθ)·Rot2+sinθ·Rot (6)T
[0102] 式中，θ为旋转角度， RA＝[RAx RAy RAz]。
[0103] 计算旋转后的单元节点位置：
[0104]
[0105] 步骤4，以第n步到第n+1步经过运动和变形的单个六面体单元为例：如图4所示，在全局坐标系中，第n步对应的六面体单元各节点记作：an‑bn‑cn‑dn‑en‑fn‑gn‑hn,节点i的位置矢量记作 单元质心位置矢量记作 则第n+1步对应的六面体单元各节点记作：an+1‑bn+1‑cn+1‑dn+1‑en+1‑fn+1‑gn+1‑hn+1,节点i的位置矢量记作 单元质心位置矢量记作[0106] 逆向平移将第n+1步的单元，以单元质心 为参考点，整体逆向平移至第n步的单
元质心 上，经过逆向平移的单元如图5所示，单元节点记作：a′‑b′‑c′‑d′‑e′‑f′‑g′‑h′，单元节点位置x′i为：
[0107]
[0108] 如图6所示，第一次逆向转动首先对经过逆向平移的单元，以单元平面质心为参考点构造参考平面1′‑4′‑6'及其法向量vec′1，采用相同的方式对第n步的单元构造参考平面
1‑4‑6及其法向量vec1,进而计算出两法向量的夹角θA及其垂直向量nvecA:
[0109]
[0110]
[0111] 进而计算第一次逆向转动矩阵R′(‑θA):
[0112] R′(‑θA)＝[1‑cos(‑θA)]V′2+sin(‑θA)V′ (11)
[0113] 式中，
[0114] 以第n步单元质心 为旋转中心对各节点a′‑b′‑c′‑d′‑e′‑f′‑g′‑h′进行逆向旋转，如图7所示，获得旋转后的单元节点位置：a″‑b″‑c″‑d″‑e″‑f″‑g″‑h″，经过第一次逆向转动的单元节点位置x″i为：
[0115] x″i＝R′(‑θA)x′i (12)
[0116] 第二次逆向转动对经过第一次逆向旋转的单元，如图8所示，以其余三个平面质心为参考点建立参考平面2′‑3′‑5'，其法向量为vec′2,按照相同的方式，对第n步的单元构造参考平面2‑3‑5，其法向量为vec2。第二次逆向旋转步骤与第一次逆向旋转相同。如图9所示，经过第二次逆向转动的单元节点位置为：a″′‑b″′‑c″′‑d″′‑e″′‑f″′‑g″′‑h″′,单元节点位置x″′i为：
[0117] x″′i＝R″(‑θB)x″i (13)
[0118] 式中，R″(‑θB)为第二次逆向转动矩阵，θB为参考平面法向量夹角；
[0119] 进而根据经过逆向运动后虚拟单元节点位置与第n步的单元节点位置的差值计算出单元节点i的纯变形qi：
[0120]
[0121] 步骤5，如图10所示，以第n步的单元质心为原点建立局部坐标系 坐标轴方向与全局坐标系O‑XYZ相同，将第n+1步单元节点位置 转换到局部坐标系中，局部坐标
系下的单元节点位置
[0122]
[0123] 式中，
[0124] 如图11所示，采用八节点六面体单元形函数，对第n步的六面体单元建立O”‑rst形函数坐标系。
[0125] 对应的八节点六面体单元的形函数Ni表达式为：
[0126]
[0127] 式中，r，s，t分别为形函数表达式参数，ri＝±1，si＝±1，ti＝±1为形函数坐标系内单元节点i的坐标。
[0128] 进一步对局部坐标系内单元节点位置 通过雅克比矩阵J转换到形函数坐标系中，雅克比矩阵J为：
[0129]
[0130] 对雅克比矩阵J求逆，建立形函数关于局部坐标系的导数与形函数关于形函数坐标的导数等式为：
[0131]
[0132] 根据形函数导数计算出单元应变增量，进而计算出单元变形后的应力，单元变形后的应力 计算公式为：
[0133]
[0134] 式中， 为第n步单元的初始应力值，D为三维应力应变的材料本构矩阵；单元应变增量 几何矩阵B＝[B1 B2 B3 B4 B5 B6 B7 B8]，
[0135] j＝1,2,3...7,8；节点纯变形q＝[qa T
qb qc qd qe qf qg qh] ，其中qa,qb,qc...qg,qh分别代表单元节点a,b,c...g,h上的单元节点纯变形。
[0136] 利用单元变形描述的虚功和应力应变描述的虚功相等这一原理建立等式：
[0137]
[0138] 式中，Va为六面体单元体积；
[0139] 进一步计算出节点内力
[0140]
[0141] 式中， 其中 分别为单元节点a,b,c...g,h虚拟位置下的单元节点内力。
[0142] 对虚拟单元节点位置下的单元节点内力进行正向运动，计算出实际位置下的单元n
节点内力fi:
[0143]
[0144] 式中，R′(θA)和R″(θB)分别为第一次和第二次正向运动转动矩阵。
[0145] 单元节点外力pn计算公式为：
[0146]
[0147] 步骤6，在第n到第n+1步，经过旋转后的机械臂单元节点位置 和 分别作为下一计算步的单元节点位置 和
[0148] 更新计算时间time＝nh；
[0149] 步骤7，判断计算时间time是否达到最终时间ET；若否，转入到步骤2到步骤6继续迭代计算；若是，计算结束，输出计算结果；
[0150] 下面结合实施例对本发明做进一步的详细描述。
[0151] 实施例：
[0152] 建立广数RB08A3工业机器人模型，模型结构参数如表1所示：
[0153] 表1广数RB08A3工业机器人结构参数
[0154]
[0155] 工业机器人动力学模型所需参数如表2所示：
[0156] 表2工业机器人动力学模型参数
[0157]
[0158]
[0159] 上述模型运动状态为J2和J3关节匀速转动，其余关节保持固定，使用本发明所述的一种基于机械臂柔性变形的工业机器人动力学建模方法，基于MATLAB对上述工业机器人
模型编程，建立其动力学模型，并记录每一步单元节点位置，利用差分的形式计算出工业机器人末端加速度，加速度的差分形式为:
[0160]
[0161] 基于上述参数对工业机器人进行实验，利用加速度传感器测得末端振动加速度，实验结果出现特征频率14Hz和164Hz，而建模结果具有13.5Hz和166Hz，相对实验结果的误
差分别为‑3.6％和1.2％。
