CN109634111A - 一种高速重载机器人动态变形计算方法 - Google Patents
一种高速重载机器人动态变形计算方法 Download PDFInfo
- Publication number
- CN109634111A CN109634111A CN201710928582.3A CN201710928582A CN109634111A CN 109634111 A CN109634111 A CN 109634111A CN 201710928582 A CN201710928582 A CN 201710928582A CN 109634111 A CN109634111 A CN 109634111A
- Authority
- CN
- China
- Prior art keywords
- coordinate system
- robot
- relative
- connecting rod
- joint
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 36
- 230000001133 acceleration Effects 0.000 claims abstract description 90
- 230000033001 locomotion Effects 0.000 claims abstract description 14
- 238000013461 design Methods 0.000 claims abstract description 8
- 239000013598 vector Substances 0.000 claims description 84
- 239000011159 matrix material Substances 0.000 claims description 64
- 238000006073 displacement reaction Methods 0.000 claims description 39
- 238000004364 calculation method Methods 0.000 claims description 35
- 230000008569 process Effects 0.000 claims description 17
- 230000014509 gene expression Effects 0.000 claims description 8
- 230000003068 static effect Effects 0.000 claims description 8
- 238000009795 derivation Methods 0.000 claims description 6
- 230000005540 biological transmission Effects 0.000 claims description 4
- 239000003638 chemical reducing agent Substances 0.000 claims description 4
- 230000005484 gravity Effects 0.000 claims description 4
- 230000010354 integration Effects 0.000 claims description 4
- 239000000463 material Substances 0.000 claims description 4
- 230000000704 physical effect Effects 0.000 claims description 4
- 238000005457 optimization Methods 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 5
- IAZDPXIOMUYVGZ-UHFFFAOYSA-N Dimethylsulphoxide Chemical compound CS(C)=O IAZDPXIOMUYVGZ-UHFFFAOYSA-N 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 238000003466 welding Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000005507 spraying Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
- G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
- G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
- G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Manipulator (AREA)
Abstract
本发明涉及一种高速重载机器人动态变形计算方法，根据机器人本体设计特征，在机器人基座处建立基坐标系，在机器人各关节处建立机器人关节坐标系，获取各关节处D‑H参数；计算末端装置的运动学方程，获取机器人的末端位姿，依次进行角速度与角加速度递推以及依次进行速度和加速度递推，计算各连杆在基坐标系下的运动状态；求解各连杆所受惯性力与惯性力矩；建立机器人有限元模型；求出机器人总的动态变形、应力及应变。本发明将机器人的动力学模型运用到有限元模型中，分析了由机器人自重、负载及各项惯性力与惯性力矩引起的变形情况，研究了末端点变形的各项影响因素，为机器人的造型设计、参数优化提供了依据。
Description
技术领域
本发明涉及高速重载机器人领域，具体地说是高速重载机器人动态变形计算方法。
背景技术
高速重载机器人被定义为其末端载荷为100kg以上，末端运行的线速度最大值在1m/s以上的工业机器人，广泛应用于汽车的喷涂、焊接，物流的搬运、堆垛以及重型制造业的自动化生产等行业。高速重载机器人因作业环境和任务的需要，具有高速重载、高刚度、高平稳性、体积重量大等特点。除了由末端负载和机器人自重引起的静变形外，机器人运行过程中的惯性力和惯性力矩也会引起机器人的变形，由此将会带来末端轨迹偏移、机器人振动、零件受损、精度下降等一系列问题。
为解决上述问题，很有必要建立对机器人动态变形的评估计算方法，以便为机器人的结构设计和运动控制提供依据。在实际工程应用中，常将计算问题简化，采用弹性有限元方法分析静变形，再通过类比法估计机器人处于运动状态下的变形[Zienkiewicz OC.The Finite Elements Methods Third Edition[M].New York：McGraw-Hill，inc.1977]。然而随着对工业机器人运动特性的要求不断提高，惯性力的影响已经不可忽视，现有的计算方法通常利用整机柔度矩阵来计算机器人的末端变形[D.A.Fresonke,E.Hernandez,and D.Tesar.Deflection Predictions for Serial Manipulators[C]//IEEE Conference on Robotics and Automation.Philadelphia,PA:IEEE,1993:pp482-487]，但这种方法没有考虑机器人静态时连杆变形对柔性矩阵的影响，建立的柔度矩阵通常也只能计算末端点变形。在整机动态变形计算方面，一直没有对各关节运动耦合对末端造成的影响作出解释，也没有分析机器人运动过程中的每一项力对总变形造成的影响。如果计算方法与参数设计不当，很容易对机器人的使用产生不良影响，反之，则可以对机器人的设计、制造与使用起到参考意义。
发明内容
针对现有技术的不足，本发明提供一种高速重载机器人动态变形计算方法，目的在于针对高速重载机器人在工作过程中的动态变形大大超出静变形，容易造成机器人振动幅度过大、末端轨迹偏移等问题，。
本发明为实现上述目的所采用的技术方案是：
一种高速重载机器人动态变形计算方法，包括以下步骤：
步骤1：根据机器人本体设计特征，在机器人基座处建立基坐标系{0}，在机器人各关节处建立机器人关节坐标系，获取各关节处D-H参数；
步骤2：计算末端装置的运动学方程，获取机器人的末端位姿，依次进行角速度与角加速度递推以及依次进行速度和加速度递推，计算各连杆在基坐标系下的运动状态；
步骤3：根据各连杆在基坐标系下的运动状态，求解各连杆所受惯性力与惯性力矩；
步骤4：根据机器人的三维模型建立机器人有限元模型；
步骤5：分别计算由机器人自重、负载及各项惯性力与惯性力矩引起的变形情况，将每项载荷引起的变形线性叠加，求出机器人总的动态变形、应力及应变。
所述依次进行角速度与角加速度递推，包括以下步骤：
递推从连杆1开始，在基座标系{0}下，0ω0， 0v0与均为0，作为递推初值；
建立连杆i在基坐标系下的角速度与角加速度递推公式。
所述建立连杆i在基坐标系下的角速度与角加速度递推公式包括以下过程：
步骤1：设参考系为{A}，坐标系{B}相对于坐标系{A}的角速度为AωB，坐标系{C}相对于坐标系{B}的角速度为CωB，若将BωC也转换为参考系{A}中的描述，则坐标系{C}相对于坐标系{A}的角速度为两相对角速度之矢量和，即：
其中，AωC为坐标系{C}相对于坐标系{A}的角速度，为坐标系{B}相对于坐标系{A}的旋转矩阵；
步骤2：对公式(1)求导，得：
其中，为坐标系{C}相对于坐标系{A}的角加速度，为坐标系{B}相对于坐标系{A}的角加速度，为坐标系{C}相对于坐标系{B}的角加速度；
步骤3：令A＝0，B＝i-1，C＝i，则连杆i在基坐标系{0}下的角速度与角加速度为：
其中，0ωi为连杆i在基坐标系{0}下的角速度，0ωi-1为连杆i-1在基坐标系{0}下的角速度，i-1ωi为连杆i在坐标系{i-1}下的角速度，为坐标系{i-1}相对于基坐标系{0}的旋转矩阵，为坐标系{i}相对于基坐标系{0}的角加速度，为坐标系{i-1}相对于基坐标系{0}的角加速度，为坐标系{i}相对于坐标系{i-1}的角加速度。
所述依次进行速度和加速度递推包括以下过程：
步骤1：导出连杆i上任意一点的速度计算公式；
步骤2：导出在基坐标系下，各连杆关节处与质心处的速度递推公式；
步骤3：导出在基坐标系下，各连杆关节处与质心处的加速度递推公式。
所述导出连杆i上任意一点的速度计算公式，包括：
将连杆i作为空间中的自由刚体，在关节坐标系{i}的原点O'处创建辅助坐标系{O'}，使该坐标系的原点始终与关节坐标系{i}的原点重合，各坐标轴的方向与基座标系{0}保持同向；
连杆i上任意一点M的速度为：
vM＝vo'+vr＝vo'+ωr×r' (5)
其中，vM为点M相对于基坐标系{0}的速度，vo’为点O'相对于基坐标系{0}的速度，vr为点M相对于辅助坐标系{O'}的速度，ωr为连杆i绕原点O'的瞬时角速度，r'为连杆i上任意点M相对于辅助坐标系{O'}的位置矢量。
所述各连杆关节处与质心处的速度递推公式为：
其中，ipi+1为坐标系{i+1}相对于坐标系{i}的位置矢量，0vi为连杆i关节处相对于基坐标系{0}的速度，0vi+1为连杆i+1关节处相对于基坐标系{0}的速度，iωi为连杆i相对于坐标系{i}的角速度，为坐标系{i}相对于基坐标系{0}的旋转矩阵，ipci为连杆i的质心相对于坐标系{i}的位置矢量，0vci为连杆i的质心相对于基坐标系{0}的速度。
所述各连杆关节处与质心处的加速度递推公式为：
其中，为连杆i+1关节处相对于基坐标系{0}的加速度，为连杆i相对于基坐标系{0}的加速度，为坐标系{i}相对于基坐标系{0}的旋转矩阵，为连杆i相对于关节坐标系{i}的角加速度，ipi+1为坐标系{i+1}相对于坐标系{i}的位置矢量，iωi为连杆i相对于坐标系{i}的角速度，为连杆i的质心相对于基坐标系{0}的加速度，为连杆i关节处相对于基坐标系{0}的加速度，ipci为坐标系{i+1}相对于坐标系{i}的位置矢量。
所述各连杆所受惯性力为：
其中，为连杆i质心所受惯性力，mi为连杆i的质量，为连杆i的质心相对于基坐标系{0}的加速度；
各连杆所受惯性力矩为：
其中，Mi为连杆i所受惯性力矩，和分别为由于连杆i的角加速度和角速度产生的惯性力，ic为连杆i在坐标系{C}中的惯性张量，坐标系{C}的原点位于连杆i质心，各坐标轴方位与原连杆坐标系{i}方位相同，为连杆i相对于基坐标系{0}的角加速度，0ωi为连杆i相对于基坐标系{0}的角速度。
所述建立机器人有限元模型，包括以下步骤：
步骤1：建立机器人三维模型，并对模型进行简化，将减速器、齿轮、传动轴、轴承等与相应杆件固连，将没有相对运动的零件连接重组为一个构件，去掉构件的倒角、圆角和一些沟槽、台阶；
步骤2：确定各部分材料物理属性，包括弹性模量、泊松比、剪切模量、抗剪模量、质量密度、屈服强度；
步骤3：添加约束条件：在机器人各转动关节添加铰链约束，各移动关节添加移动副，机器人基座与地面固定，并施加引力场；
步骤4：采用四面体结构对机器人模型划分单元网格，计算各单元刚度矩阵：
Ke＝∫∫∫BTDBdxdydz＝BTDBV (12)
其中，B为几何矩阵，D为弹性矩阵，V为单元体的体积，Ke为单元刚度矩阵；
步骤5：采用刚度集成法计算整体刚度矩阵K；
步骤6：根据弹性有限元法建立有限元方程，得出节点位移向量与节点载荷向量的关系表达式，有限元方程为：
FL＝K·{δ}e (13)
其中，FL为节点载荷向量，K为整体刚度矩阵，{δ}e＝{δ1,δ2……δn}T为单元节点位移向量；
步骤7：分别建立几何方程和物理方程，得出单元内应变向量、单元内应力向量同节点位移向量的表达式；
几何方程为：
ε＝B·{δ}e (14)
其中，ε＝{εx,εy,εz}T为单元内应变向量，B为几何矩阵，{δ}e为单元节点位移向量；
物理方程为：
σ＝Dε＝DB{δ}e＝S{δ}e (15)
其中，σ为单元内应力向量，ε为单元内应变向量，D为弹性矩阵，B为几何矩阵，S为应力应变矩阵，{δ}e为单元节点位移向量。
对于含n个连杆的串联机器人，所述机器人总的动态变形计算过程为：
Fi＝K·δi (16)
FS＝K·δ3n+1 (19)
其中，为连杆i质心所受惯性力，和分别表示由于连杆i的角加速度和角速度产生的惯性力矩，FS为机器人所受静力负载，包括重力、末端负载和弹簧力，δi、δn+i、δ2n+i、δ3n+1分别为由FS引起的单元节点位移，i＝1，2……n，δ为机器人总的动态变形；
将连杆i质心所受惯性力代入公式(16)，解线性方程组，可求得δi；同理，求得δn+i、δ2n+i、δ3n+1，将所有位移按公式(20)叠加，即求得机器人总的动态变形。
机器人所受应变计算过程为：
由公式(14)求得机器人各单元节点位移所对应的单元内应变向量，然后按下式进行叠加：
ε＝∑εi (21)
其中，εi为各单元节点位移所对应的单元内应变向量，ε为机器人总的应变。
机器人所受应力计算过程为：
由公式(15)求得机器人各单元节点位移所对应的单元内应力向量，然后按下式进行叠加：
σ＝∑σi (22)
其中，σi为各单元节点位移所对应的单元内应变向量，σ为机器人总的应变。
本发明具有以下有益效果及优点：
1.本发明根据各关节的运动耦合关系，推导出了各连杆在基坐标系下的运动状态、所受惯性力与惯性力矩，并分别计算各力产生的动态变形，弥补了传统分析方法的不足；
2.本发明采用有限元法，对机器人自重、负载及各连杆惯性力与惯性力矩对动态变形的影响分别进行计算，便于逐个分析各项力对机器人动态变形的影响大小，为机器人的结构设计与运动控制提供了有效依据；
3.本发明可推广到其它多刚体串联系统，为多刚体串联系统的设计提供了有效的分析理论和仿真手段；
4.本发明已针对200kg点焊机器人进行了仿真实验，试验结果表明，依据本发明所做的仿真实验结果与现场试验结果误差小于5％，证明了计算方法的有效性。
附图说明
图1是本发明的方法流程图；
图2是关节坐标系建立图；
图3是连杆速度递推图；
图4是连杆i受力分析图。
具体实施方式
下面结合附图及实施例对本发明做进一步的详细说明。
如图1所示为本发明的方法流程图。
一种高速重载机器人动态变形计算方法，包括以下步骤：
(1)机器人运动学建模：根据机器人本体设计特征，建立机器人关节坐标系，获取机器人的运动学方程；
(2)机器人动力学建模：根据各连杆在基坐标系下的运动状态，求解各连杆所受惯性力与惯性力矩；
(3)建立机器人有限元模型：根据机器人的三维模型建立其有限元模型，分别求出其几何矩阵、应力应变矩阵和整体刚度矩阵，列出有限元方程、几何方程和物理方程。
(4)求解机器人总的动态变形：分别分析由机器人自重、负载及各项惯性力与惯性力矩引起的变形情况，将每项载荷引起的变形线性叠加，求出机器人总的动态变形。
所述机器人运动学建模包括以下步骤：
(1.1)根据机器人本体的结构与尺寸，在机器人基座处建立基坐标系{0}，在机器人各关节处建立关节坐标系，获取各关节处D-H参数。
(1.2)计算机器人各关节坐标系相对于基坐标系{0}的变换矩阵，获取机器人的末端位姿。
(1.3)依次进行角速度与角加速度递推：
递推从连杆1开始，在基座标系{0}下，0ω0， 0v0与均为0，作为递推初值。
设参考系为{A}，坐标系{B}相对于坐标系{A}的角速度为AωB，坐标系{C}相对于坐标系{B}的角速度为CωB，若将BωC也转换为参考系{A}中的描述，则坐标系{C}相对于坐标系{A}的角速度为两相对角速度之矢量和，即：
其中，AωC为坐标系{C}相对于坐标系{A}的角速度，为坐标系{B}相对于坐标系{A}的旋转矩阵。
对公式(1)求导，得：
其中，为坐标系{C}相对于坐标系{A}的角加速度，为坐标系{B}相对于坐标系{A}的角加速度，为坐标系{C}相对于坐标系{B}的角加速度。
令A＝0，B＝i-1，C＝i，则连杆i在基坐标系{0}下的角速度与角加速度为：
其中，0ωi为连杆i在基坐标系{0}下的角速度，0ωi-1为连杆i-1在基坐标系{0}下的角速度，i-1ωi为连杆i在坐标系{i-1}下的角速度，为坐标系{i-1}相对于基坐标系{0}的旋转矩阵，为坐标系{i}相对于基坐标系{0}的角加速度，为坐标系{i-1}相对于基坐标系{0}的角加速度，为坐标系{i}相对于坐标系{i-1}的角加速度。
(1.4)依次进行速度与加速度递推：
将连杆i作为空间中的自由刚体，在关节坐标系{i}的原点O'处创建辅助坐标系{O'}，使该坐标系的原点始终与关节坐标系{i}的原点重合，各坐标轴的方向与基座标系{0}保持同向；
连杆i上任意一点M的速度为：
vM＝vo'+vr＝vo'+ωr×r' (5)
其中，vM为点M相对于基坐标系{0}的速度，vo’为点O'相对于基坐标系的速度，vr为点M相对于辅助坐标系{O'}的速度，ωr为连杆i绕原点O'的瞬时角速度，r'为连杆i上任意点M相对于辅助坐标系{O'}的位置矢量。
按上式可以写出在基坐标系{0}下，各连杆关节处与质心处的速度递推公式：
其中，ipi+1为坐标系{i+1}相对于坐标系{i}的位置矢量，0vi为连杆i关节处相对于基坐标系{0}的速度，0vi+1为连杆i+1关节处相对于基坐标系{0}的速度，iωi为连杆i相对于坐标系{i}的角速度，为坐标系{i}相对于基坐标系{0}的旋转矩阵，ipci为连杆i的质心相对于坐标系{i}的位置矢量，0vci为连杆i的质心相对于基坐标系{0}的速度。
分别对公式(6)(7)求导，得到各连杆关节处与质心处的加速度递推公式为：
其中，为连杆i+1关节处相对于基坐标系{0}的加速度，为连杆i相对于基坐标系{0}的加速度，为坐标系{i}相对于基坐标系{0}的旋转矩阵，为连杆i相对于关节坐标系{i}的角加速度，ipi+1为坐标系{i+1}相对于坐标系{i}的位置矢量，iωi为连杆i相对于坐标系{i}的角速度，为连杆i的质心相对于基坐标系{0}的加速度，为连杆i关节处相对于基坐标系{0}的加速度，ipci为坐标系{i+1}相对于坐标系{i}的位置矢量。
所述机器人动力学建模包括以下步骤：
(2.1)连杆i质心所受惯性力为：
其中，为连杆i质心所受惯性力，mi为连杆i的质量，为连杆i的质心相对于基坐标系{0}的加速度；
各连杆所受惯性力矩为：
其中，Mi为连杆i所受惯性力矩，和分别为由于连杆i的角加速度和角速度产生的惯性力，ic为连杆i在坐标系{C}中的惯性张量，坐标系{C}的原点位于连杆i质心，各坐标轴方位与原连杆坐标系{i}方位相同，为连杆i相对于基坐标系{0}的角加速度，0ωi为连杆i相对于基坐标系{0}的角速度。
所述建立机器人有限元模型，包括以下步骤：
(3.1)建立机器人三维模型，并对模型进行简化，将减速器、齿轮、传动轴、轴承等与相应杆件固连，将没有相对运动的零件连接重组为一个构件，去掉构件的倒角、圆角和一些沟槽、台阶。
(3.2)确定各部分材料物理属性，包括弹性模量、泊松比、剪切模量、抗剪模量、质量密度、屈服强度；
(3.3)添加约束条件：在机器人各转动关节添加铰链约束，各移动关节添加移动副，机器人基座与地面固定，并施加引力场；
(3.4)采用四面体结构对机器人模型划分单元网格，计算各单元刚度矩阵：
Ke＝∫∫∫BTDBdxdydz＝BTDBV (12)
其中，B为几何矩阵，D为弹性矩阵，V为单元体的体积，Ke为单元刚度矩阵；
(3.5)采用刚度集成法计算整体刚度矩阵K。
(3.6)根据弹性有限元法建立有限元方程，得出节点位移向量与节点载荷向量的关系表达式，有限元方程为：
FL＝K·{δ}e (13)
其中，FL为节点载荷向量，K为整体刚度矩阵，{δ}e＝{δ1,δ2……δn}T为节点位移向量；
(3.7)分别建立几何方程和物理方程，得出单元内应变向量、单元内应力向量与节点位移向量的关系表达式。
几何方程为：
ε＝B·{δ}e (14)
其中，ε＝{εx,εy,εz}T为单元内应变向量，B为几何矩阵，{δ}e为节点位移向量；
物理方程为：
σ＝Dε＝DB{δ}e＝S{δ}e (15)
其中，σ为单元内应力向量，ε为单元内应变向量，D为弹性矩阵，B为几何矩阵，S为应力应变矩阵，{δ}e为节点位移向量。
所述求解机器人总的动态变形包括以下步骤：
(4.1)对于含n个连杆的串联机器人，所述机器人总的动态变形计算过程为：
Fi＝K·δi (16)
FS＝K·δ3n+1 (19)
其中，为连杆i质心所受惯性力，和分别表示由于连杆i的角加速度和角速度产生的惯性力矩，FS为机器人所受静力负载，包括重力、末端负载和弹簧力，δi、δn+i、δ2n+i、δ3n+1分别为由FS引起的位移，i＝1，2……n，δ为机器人总的动态变形。
将连杆i质心所受惯性力代入公式(16)，解线性方程组，可求得δi。同理，求得δn+i、δ2n+i、δ3n+1。将所有位移按公式(20)叠加，即求得机器人总的动态变形。
(4.2)机器人所受应变计算过程为：
由公式(14)求得机器人各单元节点位移所对应的单元内应变向量，然后按下式进行叠加：
ε＝∑εi (21)
其中，εi为各单元节点位移所对应的单元内应变向量，ε为机器人总的应变。
(4.3)机器人所受应力计算过程为：
由公式(15)求得机器人各单元节点位移所对应的单元内应力向量，然后按下式进行叠加：
σ＝∑σi (22)
其中，σi为各单元节点位移所对应的单元内应变向量，σ为机器人总的应变。
实施例1：
具体流程如下：
(1)机器人运动学建模步骤
(1.1)如图2所示为本发明的关节坐标系建立图。
根据机器人本体的结构与尺寸，在机器人基座处建立基坐标系{0}，在机器人各关节处建立关节坐标系，获取各关节处D-H参数。
(1.2)计算机器人各关节坐标系相对于基坐标系{0}的变换矩阵，获取机器人的末端位姿。
(1.3)依次进行角速度与角加速度递推：
递推从连杆1开始，在基座标系{0}下，0ω0， 0v0与均为0，作为递推初值。
设参考系为{A}，坐标系{B}相对于坐标系{A}的角速度为AωB，坐标系{C}相对于坐标系{B}的角速度为CωB，若将BωC也转换为参考系{A}中的描述，则坐标系{C}相对于坐标系{A}的角速度为两相对角速度之矢量和，即：
其中，AωC为坐标系{C}相对于坐标系{A}的角速度，为坐标系{B}相对于坐标系{A}的旋转矩阵。
对公式(1)求导，得：
其中，为坐标系{C}相对于坐标系{A}的角加速度，为坐标系{B}相对于坐标系{A}的角加速度，为坐标系{C}相对于坐标系{B}的角加速度。
令A＝0，B＝i-1，C＝i，则连杆i在基坐标系{0}下的角速度与角加速度为：
其中，0ωi为连杆i在基坐标系{0}下的角速度，0ωi-1为连杆i-1在基坐标系{0}下的角速度，i-1ωi为连杆i在坐标系{i-1}下的角速度，为坐标系{i-1}相对于基坐标系{0}的旋转矩阵，为坐标系{i}相对于基坐标系{0}的角加速度，为坐标系{i-1}相对于基坐标系{0}的角加速度，为坐标系{i}相对于坐标系{i-1}的角加速度。
(1.4)依次进行速度与加速度递推：
如图3所示为本发明的连杆速度递推图。
将连杆i作为空间中的自由刚体，在关节坐标系{i}的原点O'处创建辅助坐标系{O'}，使该坐标系的原点始终与关节坐标系{i}的原点重合，各坐标轴的方向与基座标系{0}保持同向；
连杆i上任意一点M的速度为：
vM＝vo'+vr＝vo'+ωr×r' (5)
其中，vM为点M相对于基坐标系{0}的速度，vo’为点O'相对于基坐标系的速度，vr为点M相对于辅助坐标系{O'}的速度，ωr为连杆i绕原点O'的瞬时角速度，r'为连杆i上任意点M相对于辅助坐标系{O'}的位置矢量。
按上式可以写出在基坐标系{0}下，各连杆关节处与质心处的速度递推公式：
其中，ipi+1为坐标系{i+1}相对于坐标系{i}的位置矢量，0vi为连杆i关节处相对于基坐标系{0}的速度，0vi+1为连杆i+1关节处相对于基坐标系{0}的速度，iωi为连杆i相对于坐标系{i}的角速度，为坐标系{i}相对于基坐标系{0}的旋转矩阵，ipci为连杆i的质心相对于坐标系{i}的位置矢量，0vci为连杆i的质心相对于基坐标系{0}的速度。
分别对公式(6)、(7)求导，得到各连杆关节处与质心处的加速度递推公式为：
其中，为连杆i+1关节处相对于基坐标系{0}的加速度，为连杆i相对于基坐标系{0}的加速度，为坐标系{i}相对于基坐标系{0}的旋转矩阵，为连杆i相对于关节坐标系{i}的角加速度，ipi+1为坐标系{i+1}相对于坐标系{i}的位置矢量，iωi为连杆i相对于坐标系{i}的角速度，为连杆i的质心相对于基坐标系{0}的加速度，为连杆i关节处相对于基坐标系{0}的加速度，ipci为坐标系{i+1}相对于坐标系{i}的位置矢量。
(2)如图4所示为本发明的连杆i受力分析图。
机器人动力学建模步骤
(2.1)连杆i质心所受惯性力为：
其中，为连杆i质心所受惯性力，mi为连杆i的质量，为连杆i的质心相对于基坐标系{0}的加速度；
(2.2)各连杆所受惯性力矩为：
其中，Mi为连杆i所受惯性力矩，和分别为由于连杆i的角加速度和角速度产生的惯性力，ic为连杆i在坐标系{C}中的惯性张量，坐标系{C}的原点位于连杆i质心，各坐标轴方位与原连杆坐标系{i}方位相同，为连杆i相对于基坐标系{0}的角加速度，0ωi为连杆i相对于基坐标系{0}的角速度。
(3)建立机器人有限元模型步骤
(3.1)建立机器人三维模型，并对模型进行简化，将减速器、齿轮、传动轴、轴承等与相应杆件固连，将没有相对运动的零件连接重组为一个构件，去掉构件的倒角、圆角和一些沟槽、台阶；
(3.2)确定各部分材料物理属性，包括弹性模量、泊松比、剪切模量、抗剪模量、质量密度、屈服强度；
(3.3)添加约束条件：在机器人各转动关节添加铰链约束，各移动关节添加移动副，机器人基座与地面固定，并施加引力场；
(3.4)采用四面体结构对机器人模型划分单元网格，计算各单元刚度矩阵：
Ke＝∫∫∫BTDBdxdydz＝BTDBV (12)
其中，B为几何矩阵，D为弹性矩阵，V为单元体的体积，Ke为单元刚度矩阵；
(3.5)采用刚度集成法计算整体刚度矩阵K。
(3.6)根据弹性有限元法建立有限元方程，得出节点位移向量与节点载荷向量的关系表达式，有限元方程为：
FL＝K·{δ}e (13)
其中，FL为节点载荷向量，K为整体刚度矩阵，{δ}e＝{δ1,δ2……δn}T为节点位移向量；
(3.7)分别建立几何方程和物理方程，得出单元内应变向量、单元内应力向量与节点位移向量的关系表达式。
几何方程为：
ε＝B·{δ}e (14)
其中，ε＝{εx,εy,εz}T为单元内应变向量，B为几何矩阵，{δ}e为节点位移向量；
物理方程为：
σ＝Dε＝DB{δ}e＝S{δ}e (15)
其中，σ为单元内应力向量，ε为单元内应变向量，D为弹性矩阵，B为几何矩阵，S为应力应变矩阵，{δ}e为节点位移向量。
(4)求解机器人总的动态变形步骤
对于含n个连杆的串联机器人，所述机器人总的动态变形计算过程为：
Fi＝K·δi (16)
FS＝K·δ3n+1 (19)
其中，为连杆i质心所受惯性力，和分别表示由于连杆i的角加速度和角速度产生的惯性力矩，FS为机器人所受静力负载，包括重力、末端负载和弹簧力，δi、δn+i、δ2n+i、δ3n+1分别为由FS引起的位移，i＝1，2……n，δ为机器人总的动态变形。
将连杆i质心所受惯性力代入公式(16)，解线性方程组，可求得δi。同理，求得δn+i、δ2n+i、δ3n+1。将所有位移按公式(20)叠加，即求得机器人总的动态变形。
由公式(14)求得机器人各单元节点位移所对应的单元内应变向量，然后按下式进行叠加：
ε＝∑εi (21)
其中，εi为各单元节点位移所对应的单元内应变向量，ε为机器人总的应变。
机器人所受应力计算过程为：
由公式(15)求得机器人各单元节点位移所对应的单元内应力向量，然后按下式进行叠加：
σ＝∑σi (22)
其中，σi为各单元节点位移所对应的单元内应变向量，σ为机器人总的应变。
Claims (10)
1.一种高速重载机器人动态变形计算方法，其特征在于：包括以下步骤：
步骤1：根据机器人本体设计特征，在机器人基座处建立基坐标系{0}，在机器人各关节处建立机器人关节坐标系，获取各关节处D-H参数；
步骤2：计算末端装置的运动学方程，获取机器人的末端位姿，依次进行角速度与角加速度递推以及依次进行速度和加速度递推，计算各连杆在基坐标系下的运动状态；
步骤3：根据各连杆在基坐标系下的运动状态，求解各连杆所受惯性力与惯性力矩；
步骤4：根据机器人的三维模型建立机器人有限元模型；
步骤5：分别计算由机器人自重、负载及各项惯性力与惯性力矩引起的变形情况，将每项载荷引起的变形线性叠加，求出机器人总的动态变形、应力及应变。
2.根据权利要求1所述的高速重载机器人动态变形计算方法，其特征在于：所述依次进行角速度与角加速度递推，包括以下步骤：
递推从连杆1开始，在基座标系{0}下，0ω0， 0v0与均为0，作为递推初值；
建立连杆i在基坐标系下的角速度与角加速度递推公式。
3.根据权利要求2所述的高速重载机器人动态变形计算方法，其特征在于：所述建立连杆i在基坐标系下的角速度与角加速度递推公式包括以下过程：
步骤1：设参考系为{A}，坐标系{B}相对于坐标系{A}的角速度为AωB，坐标系{C}相对于坐标系{B}的角速度为CωB，若将BωC也转换为参考系{A}中的描述，则坐标系{C}相对于坐标系{A}的角速度为两相对角速度之矢量和，即：
其中，AωC为坐标系{C}相对于坐标系{A}的角速度，为坐标系{B}相对于坐标系{A}的旋转矩阵；
步骤2：对公式(1)求导，得：
其中，为坐标系{C}相对于坐标系{A}的角加速度，为坐标系{B}相对于坐标系{A}的角加速度，为坐标系{C}相对于坐标系{B}的角加速度；
步骤3：令A＝0，B＝i-1，C＝i，则连杆i在基坐标系{0}下的角速度与角加速度为：
其中，0ωi为连杆i在基坐标系{0}下的角速度，0ωi-1为连杆i-1在基坐标系{0}下的角速度，i-1ωi为连杆i在坐标系{i-1}下的角速度，为坐标系{i-1}相对于基坐标系{0}的旋转矩阵，为坐标系{i}相对于基坐标系{0}的角加速度，为坐标系{i-1}相对于基坐标系{0}的角加速度，为坐标系{i}相对于坐标系{i-1}的角加速度。
4.根据权利要求1所述的高速重载机器人动态变形计算方法，其特征在于：所述依次进行速度和加速度递推包括以下过程：
步骤1：导出连杆i上任意一点的速度计算公式；
步骤2：导出在基坐标系下，各连杆关节处与质心处的速度递推公式；
步骤3：导出在基坐标系下，各连杆关节处与质心处的加速度递推公式。
5.根据权利要求4所述的高速重载机器人动态变形计算方法，其特征在于：所述导出连杆i上任意一点的速度计算公式，包括：
将连杆i作为空间中的自由刚体，在关节坐标系{i}的原点O'处创建辅助坐标系{O'}，使该坐标系的原点始终与关节坐标系{i}的原点重合，各坐标轴的方向与基座标系{0}保持同向；
连杆i上任意一点M的速度为：
vM＝vo'+vr＝vo'+ωr×r' (5)
其中，vM为点M相对于基坐标系{0}的速度，vo’为点O'相对于基坐标系{0}的速度，vr为点M相对于辅助坐标系{O'}的速度，ωr为连杆i绕原点O'的瞬时角速度，r'为连杆i上任意点M相对于辅助坐标系{O'}的位置矢量。
6.根据权利要求4所述的高速重载机器人动态变形计算方法，其特征在于：所述各连杆关节处与质心处的速度递推公式为：
其中，ipi+1为坐标系{i+1}相对于坐标系{i}的位置矢量，0vi为连杆i关节处相对于基坐标系{0}的速度，0vi+1为连杆i+1关节处相对于基坐标系{0}的速度，iωi为连杆i相对于坐标系{i}的角速度，为坐标系{i}相对于基坐标系{0}的旋转矩阵，ipci为连杆i的质心相对于坐标系{i}的位置矢量，0vci为连杆i的质心相对于基坐标系{0}的速度。
7.根据权利要求4所述的高速重载机器人动态变形计算方法，其特征在于：所述各连杆关节处与质心处的加速度递推公式为：
其中，为连杆i+1关节处相对于基坐标系{0}的加速度，为连杆i相对于基坐标系{0}的加速度，为坐标系{i}相对于基坐标系{0}的旋转矩阵，为连杆i相对于关节坐标系{i}的角加速度，ipi+1为坐标系{i+1}相对于坐标系{i}的位置矢量，iωi为连杆i相对于坐标系{i}的角速度，为连杆i的质心相对于基坐标系{0}的加速度，为连杆i关节处相对于基坐标系{0}的加速度，ipci为坐标系{i+1}相对于坐标系{i}的位置矢量。
8.根据权利要求1所述的高速重载机器人动态变形计算方法，其特征在于：
所述各连杆所受惯性力为：
其中，Fi i为连杆i质心所受惯性力，mi为连杆i的质量，为连杆i的质心相对于基坐标系{0}的加速度；
各连杆所受惯性力矩为：
其中，Mi为连杆i所受惯性力矩，和分别为由于连杆i的角加速度和角速度产生的惯性力，ic为连杆i在坐标系{C}中的惯性张量，坐标系{C}的原点位于连杆i质心，各坐标轴方位与原连杆坐标系{i}方位相同，为连杆i相对于基坐标系{0}的角加速度，0ωi为连杆i相对于基坐标系{0}的角速度。
9.根据权利要求1所述的高速重载机器人动态变形计算方法，其特征在于：所述建立机器人有限元模型，包括以下步骤：
步骤1：建立机器人三维模型，并对模型进行简化，将减速器、齿轮、传动轴、轴承等与相应杆件固连，将没有相对运动的零件连接重组为一个构件，去掉构件的倒角、圆角和一些沟槽、台阶；
步骤2：确定各部分材料物理属性，包括弹性模量、泊松比、剪切模量、抗剪模量、质量密度、屈服强度；
步骤3：添加约束条件：在机器人各转动关节添加铰链约束，各移动关节添加移动副，机器人基座与地面固定，并施加引力场；
步骤4：采用四面体结构对机器人模型划分单元网格，计算各单元刚度矩阵：
Ke＝∫∫∫BTDBdxdydz＝BTDBV (12)
其中，B为几何矩阵，D为弹性矩阵，V为单元体的体积，Ke为单元刚度矩阵；
步骤5：采用刚度集成法计算整体刚度矩阵K；
步骤6：根据弹性有限元法建立有限元方程，得出节点位移向量与节点载荷向量的关系表达式，有限元方程为：
FL＝K·{δ}e (13)
其中，FL为节点载荷向量，K为整体刚度矩阵，{δ}e＝{δ1,δ2……δn}T为单元节点位移向量；
步骤7：分别建立几何方程和物理方程，得出单元内应变向量、单元内应力向量同节点位移向量的表达式；
几何方程为：
ε＝B·{δ}e (14)
其中，ε＝{εx,εy,εz}T为单元内应变向量，B为几何矩阵，{δ}e为单元节点位移向量；
物理方程为：
σ＝Dε＝DB{δ}e＝S{δ}e (15)
其中，σ为单元内应力向量，ε为单元内应变向量，D为弹性矩阵，B为几何矩阵，S为应力应变矩阵，{δ}e为单元节点位移向量。
10.根据权利要求1所述的高速重载机器人动态变形计算方法，其特征在于：所述机器人总的动态变形计算过程为：
Fi＝K·δi (16)
FS＝K·δ3n+1 (19)
其中，Fi i为连杆i质心所受惯性力，和分别表示由于连杆i的角加速度和角速度产生的惯性力矩，FS为机器人所受静力负载，包括重力、末端负载和弹簧力，δi、δn+i、δ2n+i、δ3n+1分别为由Fi i、FS引起的单元节点位移，i＝1，2……n，δ为机器人总的动态变形；
所述机器人总的应变的计算过程为：
由公式(14)求得机器人各单元节点位移所对应的单元内应变向量，然后按下式进行叠加：
ε＝∑εi (21)
其中，εi为各单元节点位移所对应的单元内应变向量，ε为机器人总的应变；
所述机器人总的应力的计算过程为：
由公式(15)求得机器人各单元节点位移所对应的单元内应力向量，然后按下式进行叠加：
σ＝∑σi (22)
其中，σi为各单元节点位移所对应的单元内应力向量，σ为机器人总的应力。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710928582.3A CN109634111B (zh) | 2017-10-09 | 2017-10-09 | 一种高速重载机器人动态变形计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710928582.3A CN109634111B (zh) | 2017-10-09 | 2017-10-09 | 一种高速重载机器人动态变形计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109634111A true CN109634111A (zh) | 2019-04-16 |
| CN109634111B CN109634111B (zh) | 2020-08-04 |
Family
ID=66051193
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710928582.3A Active CN109634111B (zh) | 2017-10-09 | 2017-10-09 | 一种高速重载机器人动态变形计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109634111B (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110348140A (zh) * | 2019-07-15 | 2019-10-18 | 清华大学 | 基于拖曳距的双轮机器人建模与静止平衡方法及装置 |
| CN110414071A (zh) * | 2019-07-03 | 2019-11-05 | 上海航天设备制造总厂有限公司 | 一种重载机器人机构的参数同步优化方法 |
| CN111055285A (zh) * | 2020-01-08 | 2020-04-24 | 山东理工大学 | 一种仿人柔性关节手臂变负载工况下的振动抑制方法 |
| CN112818481A (zh) * | 2021-01-19 | 2021-05-18 | 中国科学院沈阳自动化研究所 | 弹性平面约束的细长软体机器人的建模与控制方法 |
| CN112959354A (zh) * | 2019-12-13 | 2021-06-15 | 中国科学院沈阳自动化研究所 | 一种引入弹性变形的机械臂标定方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2006040066A (ja) * | 2004-07-28 | 2006-02-09 | Nissan Motor Co Ltd | 歪み量解析装置およびその方法 |
| CN103970588A (zh) * | 2014-05-08 | 2014-08-06 | 中国人民解放军军事交通学院 | 一种Mecanum万向轮滚子轮廓补偿和修正建模仿真的方法 |
| CN107145662A (zh) * | 2017-05-04 | 2017-09-08 | 北京理工大学 | 一种介观尺度车削变形预测方法 |
| CN107169196A (zh) * | 2017-05-11 | 2017-09-15 | 西北工业大学 | 空间机器人由末端执行器向基座的动力学建模方法 |
-
2017
- 2017-10-09 CN CN201710928582.3A patent/CN109634111B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2006040066A (ja) * | 2004-07-28 | 2006-02-09 | Nissan Motor Co Ltd | 歪み量解析装置およびその方法 |
| CN103970588A (zh) * | 2014-05-08 | 2014-08-06 | 中国人民解放军军事交通学院 | 一种Mecanum万向轮滚子轮廓补偿和修正建模仿真的方法 |
| CN107145662A (zh) * | 2017-05-04 | 2017-09-08 | 北京理工大学 | 一种介观尺度车削变形预测方法 |
| CN107169196A (zh) * | 2017-05-11 | 2017-09-15 | 西北工业大学 | 空间机器人由末端执行器向基座的动力学建模方法 |
Non-Patent Citations (3)
| Title |
|---|
| 叶建辉: "SCARA机器人的运动分析及位姿误差研究", 《中国优秀硕士学位论文全文数据库·信息科技辑》 * |
| 徐稀文等: "弧焊机器人大臂有限元静态分析与优化设计", 《设计与工艺》 * |
| 杨国良: "工业机器人动力学仿真及有限元分析", 《中国优秀硕士学位论文全文数据库·信息科技辑》 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110414071A (zh) * | 2019-07-03 | 2019-11-05 | 上海航天设备制造总厂有限公司 | 一种重载机器人机构的参数同步优化方法 |
| CN110348140A (zh) * | 2019-07-15 | 2019-10-18 | 清华大学 | 基于拖曳距的双轮机器人建模与静止平衡方法及装置 |
| CN112959354A (zh) * | 2019-12-13 | 2021-06-15 | 中国科学院沈阳自动化研究所 | 一种引入弹性变形的机械臂标定方法 |
| CN112959354B (zh) * | 2019-12-13 | 2022-03-15 | 中国科学院沈阳自动化研究所 | 一种引入弹性变形的机械臂标定方法 |
| CN111055285A (zh) * | 2020-01-08 | 2020-04-24 | 山东理工大学 | 一种仿人柔性关节手臂变负载工况下的振动抑制方法 |
| CN111055285B (zh) * | 2020-01-08 | 2022-11-25 | 山东理工大学 | 一种仿人柔性关节手臂变负载工况下的振动抑制方法 |
| CN112818481A (zh) * | 2021-01-19 | 2021-05-18 | 中国科学院沈阳自动化研究所 | 弹性平面约束的细长软体机器人的建模与控制方法 |
| CN112818481B (zh) * | 2021-01-19 | 2023-11-07 | 中国科学院沈阳自动化研究所 | 弹性平面约束的细长软体机器人的建模与控制方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109634111B (zh) | 2020-08-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109634111B (zh) | 一种高速重载机器人动态变形计算方法 | |
| CN107703748B (zh) | 一种基于偏置板设计的重载机器人静态刚度辨识方法 | |
| CN110539302B (zh) | 一种工业机器人整体动力学建模及动力学参数辨识方法 | |
| CN109940613B (zh) | 一种计算含压电材料机械臂动力学响应及控制的仿真方法 | |
| CN108466289A (zh) | 一种考虑关节摩擦的并联机器人的动力学建模方法 | |
| CN108638070A (zh) | 基于动态平衡的机器人负载重力参数辨识方法 | |
| CN102495550A (zh) | 并联机器人的正、逆动力学响应分析与控制方法 | |
| Liu et al. | Dynamic modeling and analysis of 3-R RS parallel manipulator with flexible links | |
| Moldovan et al. | Dynamic numerical simulation of the 6-PGK parallel robot manipulator | |
| CN110549340A (zh) | 重载码垛机器人频响特性分析方法及系统 | |
| Fadaei et al. | Dynamics modeling of a stewart platform in Simulink MSC ADAMS | |
| Li et al. | Modeling and simulation research of Kane dynamics method for the 5-DOF modular industrial robot | |
| Zhang et al. | Modal and natural frequency sensitivity analysis of electrohydraulic Stewart platform | |
| Arora et al. | Finite element analysis and multibody dynamics of 6-dof industrial robot | |
| Qin et al. | Analyzing trajectory tracking accuracy of a flexible multi-purpose deployer | |
| Ye et al. | Motion optimization of an omnidirectional mobile robot with MY wheel based on contact mechanics | |
| CN114888812A (zh) | 一种基于刚度性能优化的移动机械臂站位规划方法 | |
| Luo et al. | Rigid-flexible coupling dynamics analysis of a spot-welding robot | |
| Wang et al. | Simulation and analysis of mechanical characteristics of a 6-dof spray-painting robot | |
| Liu et al. | Rigid-Flexible Coupling Simulation and Vibration Analysis of Flexible Robot | |
| Posiadala et al. | Modeling and analysis of the dynamics of load carrying system | |
| Zhang et al. | Analysis and Simulation of Mechanical Arm Dynamics Model Based on Simulink | |
| CN114434429B (zh) | 一种工业机器人动态精度设计方法和系统 | |
| Kaida et al. | Development of space-saving mechanism in shaking table systems for earthquake simulators | |
| Rueda et al. | Structural analysis of a delta-type parallel industrial robot using flexible dynamic of ANSYS 11.0 |
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 |