CN112001087B - 一种旋转关节型工业机器人非线性动力学建模分析方法 - Google Patents

一种旋转关节型工业机器人非线性动力学建模分析方法 Download PDF

Info

Publication number
CN112001087B
CN112001087B CN202010875365.4A CN202010875365A CN112001087B CN 112001087 B CN112001087 B CN 112001087B CN 202010875365 A CN202010875365 A CN 202010875365A CN 112001087 B CN112001087 B CN 112001087B
Authority
CN
China
Prior art keywords
point
industrial robot
formula
matrix
hinge
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010875365.4A
Other languages
English (en)
Other versions
CN112001087A (zh
Inventor
李波
田威
崔光裕
廖文和
李宇飞
张霖
胡俊山
石瀚斌
梁爽
王俊
张苇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Changjiang Industrial Technology Research Institute Co ltd
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202010875365.4A priority Critical patent/CN112001087B/zh
Publication of CN112001087A publication Critical patent/CN112001087A/zh
Application granted granted Critical
Publication of CN112001087B publication Critical patent/CN112001087B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Manipulator (AREA)

Abstract

本发明公开了一种旋转关节型工业机器人非线性动力学建模分析方法，包括如下步骤：将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的系统动力学模型；选取工业机器人各关节转角值作为系统广义坐标，确定积分变量x(ti)及其初始条件x(t0)；计算工业机器人体元件和铰元件的传递矩阵和传递方程；根据元件传递矩阵和机器人系统拓扑结构，计算总传递矩阵，结合边界条件，求解总传递方程，获得边界点状态矢量中的未知状态变量；根据动力学模型中各连接点处的加速度、角加速度，对动力学模型进行运动学分析，获得广义坐标的二阶导数；结合数值积分方法，获得ti+1时刻的积分变量x(ti+1)，当到达期望的计算时间则结束，否则返回重新确定初始条件并再次分析计算。

Description

一种旋转关节型工业机器人非线性动力学建模分析方法
技术领域
本发明属于工业机器人动力学建模技术领域，具体涉及一种旋转关节型工业机器人非线性动力学建模分析方法。
背景技术
工业机器人在机械加工、汽车制造、航空航天制造等自动化生产领域的广泛应用，对其运动控制精度和动态特性等核心性能提出了更高要求。传统的基于运动学的控制与设计研究已经无法满足应用需求，必须从动力学层面入手，通过对机器人加工系统进行动力学研究，获取机器人运动过程中的动力响应，作为机器人高精度控制和优化设计的基础。机器人的动力学研究对于提高机器人运动精度与动态特性至关重要，一方面必须建立能够准确描述机器人系统动态特性的动力学方程，另一方面在保证动力学建模精度的前提下实现动力学快速计算，满足实时控制要求。
严格来说工业机器人动力学系统是一个非线性多体系统，其动力学方程相当复杂。动力学研究方法主要有牛顿-欧拉法、拉格朗日法、约丹速度变分法、Kane法、Wittenburg法、Schiehlen法、Shabana法等。已公开发明专利“一种工业机器人整体动力学建模及动力学参数辨识方法，CN 110539302 A”采用拉格朗日方法建立工业机器人整体动力学模型，并进行了动力学参数辨识。现有技术存在的问题是：现有动力学建模方法在计算工业机器人动力学响应时，需要建立系统的总体动力学方程，虽然能够较为准确地描述工业机器人的实际动态特性，但其系统矩阵阶次随着自由度的增加不断提高，导致计算速度下降、难以满足机器人运动过程中的实时计算与控制需求。
发明内容
本发明所要解决的技术问题是针对上述现有技术的不足，提供一种旋转关节型工业机器人非线性动力学建模分析方法。
为实现上述技术目的，本发明采取的技术方案为：
一种旋转关节型工业机器人非线性动力学建模分析方法，其中：包括如下步骤：
步骤S1：将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型，并建立拓扑结构图；
步骤S2：选取工业机器人各关节转角值作为系统广义坐标，确定积分变量x(ti)及其初始参数；
步骤S3：对于ti时刻，根据积分变量进行对铰元件和体元件进行运动学分析，得到各体元件的位置、姿态、速度、角速度变量；
步骤S4：计算工业机器人体元件和铰元件的传递矩阵和传递方程；
步骤S5：根据元件传递矩阵和机器人系统拓扑结构，计算总传递矩阵，结合边界条件，求解总传递方程，获得边界点状态矢量中的未知状态变量；
步骤S6：根据动力学模型中各连接点处的加速度、角加速度，对动力学模型进行运动学分析，获得广义坐标的二阶导数，结合数值积分方法，获得ti+1时刻的积分变量x(ti+1)，当到达期望的计算时间，则结束，否则返回步骤S3继续计算。
为优化上述技术方案，采取的具体措施还包括：
进一步地，步骤S1具体为：
S11：根据旋转型工业机器人的几何结构、装配关系，以及各个部件间的连接方式，将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型；
S12：将工业机器人的连杆和连接件等价于体元件，将工业机器人用于连接的关节等价于铰元件，其中体元件包含质量，铰元件不包含质量；
S13：将工业机器人的基座边界编号为0，铰元件编号为奇数，体元件编号为偶数。
进一步地，步骤S2具体为：选取工业机器人各关节转角值作为系统广义坐标，即q＝[θ1 θ3 θ5 θ7 θ9 θ11]T，确定积分变量
Figure GDA0002705925310000021
和每个铰元件的初始转动角度和角速度。
进一步地，步骤S3具体为：
S31：以工业机器人的基座为原点建立oxyz惯性坐标系，以体元件的力的输入点I为中心建立IxIyIzI连体坐标系，建立各部件连接点的状态矢量z为
Figure GDA0002705925310000022
其中，
Figure GDA0002705925310000023
为连接点的绝对加速度在惯性坐标系oxyz的坐标列阵，
Figure GDA0002705925310000024
为连接点处对应坐标系的绝对角加速度在惯性系oxyz的坐标列阵，m＝[mx my mz]T和q＝[qx qy qz]T分别为连接点处内力矩、内力在惯性系oxyz的坐标列阵；
S32：对体元件进行运动学分析
由体元件的输入点I和输出点O的绝对角速度、角加速度相等得：
Figure GDA00027059253100000312
式中，ΩI、
Figure GDA0002705925310000032
表示输入点I的绝对角速度和绝对角加速度，ΩO、
Figure GDA0002705925310000033
表示输出点O的绝对角速度和绝对角加速度；
输出点O相对于惯性坐标系原点的位置矢量在惯性坐标系中的投影rO表示为
rO＝rI+rIO＝rI+AIlIO (3)
ΩI＝AIωI
式中，rI为输入点I相对于惯性坐标系原点的位置矢量，rIO和lIO分别表示输出点相对于输入点的位置矢量在惯性坐标系和连体坐标系IxIyIzI的投影，AI为输入点连体系IxIyIzI相对于惯性系oxyz的方向余弦矩阵，ΩI、ωI分别为输入点连体坐标系在惯性坐标系中的角速度矢量在惯性坐标系oxyz和输入端连体坐标系IxIyIzI中的投影；
计算体元件的速度
Figure GDA0002705925310000034
S33：对铰元件进行动力学分析，则铰元件输入点与输出点状态矢量为
Figure GDA0002705925310000035
式中，Aj,O为铰元件输出点连体系相对于惯性坐标系的方向余弦矩阵，
Figure GDA0002705925310000036
表示铰元件外接刚体相对于内接刚体的角速度。
进一步地，步骤S4具体为：对式(3)关于时间求导可得
Figure GDA0002705925310000037
Figure GDA0002705925310000038
Figure GDA0002705925310000039
Figure GDA00027059253100000310
其中
Figure GDA00027059253100000311
表示任一向量a对应的反对称矩阵；
将式(9)代入式(6)可得
Figure GDA0002705925310000041
同理，对于质心C，有
Figure GDA0002705925310000042
根据刚体质心运动定律，可得体元件随输入点平动的动力学方程为
Figure GDA0002705925310000043
式中，m为体元件质量，qI,qO分别为输入点I、输出点O的内力，fC为作用在质心C处的外力；
根据体元件相对于任意动点的相对动量矩定理，选择输入点I为动点，并投影至惯性坐标系，可得
Figure GDA0002705925310000044
式中，JI为体元件相对于输入点的转动惯量在输入点连体坐标系IxIyIzI中的矩阵表示，mC表示作用在体元件质心处的外力矩；
将式(12)代入式(13)，并整理式(10)、(4)、(13)、(12)可得
Figure GDA0002705925310000045
将式(14)写成矩阵的形式
zO＝UzI (15)
可得体元件的传递矩阵为
Figure GDA0002705925310000046
式中，I3为3阶单位矩阵，O3为3×3全零矩阵，O3×1为3×1全零矩阵，O1×3为1×3全零矩阵；
得到各体元件的传递方程为
Figure GDA0002705925310000051
式中，Ui(i＝2,4,6,8,10,12)为体元件的传递矩阵，z2,1,z2,3,L,z12,11为工业机器人系统内各连接点状态矢量，z12,0为工业机器人系统自由端边界点状态矢量；
S42：计算铰元件的传递矩阵和传递方程
由铰元件的输入点与输出点状态矢量式(5)第四项等式对时间求一阶导数可得
Figure GDA0002705925310000052
式(18)两边同时乘以
Figure GDA0002705925310000053
得
Figure GDA0002705925310000054
式中，
Figure GDA0002705925310000055
Figure GDA0002705925310000056
为方向余弦矩阵Aj,O的转置。
铰元件绕z轴方向的内力矩
Figure GDA0002705925310000057
可得
Figure GDA0002705925310000058
式中，H2＝[0 0 1]。
当外接体元件j+1后再接一个铰元件，可得
Figure GDA0002705925310000059
由于外接体元件j+1的传递方程和传递矩阵分别为
zj+1,O＝Uj+1zj+1,I (22)
Figure GDA00027059253100000510
式中，u3,1、u3,2、u3,3、u3,4、u3,5为铰元件j外接体元件j+1传递矩阵的分块矩阵中的元素。
则有
Figure GDA0002705925310000061
等式两边同时左乘
Figure GDA0002705925310000062
并结合式(21)，整理可得
Figure GDA0002705925310000063
由于第j个铰元件输出点状态矢量与第j+1个体元件输入点状态矢量相等，即
Figure GDA00027059253100000610
将式(26)代入式(25)，得
Figure GDA0002705925310000065
联立式(19)与式(27)，可得铰元件j的传递矩阵为
Figure GDA0002705925310000066
Figure GDA0002705925310000067
Figure GDA0002705925310000068
则各铰元件的传递方程为
Figure GDA0002705925310000069
式中，Uj(j＝1,3,5,7,9,11)为铰元件的传递矩阵，z1,0为工业机器人系统固定端边界点状态矢量。
进一步地，步骤S5具体为：
根据各体元件传递方程(17)和各铰元件传递方程(29)，得到机器人系统总传递方程
z12,0＝U1-12z1,0 (30)
式中，U1-12＝U12…U2U1为系统总传递矩阵；
工业机器人固定端边界点和自由端边界点的边界条件为
Figure GDA0002705925310000071
将边界条件(31)代入式(30)，得到分块矩阵形式为
Figure GDA0002705925310000072
Figure GDA0002705925310000073
优化式(33)，计算边界点状态矢量中的未知状态变量
Figure GDA0002705925310000074
根据计算得到的未知状态变量获得系统固定端状态矢量，然后依次利用各元件的传递方程(17)和(29)，获得系统中各连接点的状态矢量。
进一步地，所述步骤S6中铰元件θj(j＝1,3,5,7,9,11)的广义坐标的二阶导数的计算过程为：
在式(18)两边同时左乘
Figure GDA0002705925310000075
得到
Figure GDA0002705925310000076
本发明的有益效果：
本发明提供了一种旋转关节型工业机器人非线性动力学建模分析方法，无需系统总体动力学方程、系统矩阵阶次不随系统自由度的增加而增加，可大幅度提高机器人空间大运动非线性动力学建模与计算效率。
附图说明
图1是本发明的旋转关节型工业机器人多体系统动力学模型示意图；
图2是本发明的旋转关节型工业机器人动力学模型拓扑结构示意图；
图3是本发明的空间运动刚体及其坐标系示意图；
图4是本发明的空间运动柱铰及其坐标系示意图；
图5是本发明的旋转关节型工业机器人多体系统动力学数值计算流程图。
具体实施方式
以下结合附图对本发明的实施例作进一步详细描述，本发明采用KUKA KR500-2830MT型工业机器人多体系统动力学模型。
如图1所示，本发明为一种旋转关节型工业机器人非线性动力学建模分析方法，其中：包括如下步骤：
步骤S1：将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型，并建立拓扑结构图；
S11：根据旋转型工业机器人的几何结构、装配关系，以及各个部件间的连接方式，将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型；
S12：将工业机器人的连杆和连接件等价于体元件，将工业机器人用于连接的关节等价于铰元件，其中体元件包含质量，铰元件不包含质量；
S13：将工业机器人的基座边界编号为0，铰元件编号为奇数，体元件编号为偶数。
机器人地基编号为0；6个关节依次编号为1、3、5、7、9、11；6个连杆依次编号为2、4、6、8、10、12。需要指出的是，KUKA机器人的后三个电机在第3连杆上，因此体元件6的质量包括第3连杆及后三个电机伺服机构的质量。对于其他型号的机器人，只需根据其电机的具体布局方式，选择对应的铰接形式即可。
将机器人各连杆视为空间运动刚体，各关节视为空间光滑柱铰，电机产生的柱铰轴方向的控制力矩分别作用在相邻的两个连杆上，且锁止机构工作时的弹性阻尼作用等效为与柱铰轴方向并联的扭簧和旋转阻尼器。将机器人加工时外界对末端执行器的力处理为作用在机器人系统上的外力。
如图2所示，KUKA KR500-2830MT型工业机器人多体动力学模型即为：在电机驱动力和外界负载加工力作用下的由6个空间运动刚体和6个空间柱铰组成的多刚体系统。根据KUKA KR500-2830MT型工业机器人多体系统动力学模型，建立其动力学模型拓扑结构，其中封闭曲线代表体元件；箭头代表铰元件，箭头方向代表机器人系统连接点状态矢量的传递方向；数字为元件编号；边界统一编号为0。
步骤S2：选取工业机器人各关节转角值作为系统广义坐标，确定积分变量x(ti)及其初始参数；
步骤S2具体为：选取工业机器人各关节转角值作为系统广义坐标，即q＝[θ1 θ3 θ5θ7 θ9 θ11]T，确定积分变量
Figure GDA0002705925310000091
和每个铰元件的初始转动角度和角速度。
步骤S3：于ti时刻，根据积分变量进行对铰元件和体元件进行运动学分析，得到各体元件的位置、姿态、速度、角速度变量。
步骤S3具体为：
S31：以工业机器人的基座为原点建立oxyz惯性坐标系，以体元件的力的输入点I为中心建立IxIyIzI连体坐标系，建立各部件连接点的状态矢量z为
Figure GDA0002705925310000092
其中，
Figure GDA0002705925310000093
为连接点的绝对加速度在惯性坐标系oxyz的坐标列阵，
Figure GDA0002705925310000094
为连接点处对应坐标系的绝对角加速度在惯性系oxyz的坐标列阵，m＝[mx my mz]T和q＝[qx qy qz]T分别为连接点处内力矩、内力在惯性系oxyz的坐标列阵；
S32：对体元件进行运动学分析
由体元件的输入点I和输出点O的绝对角速度、角加速度相等得：
Figure GDA0002705925310000098
式中，ΩI、
Figure GDA0002705925310000096
表示输入点I的绝对角速度和绝对角加速度，ΩO、
Figure GDA0002705925310000097
表示输出点O的绝对角速度和绝对角加速度；
输出点O相对于惯性坐标系原点的位置矢量在惯性坐标系中的投影rO表示为
rO＝rI+rIO＝rI+AIlIO (3)
ΩI＝AIωI
式中，rI为输入点I相对于惯性坐标系原点的位置矢量，rIO和lIO分别表示输出点相对于输入点的位置矢量在惯性坐标系和连体坐标系IxIyIzI的投影，AI为输入点连体系IxIyIzI相对于惯性系oxyz的方向余弦矩阵，ΩI、ωI分别为输入端连体坐标系在惯性坐标系中的角速度矢量在惯性坐标系oxyz和输入端连体坐标系IxIyIzI中的投影；
计算体元件的速度
Figure GDA0002705925310000101
S33：对铰元件进行动力学分析，则铰元件输入点与输出点状态矢量为
Figure GDA0002705925310000102
式中，Aj,O为铰元件输出点连体系相对于惯性坐标系的方向余弦矩阵，
Figure GDA0002705925310000103
表示铰元件外接刚体相对于内接刚体的角速度。
步骤S4：计算工业机器人体元件和铰元件的传递矩阵和传递方程；
步骤S4具体为：
S41：计算体元件的传递矩阵和传递方程
对式(3)关于时间求导可得
Figure GDA0002705925310000104
Figure GDA0002705925310000105
Figure GDA0002705925310000106
Figure GDA0002705925310000107
其中
Figure GDA0002705925310000108
表示任一向量a对应的反对称矩阵；
将式(9)代入式(6)可得
Figure GDA0002705925310000109
同理，对于质心C，有
Figure GDA00027059253100001010
根据刚体质心运动定律，可得体元件随输入点平动的动力学方程为
Figure GDA00027059253100001011
式中，m为体元件质量，qI,qO分别为输入点I、输出点O的内力，fC为作用在质心C处的外力；
根据体元件相对于任意动点的相对动量矩定理，选择输入点I为动点，并投影至惯性坐标系，可得
Figure GDA0002705925310000111
式中，JI为体元件相对于输入点的转动惯量在输入点连体坐标系IxIyIzI中的矩阵表示，mC表示作用在体元件质心处的外力矩；
将式(12)代入式(13)，并整理式(10)、(4)、(13)、(12)可得
Figure GDA0002705925310000112
将式(14)写成矩阵的形式
zO＝UzI (15)
可得体元件的传递矩阵为
Figure GDA0002705925310000113
式中，I3为3阶单位矩阵，O3为3×3全零矩阵，O3×1为3×1全零矩阵，O1×3为1×3全零矩阵；
得到各体元件的传递方程为
Figure GDA0002705925310000114
式中，Ui(i＝2,4,6,8,10,12)为体元件的传递矩阵，z2,1,z2,3,L,z12,11为工业机器人系统内各连接点状态矢量，z12,0为工业机器人系统自由端边界点状态矢量；
S42：计算铰元件的传递矩阵和传递方程
由铰元件的输入点与输出点状态矢量式(5)第四项等式对时间求一阶导数可得
Figure GDA0002705925310000115
式(18)两边同时乘以
Figure GDA0002705925310000116
得
Figure GDA0002705925310000121
式中，
Figure GDA0002705925310000122
Figure GDA0002705925310000123
为方向余弦矩阵Aj,O的转置；
铰元件绕z轴方向的内力矩
Figure GDA0002705925310000124
可得
Figure GDA0002705925310000125
式中，H2＝[0 0 1]。
当外接体元件j+1后再接一个铰元件，可得
Figure GDA0002705925310000126
由于外接体元件j+1的传递方程和传递矩阵分别为
zj+1,O＝Uj+1zj+1,I (22)
Figure GDA0002705925310000127
式中，u3,1、u3,2、u3,3、u3,4、u3,5为铰元件j外接体元件j+1传递矩阵的分块矩阵中的元素；
则有
Figure GDA0002705925310000128
等式两边同时左乘
Figure GDA0002705925310000129
并结合式(21)，整理可得
Figure GDA00027059253100001210
由于第j个铰元件输出点状态矢量与第j+1个体元件输入点状态矢量相等，即
Figure GDA00027059253100001213
将式(26)代入式(25)，得
Figure GDA00027059253100001212
联立式(19)与式(27)，可得铰元件j的传递矩阵为
Figure GDA0002705925310000131
Figure GDA0002705925310000132
Figure GDA0002705925310000133
则各铰元件的传递方程为
Figure GDA0002705925310000134
式中，Uj(j＝1,3,5,7,9,11)为铰元件的传递矩阵，z1,0为工业机器人系统固定端边界点状态矢量。
步骤S5：根据元件传递矩阵和机器人系统拓扑结构，计算总传递矩阵，结合边界条件，求解总传递方程，获得边界点状态矢量中的未知状态变量；
步骤S5具体为:
根据各体元件传递方程(17)和各铰元件传递方程(29)，得到机器人系统总传递方程
z12,0＝U1-12z1,0 (30)
式中，U1-12＝U12…U2U1为系统总传递矩阵；
工业机器人固定端边界点和自由端边界点的边界条件为
Figure GDA0002705925310000135
将边界条件(31)代入式(30)，得到分块矩阵形式为
Figure GDA0002705925310000141
Figure GDA0002705925310000142
优化式(33)，计算边界点状态矢量中的未知状态变量
Figure GDA0002705925310000143
根据计算得到的未知状态变量获得系统固定端状态矢量，然后依次利用各元件的传递方程(17)和(29)，获得系统中各连接点的状态矢量。
步骤S5：根据动力学模型中各连接点处的加速度、角加速度，对动力学模型进行运动学分析，获得广义坐标的二阶导数；
所述步骤S6中铰元件θj(j＝1,3,5,7,9,11)的广义坐标的二阶导数的计算过程为：在式(18)两边同时左乘
Figure GDA0002705925310000144
得到
Figure GDA0002705925310000145
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (6)

1.一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于，包括如下步骤：
步骤S1：将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型，并建立拓扑结构图；
步骤S2：选取工业机器人各铰元件作为系统广义坐标，确定积分变量x(ti)及其初始条件；
步骤S3：对于ti时刻，根据积分变量进行对铰元件和体元件进行运动学分析，得到各体元件的位置、姿态、速度、角速度变量；具体包括如下子步骤：
S31：以工业机器人的基座为原点建立oxyz惯性坐标系，以体元件的力的输入点I为中心建立IxIyIzI连体坐标系，建立各部件连接点的状态矢量z为
Figure FDA0003249931780000011
其中，
Figure FDA0003249931780000012
为连接点的绝对加速度在惯性坐标系oxyz的坐标列阵，
Figure FDA0003249931780000013
为连接点处对应坐标系的绝对角加速度在惯性系oxyz的坐标列阵，m＝[mx my mz]T和q＝[qx qy qz]T分别为连接点处内力矩、内力在惯性系oxyz的坐标列阵；
S32：对体元件进行运动学分析
由体元件的输入点I和输出点O的绝对角速度、角加速度相等得：
Figure FDA0003249931780000014
式中，ΩI、
Figure FDA0003249931780000015
表示输入点I的绝对角速度和绝对角加速度，ΩO、
Figure FDA0003249931780000016
表示输出点O的绝对角速度和绝对角加速度；
输出点O相对于惯性坐标系原点的位置矢量在惯性坐标系中的投影rO表示为
rO＝rI+rIO＝rI+AIlIO (3)
ΩI＝AIωI
式中，rI为输入点I相对于惯性坐标系原点的位置矢量，rIO和lIO分别表示输出点相对于输入点的位置矢量在惯性坐标系和连体坐标系IxIyIzI的投影，AI为输入点连体系IxIyIzI相对于惯性系oxyz的方向余弦矩阵，ΩI、ωI分别为输入端连体坐标系在惯性坐标系中的角速度矢量在惯性坐标系oxyz和输入端连体坐标系IxIyIzI中的投影；
计算体元件的速度
Figure FDA0003249931780000021
式中，
Figure FDA0003249931780000022
表示输出点O的速度，
Figure FDA0003249931780000023
表示输入点I的速度，
Figure FDA0003249931780000024
表示AI对时间的一阶导数；
S33：对铰元件进行动力学分析，则铰元件输入点与输出点状态矢量为
Figure FDA0003249931780000025
式中，
Figure FDA0003249931780000026
mj,O、qj,O、Ωj,O分别为第j个铰元件输出点的绝对加速度、内力矩、内力和绝对角速度的坐标列阵；
Figure FDA0003249931780000027
mj,I、qj,I、Ωj,I分别为第j个铰元件输入点的绝对加速度、内力矩、内力和绝对角速度的坐标列阵；Aj,O为铰元件输出点连体系相对于惯性坐标系的方向余弦矩阵，
Figure FDA0003249931780000028
表示铰元件外接刚体相对于内接刚体的角速度，θj为旋转角度；
步骤S4：计算工业机器人体元件和铰元件的传递矩阵和传递方程；
步骤S5：根据元件传递矩阵和机器人系统拓扑结构，计算总传递矩阵，结合边界条件，求解总传递方程，获得边界点状态矢量中的未知状态变量；
步骤S6：根据动力学模型中各连接点处的加速度、角加速度，对动力学模型进行运动学分析，获得广义坐标的二阶导数,结合数值积分方法，获得ti+1时刻的积分变量x(ti+1)，当到达期望的计算时间，则结束，否则返回步骤S3继续计算。
2.根据权利要求1所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S1具体为：
S11：根据旋转型工业机器人的几何结构、装配关系，以及各个部件间的连接方式，将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型；
S12：将工业机器人的连杆和连接件等价于体元件，将工业机器人用于连接的关节等价于铰元件，其中体元件包含质量，铰元件不包含质量；
S13：将工业机器人的基座边界编号为0，铰元件编号为奇数，体元件编号为偶数。
3.根据权利要求1所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S2具体为：选取工业机器人各关节转角值作为系统广义坐标，即q＝[θ1θ3 θ5 θ7 θ9 θ11]T，确定积分变量
Figure FDA0003249931780000029
和每个铰元件的初始转动角度和角速度。
4.根据权利要求1所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S4具体为：
S41：计算体元件的传递矩阵和传递方程
对式(3)关于时间求导可得
Figure FDA0003249931780000031
Figure FDA0003249931780000032
Figure FDA0003249931780000033
Figure FDA0003249931780000034
其中
Figure FDA0003249931780000035
表示任一向量a对应的反对称矩阵；
将式(9)代入式(6)可得
Figure FDA0003249931780000036
同理，对于质心C，有
Figure FDA0003249931780000037
根据刚体质心运动定律，可得体元件随输入点平动的动力学方程为
Figure FDA0003249931780000038
式中，m为体元件质量，qI,qO分别为输入点I、输出点O的内力，fC为作用在质心C处的外力；
根据体元件相对于任意动点的相对动量矩定理，选择输入点I为动点，并投影至惯性坐标系，可得
Figure FDA0003249931780000039
式中，JI为体元件相对于输入点的转动惯量在输入点连体坐标系IxIyIzI中的矩阵表示，mC表示作用在体元件质心处的外力矩；
将式(12)代入式(13)，并整理式(10)、(4)、(13)、(12)可得
Figure FDA00032499317800000310
将式(14)写成矩阵的形式
zO＝UzI (15)
可得体元件的传递矩阵为
Figure FDA0003249931780000041
式中，I3为3阶单位矩阵，O3为3×3全零矩阵，O3×1为3×1全零矩阵，O1×3为1×3全零矩阵；
得到各体元件的传递方程为
Figure FDA0003249931780000042
式中，Ui(i＝2,4,6,8,10,12)为体元件的传递矩阵，z2,1,z2,3,…,z12,11为工业机器人系统内各连接点状态矢量，z12,0为工业机器人系统自由端边界点状态矢量；
S42：计算铰元件的传递矩阵和传递方程
由铰元件的输入点与输出点状态矢量式(5)中第四项等式对时间求一阶导数可得
Figure FDA0003249931780000043
式(18)两边同时乘以
Figure FDA0003249931780000044
得
Figure FDA0003249931780000045
式中，
Figure FDA0003249931780000046
Figure FDA0003249931780000047
为方向余弦矩阵Aj,O的转置；
铰元件绕z轴方向的内力矩
Figure FDA0003249931780000048
可得
Figure FDA0003249931780000049
式中，H2＝[0 0 1]；
当外接体元件j+1后再接一个铰元件，可得
Figure FDA00032499317800000410
由于外接体元件j+1的传递方程和传递矩阵分别为
zj+1,O＝Uj+1zj+1,I (22)
Figure FDA0003249931780000051
式中，u3,1、u3,2、u3,3、u3,4、u3,5为铰元件j外接体元件j+1传递矩阵的分块矩阵中的元素；
则有
Figure FDA0003249931780000052
等式两边同时左乘
Figure FDA0003249931780000053
并结合式(21)，整理可得
Figure FDA0003249931780000054
由于第j个铰元件输出点状态矢量与第j+1个体元件输入点状态矢量相等，即
Figure FDA0003249931780000055
将式(26)代入式(25)，得
Figure FDA0003249931780000056
联立式(19)与式(27)，可得铰元件j的传递矩阵为
Figure FDA0003249931780000057
Figure FDA0003249931780000058
Figure FDA0003249931780000059
则各铰元件的传递方程为
Figure FDA00032499317800000510
式中，Uj(j＝1,3,5,7,9,11)为铰元件的传递矩阵，z1,0为工业机器人系统固定端边界点状态矢量。
5.根据权利要求1所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S5具体为：
根据各体元件传递方程(17)和各铰元件传递方程(29)，得到机器人系统总传递方程
z12,0＝U1-12z1,0 (30)
式中，U1-12＝U12…U2U1为系统总传递矩阵；
工业机器人固定端边界点和自由端边界点的边界条件为
Figure FDA0003249931780000061
将边界条件(31)代入式(30)，得到分块矩阵形式为
Figure FDA0003249931780000062
Figure FDA0003249931780000063
优化式(33)，计算边界点状态矢量中的未知状态变量
Figure FDA0003249931780000064
根据计算得到的未知状态变量获得系统固定端状态矢量，然后依次利用各元件的传递方程(17)和(29)，获得系统中各连接点的状态矢量。
6.根据权利要求1所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S6中铰元件θj(j＝1,3,5,7,9,11)的广义坐标的二阶导数的计算过程为：在式(18)两边同时左乘
Figure FDA0003249931780000065
得到
Figure FDA0003249931780000066
CN202010875365.4A 2020-08-27 2020-08-27 一种旋转关节型工业机器人非线性动力学建模分析方法 Active CN112001087B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010875365.4A CN112001087B (zh) 2020-08-27 2020-08-27 一种旋转关节型工业机器人非线性动力学建模分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010875365.4A CN112001087B (zh) 2020-08-27 2020-08-27 一种旋转关节型工业机器人非线性动力学建模分析方法

Publications (2)

Publication Number Publication Date
CN112001087A CN112001087A (zh) 2020-11-27
CN112001087B true CN112001087B (zh) 2021-10-19

Family

ID=73470459

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010875365.4A Active CN112001087B (zh) 2020-08-27 2020-08-27 一种旋转关节型工业机器人非线性动力学建模分析方法

Country Status (1)

Country Link
CN (1) CN112001087B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113954079B (zh) * 2021-11-23 2023-03-24 北京邮电大学 一种同质模块化机器人通用数学表征方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104793497A (zh) * 2015-04-20 2015-07-22 天津理工大学 基于多体系统离散时间传递矩阵法的机器人动力学建模方法
CN106708078A (zh) * 2017-02-21 2017-05-24 西北工业大学 一种适用于空间机器人执行器故障下的快速姿态稳定方法
CN109543264A (zh) * 2018-11-12 2019-03-29 天津理工大学 一种基于多维度重构校正的柔性多体机器人建模与求解方法
CN110340888A (zh) * 2018-10-30 2019-10-18 大连理工大学 一种空间机器人抓捕控制系统、强化学习方法及动力学建模方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104793497A (zh) * 2015-04-20 2015-07-22 天津理工大学 基于多体系统离散时间传递矩阵法的机器人动力学建模方法
CN106708078A (zh) * 2017-02-21 2017-05-24 西北工业大学 一种适用于空间机器人执行器故障下的快速姿态稳定方法
CN110340888A (zh) * 2018-10-30 2019-10-18 大连理工大学 一种空间机器人抓捕控制系统、强化学习方法及动力学建模方法
CN109543264A (zh) * 2018-11-12 2019-03-29 天津理工大学 一种基于多维度重构校正的柔性多体机器人建模与求解方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
DEVELOPMENTS IN TRANSFER MATRIX METHOD FOR MULTIBODY SYSTEMS (RUI METHOD) AND ITS APPLICATIONS;Xiaoting Rui;《IDETC/CIE 2018》;20190404;全文 *
Transfer Matrix Method for Multi-body Systems Dynamics Modeling of Rescue Robot;Yue Guo, Wei Chen and Xinhua Zhao;《IEEE》;20140828;全文 *
基于多刚体系统传递矩阵法的救援机器人动力学建模与规划;郭月;《中国优秀博硕士学位论文全文数据库(硕士).信息科技辑》;20151215;对比文件1正文 *

Also Published As

Publication number Publication date
CN112001087A (zh) 2020-11-27

Similar Documents

Publication Publication Date Title
CN104723340B (zh) 基于连接和阻尼配置的柔性关节机械臂的阻抗控制方法
CN103495977B (zh) 一种6r型工业机器人负载识别方法
CN111993417B (zh) 一种基于rbf神经网络的机械臂自适应阻抗控制方法
Shao et al. Research on the inertia matching of the Stewart parallel manipulator
CN109657282B (zh) 一种基于拉格朗日动力学的h型运动平台建模方法
CN104723341A (zh) 基于连接和阻尼配置的柔性关节机械臂的位置控制方法
CN109634111B (zh) 一种高速重载机器人动态变形计算方法
CN102962838A (zh) 具有封闭式运动学正解的六自由度并联机构及解析方法
CN112001087B (zh) 一种旋转关节型工业机器人非线性动力学建模分析方法
Xie et al. Dynamic modeling and performance analysis of a new redundant parallel rehabilitation robot
CN110909438A (zh) 一种基于动力学模型的轻载关节型并联机器人控制方法
Shao et al. Driving force analysis for the secondary adjustable system in FAST
CN110774286B (zh) 一种基于刚柔耦合动力学的五自由度机械手的控制方法
Noppeney et al. Task-space impedance control of a parallel Delta robot using dual quaternions and a neural network
CN114888812A (zh) 一种基于刚度性能优化的移动机械臂站位规划方法
CN113821935B (zh) 基于对称约束的动力学模型的建立方法及其系统
Shan et al. A trajectory planning and simulation method for welding robot
Liu et al. N-PD cross-coupling synchronization control based on adjacent coupling error analysis
Yang et al. Research on gravity compensation in motion control of multi-joint robot
Wang et al. Newton-Euler method for dynamic modeling and control of parallel polishing manipulator
Ha et al. Wireless-communicated computed-torque control of a SCARA robot and two-dimensional input shaping for a spherical pendulum
Xu et al. Kinematic analysis of a novel type of 6-DOF hybrid serial-parallel manipulator
Hovland et al. Benchmark of the 3-dof gantry-tau parallel kinematic machine
CN109571483B (zh) 一种空间机械臂任务轨迹规划域构建方法
Liu et al. Structure Design and Dynamic Simulation Analysis of Dual-arm Handling Robot

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20220130

Address after: No. 29, Qinhuai District, Qinhuai District, Nanjing, Jiangsu

Patentee after: Nanjing University of Aeronautics and Astronautics Asset Management Co.,Ltd.

Address before: No. 29, Qinhuai District, Qinhuai District, Nanjing, Jiangsu

Patentee before: Nanjing University of Aeronautics and Astronautics

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220223

Address after: 211599 No. 59, Wangqiao Road, Xiongzhou street, Liuhe District, Nanjing, Jiangsu

Patentee after: Nanjing Changjiang Industrial Technology Research Institute Co.,Ltd.

Address before: No. 29, Qinhuai District, Qinhuai District, Nanjing, Jiangsu

Patentee before: Nanjing University of Aeronautics and Astronautics Asset Management Co.,Ltd.

TR01 Transfer of patent right