机械臂铣削加工运动分析和拓扑优化设计

Analysis of Motion and Topological Optimization Design of Industrial Robotic Arms

摘要: 工业机械臂在机械制造与产品生产中扮演着非常重要的角色，通常在打磨、铣削、抛光和搬运等领域均有运用，机械臂的运用必须建立在对其运动学知识的熟练掌握。本文以KUKA机械臂为研究对象，使用D-H参数法建立机械臂的数学模型，进行正逆向运动学的推导，并且考虑该机器人用于铣削加工时，使用Ansys进行有限元分析，发现连杆2设计较为保守，可以进行一定程度的拓扑优化，优化后连杆2质量降低28.8 kg。通过理论分析建模构建的KUKA机械臂的正逆向运动学模型能够使用关节角度计算机器人末端位置，或者已知末端位置求解关节角度，对六轴机械臂运动学建模有参考意义。

Abstract: Industrial robotic arms play a very important role in mechanical manufacturing and production. They are commonly used in grinding, milling, polishing, and handling, among other areas. The use of robotic arms must be based on a good understanding of kinematics. This paper takes the KUKA robot arm as the research object and uses the D-H parameter method to establish a mathematical model of the robot arm. The forward and inverse kinematics are derived, and finite element analysis is conducted using ANSYS when considering the robot arm used for milling operations. Link 2 is designed too conservatively, and topological optimization is possible. After optimization, the mass of link 2 is reduced by 28.8 kg. The forward and inverse kinematics model of the KUKA robot arm built through theoretical analysis and modeling can be used to calculate the robot end effector’s position or to solve the joint angles given the end effector position. It provides reference for the kinematics modeling of six-axis robot arms.

1. 引言

在工业制造中，机械臂的运用十分广泛，尤其是多关节机械臂运用更为普遍。工业机械臂通常可以由人类操作，或者是按照事先编写好的程序进行工作，在给人们带来非常巨大的经济效益的同时也为我们完成了很多难度较大、危险系数较高的工作，因此非常有必要对工业机械臂进行研究。《中国制造2025》中将机械臂作为需要重点攻坚克难的细分技术领域。由于国内的工业机械臂技术起步较晚，目前市场上使用的工业机械臂多是以进口为主，比如库卡、ABB等品牌。运动学分析是机械臂基础，机器臂的运动过程是基于运动学建立的。严大亮等[1]使用DH参数法研究了重载机械臂运动学的参数辨识问题。由于实际制造过程中存在一定误差，杨源[2]等使用标定板对机械臂的DH参数进行了重新标定。现成的机械臂通常具有较好的性能，但在处理特定问题通常具有一定的局限性，并且各自厂商在设计时对很多部件留有非常大的余量，因此针对特定运用条件下整体的安全性有较多的富余，一定程度上可以对相关部件进行优化设计，以减轻整体的质量，更加适应特定的工况，降低成本。通常在汽车表面喷涂、较轻物体搬运时对机械臂刚度要求不高，但是在焊接、铣削等领域运用时通常对精度、负载、刚度要求较高。天津理工大学徐赛华等[3]运用响应面法对重载机械臂进行了优化设计，一定程度上降低了整体的重量。宁波大学张刚[4]以MOTOMAN-HP20D串联型机器人为研究对象，进行了整机的动态分析与拓扑优化设计。国外学者Park J H等[5]研究了双连杆高速机器人机械结构优化设计与控制。拓扑优化一定程度上能使得机械臂更为轻便，以及结构更加紧凑[6]。

本文对目前市场上占有量比较大的KUKA机械臂进行了运动学分析，并针对机械臂连杆2进行了拓扑优化。

2. 运动学分析

2.1. KUKA机械臂运动学模型

运动学是以运动特性为研究目标，不考虑运动时产生的力[5]。机械臂的机械机构通常是由一系列刚性连杆通过关节连接起来，而关节的个数被称为机械臂的自由度。本文研究的KUKA机械臂运动关节有六个，即六个自由度。前三个关节用来确定机械臂末端位置，后三个关节用来确定姿态。图1是机械臂的模型图，R1~R6分别为旋转关节所在位置。

对于机械臂的运动学分析，Denavit和Hartenberg于1955年在ASME Journal of Applied Mechanics上发表的一篇论文提出了一种使用连杆参数描述机构运动关系的方法，此后该方法一直被广大学者使用，通常便被称作D-H参数法。按照D-H参数法求解机械臂运动学相关问题有以下5个步骤：

Figure 1. KUKA Robotic arm’s 3D model

图1. KUKA机械臂3D模型

1) 按照顺序找出关节轴。

2) 找出关节轴i与关节轴i+ 1之间的公垂线或者交点，以关节轴交点或公垂线与关节轴i的交点作为原点。

3) Z轴沿着关节轴的指向。

4) X轴沿着公垂线的指向。

5) 遵守右手定则确定Y轴。

根据图1机械臂的几何模型和D-H参数法建模的步骤可以抽象出机械臂的坐标系与相关参数，图2为机械臂连杆坐标系。表1为D-H参数法中涉及的参数。

Figure 2. KUKA robotic arm link coordinate system

图2. KUKA机械臂连杆坐标系

其中${\theta }_1$ 到${\theta }_6$ 为关节运动过程的旋转角度变量，正向运动学就是通过控制这6个变量从而控制末端的姿态，逆向运动学就是已知机器人末端位置求解${\theta }_1$ 到${\theta }_6$ 等六个关节变量。

Table 1. D-H parameters

表1. D-H参数

关节i	${\alpha }_{i-1}$	$a_{i-1}$ /mm	$d_i$ /mm	${\theta }_i$
1	0	0	0	${\theta }_1$
2	$-\pi /2$	260	0	${\theta }_2$
3	0	680	0	${\theta }_3$
4	$-\pi /2$	−35	670	${\theta }_4$
5	$\pi /2$	0	0	${\theta }_5$
6	$-\pi /2$	0	0	${\theta }_6$

2.2. 正向运动学分析

正向运动学是指每个关节的类型、关节变量均已知，需要通过坐标的变换最终得到末端执行器的在整体坐标系下的位置与姿态。DH参数法对应的相邻连杆i到i − 1之间的变换矩阵如下：

$\begin{array}{c}{}_i{}^{i-1}T=R_x\left({\alpha }_{i-1}\right)D_x\left(a_{i-1}\right)R_z\left({\theta }_i\right)D_z\left(d_i\right)\\ =\left[\begin{array}{cccc}\text{c}{\theta }_i& -\text{s}{\theta }_i& 0& {\alpha }_{i-1}\\ \text{s}{\theta }_i\text{c}{\alpha }_{i-1}& \text{c}{\theta }_i\text{c}{\alpha }_{i-1}& -\text{s}{\alpha }_{i-1}& -\text{s}{\alpha }_{i-1}{d}_i\\ \text{s}{\theta }_i\text{s}{\alpha }_{i-1}& \text{c}{\theta }_i\text{s}{\alpha }_{i-1}& \text{c}{\alpha }_{i-1}& \text{c}{\alpha }_{i-1}{d}_i\\ 0& 0& 0& 1\end{array}\right]\end{array}$ (1)

式中$c{\theta }_i=\cos {\theta }_i$ ，$s{\theta }_i=\sin {\theta }_i$ ，$\text{c}{\alpha }_{i-1}=\text{cos}{\alpha }_{i-1}$ ，$\text{s}{\alpha }_{i-1}=\text{sin}{\alpha }_{i-1}$ 。

将所得的单个连杆变换矩阵依次相乘便可以得到末端执行件相对于基座坐标系的位姿矩阵${}_6{}^{0}T$ 。

$\begin{array}{c}{}_6{}^{0}T={}_1{}^{0}T{}_2{}^{1}T{}_3{}^{2}T{}_4{}^{3}T{}_5{}^{4}T{}_6{}^{5}T=\left[\begin{array}{cccc}{n}_x& o_x& a_x& p_x\\ n_y& o_y& a_y& p_y\\ n_z& o_z& a_z& p_z\\ 0& 0& 0& 1\end{array}\right]\end{array}$ (2)

以上便是机械臂的正向运动学分析，在已知机械臂型号参数后便可通过机械臂关节角度计算机械臂末端参数。

2.3. 逆向运动学分析

逆向运动学顾名思义就是在已知机械臂末端位置和姿态，即已知公式2矩阵中的12个参数，反过来求解满足要求的关节角数值，即${\theta }_1$ 到${\theta }_6$ 的具体数值。逆向运动学与实际加工生产更加贴近，实际加工生产通常是知道机械臂的运动轨迹(通常可以理解为无数多个点)，需要控制各个关节的连续运动时的关节角，以保证机械臂末端在正确的位置。因此与正向运动学相比较，逆向运动学更为重要。可以使用代数方法求得各关节角度，即直接在变换矩阵两边依次乘与每个变换矩阵的逆矩阵，分离出关节变量，以下为逆向运动学求解过程。

首先分离变量求解关节角${\theta }_1$ 和${\theta }_3$ ，求解的过程就是矩阵乘法的逆运算。对式2中的变换矩阵${}_6{}^{0}T$ 左乘关节1的变换矩阵的逆矩阵${\left[{}_1^{0}T\left({\theta }_1\right)\right]}^{-1}$ ，等式两边元素对应相等建立方程便可求解${\theta }_1$ 和${\theta }_3$ 。

$\begin{array}{c}{\theta }_1=\arctan \frac{p_y}{p_x}\end{array}$ (3)

$\begin{array}{c}{\theta }_3=\arctan \frac{d_4}{a_3}-\arctan \left(\pm \frac{K}{\sqrt{{\rho }^2-K^2}}\right)\end{array}$ (4)

求解${\theta }_2$ 时需要选取关节1、2、3整体的变换矩阵的逆矩阵${\left[{}_3^{0}T\left(\theta \right)\right]}^{-1}$ 左乘式2，同样按照元素对应相等原理，并结合三角和差公式整理可得到${\theta }_2$ 。

$\begin{array}{c}{\theta }_{23}=\arctan \frac{\left(a_2{S}_3-d_4\right)\left(C_1{p}_x+S_1{p}_y-a_1\right)-p_z\left(a_3+a_2{C}_3\right)}{\left(a_3+a_2{C}_3\right)\left(C_1{p}_x+S_1{p}_y-a_1\right)-p_z\left(a_2{S}_3-d_4\right)}\end{array}$ (5)

$\begin{array}{c}{\theta }_2={\theta }_{23}-{\theta }_3\end{array}$ (6)

关节4、5、6的求解方式与前三个关节求解类似，均在式2的基础上左乘关节变换矩阵的逆矩阵，然后整理对应元素相等，即可求解。

$\begin{array}{c}{\theta }_4=\arctan \frac{S_1{a}_x-C_1{a}_y}{C_1{C}_{23}{a}_x+S_1{C}_{23}{a}_y-S_{23}{a}_z}\end{array}$ (7)

$\begin{array}{c}{\theta }_5=\arctan \frac{S_5}{C_5}\end{array}$ (8)

$\begin{array}{c}{\theta }_6=\arctan \frac{S_6}{C_6}\end{array}$ (9)

公式3到9可以具体求解出关节变量，并且针对关节变量约定$C_i=\cos {\theta }_i$ ，$S_i=\sin {\theta }_i$

$C_{12}=\cos \left({\theta }_1+{\theta }_2\right)$ ，$S_{12}=\sin \left({\theta }_1+{\theta }_2\right)$ 。其中以上公式中相关参数表示如下：

$\begin{array}{c}K=\frac{p_x^2+p_y^2+p_z^2+a_1^2-2a_1{C}_1{p}_x-2a_1{S}_1{p}_y-a_3^2-d_4^2-a_2^2}{2a_2}\end{array}$ (10)

$\begin{array}{c}\rho =\sqrt{a_3^2+d_4^2}\end{array}$ (11)

为进行正逆向运动学的验证，随机设定关节1至关节6的角度为20˚，−100˚，140˚，35˚，80˚，42˚，由此根据式1和2可以得出机器人的末端姿态矩阵${}_6{}^{0}T$ 如下：

${}_6{}^{0}T=\left[\begin{array}{cccc}0.5716& -0.6039& -0.5555& 506.7534\\ -0.4540& -0.7967& 0.3989& 279.4192\\ -0.6835& 0.0241& -0.7296& 771.9144\\ 0& 0& 0& 1\end{array}\right]$

已知末端姿态矩阵${}_6{}^{0}T$ ，根据公式3至11可以反向求解机器人的关节变量${\theta }_1$ 至${\theta }_6$ 。由于机器人关节2可能存在两组解，同时翻转腕关节可以在关节4~6产生新解，以及关节4和6旋转范围为±350˚，求解得出所有满足${}_6{}^{0}T$ 的关节角共16组如下表2：

Table 2. The solution results of the inverse kinematics of the robot

表2. 机器人运动学逆解求解结果

关节角	${\theta }_1(˚)$	${\theta }_2(˚)$	${\theta }_3(˚)$	${\theta }_4(˚)$	${\theta }_5(˚)$	${\theta }_6(˚)$
第1组解	20.0000	−100.0000	140.0000	35.0000	80.0000	42.0000
第2组解	20.0000	−100.0000	140.0000	−145.0000	−80.0000	−138.0000
第3组解	20.0000	−100.0000	140.0000	−325.0000	80.0000	42.0000

续表

第4组解	20.0000	−100.0000	140.0000	215.0000	−80.0000	−138.0000
第5组解	20.0000	−100.0000	140.0000	35.0000	80.0000	−318.0000
第6组解	20.0000	−100.0000	140.0000	−145.0000	−80.0000	222.0000
第7组解	20.0000	−100.0000	140.0000	−325.0000	80.0000	−318.0000
第8组解	20.0000	−100.0000	140.0000	215.0000	−80.0000	222.0000
第9组解	−160.0000	−119.8252	−97.4304	−145.2451	82.2515	43.5882
第10组解	−160.0000	−119.8252	−97.4304	34.7549	−82.2515	−136.4118
第11组解	−160.0000	−119.8252	−97.4304	214.7549	82.2515	43.5882
第12组解	−160.0000	−119.8252	−97.4304	−325.2451	−82.2515	−136.4118
第13组解	−160.0000	−119.8252	−97.4304	−145.2451	82.2515	−316.4118
第14组解	−160.0000	−119.8252	−97.4304	34.7549	−82.2515	223.5882
第15组解	−160.0000	−119.8252	−97.4304	214.7549	82.2515	−316.4118
第16组解	−160.0000	−119.8252	−97.4304	−325.2451	−82.2515	223.5882

3. 拓扑优化设计

运动学仅仅考虑机构的运动关系，并不会研究每个连杆之间的受力情况，真实的机械臂在进行切削量较小的切削加工时受力较小，现成产品的设计有较多的设计余量，因此考虑针对实际的使用情况进行静力学分析，选取受力最小的构件进行拓扑优化设计。

3.1. 机械臂静力学分析

当机械臂被用于铣削、滚压加工时，在机械臂末端固定电主轴后便可以装夹铣削滚压刀具，进行铣削滚压加工[7]。根据加工对象、刀具种类、切削用量可以计算出最终传递到机械臂上的加工反力。铣削加工时刀具材料为硬质合金，被加工材料为6061铝合金，铣削力一般可以采用以下经验公式计算：

$\begin{array}{c}{F}_c=\frac{9.81\times 789.3}{4}{a}_{p}Z{f}_z^{0.75}{a}_e^{1.1}{d}_0^{-1.3}{n}^{-0.2}\end{array}$ (12)

公式12中，$F_c$ ——圆周铣削分力、$a_p$ ——铣削深度、$Z$ ——铣刀齿数、$f_z$ ——每齿进给量、$a_e$ ——铣削宽度、$d_0$ ——铣刀直径、$n$ ——电主轴转速。

实际加工中铣削参数为：$a_p=0.5\text{ mm}$ ；$Z=4$ ；$f_z =0.15\text{ mm}$ ；$a_e =3\text{ mm}$ ；$d_0 =6\text{ mm}$ ；$n=15600\text{ r/min}$ 。带入参数可得圆周铣削力$F_c=264.67\text{ N}$ ，将铣削力分解可得到垂直进给反力$F_v=0.55F_c=145.57\text{ N}$ ，横向进给反力$F_e=0.35F_c=92.63\text{ N}$ ，纵向进给反力$F_f=0.9F_c=238.20\text{ N}$ 。加工过程中这些力作用在刀具上并传递到机械臂上。

根据以上条件对机械臂进行有限元分析，计算得到的滚压力作用在机器人末端，机器人整体与地面连接，直接设定为底座固定。机器人关节设定在垂直于加工工件表面的角度，固定各个关节，随后进行网格划分，完成求解。

图3为分析后得到的等效应力分布云图。从图中可以看出机械臂的连杆2在此加工情况下产生的应力最大不超过21.5 MPa，远远低于机械臂材料的许用应力180 MPa，因此可以考虑对连杆2进行拓扑优化，优化其结构，以减轻质量，节约工件成本。

Figure 3. Finite element analysis result of robotic arm

图3. 机械臂有限元分析结果

3.2. 拓扑优化理论基础

拓扑优化的过程一般是将设计区域离散成足够小的子设计区域，然后对这若干个子设计区域进行结构优化，即按照某种优化策略和准则从这若干个子设计区域中删除某些区域，并使用保留下来的区域表达结构的最优拓扑[8]。

对机械臂部件的优化首先要确定优化目标，然后建立目标函数，对优化问题进行求解。本文的目标是优化机械臂的质量，下面以简单模型建立优化设计问题的数学模型。此种问题一般将设计区域离散成N个带孔洞的结构，如图4所示，对于整体设计区域的质量可以表示为：

(13)

以孔洞边长为设计变量$a_i$ ，整个结构的质量W为目标函数，位移、变形、应力的边界条件作为约束条件$g_j$ ，则求解质量最小的优化问题数学模型可以表述如下：

(14)

Figure 4. Simplified model of topology optimization theory

图4. 拓扑优化理论简化模型

3.3. 连杆拓扑优化

从图3中可知连杆2在铣削时产生的应力非常小，设计可能过于保守，可以进行拓扑优化。目前市面上常用的优化设计软件有Ansys，HyperWorks等，Ansys软件是将结构、流体、电场、磁场、声场等领域分析融合为一体的大型通用有限元分析商用软件，由美国Ansys公司开发。本次拓扑优化直接利用Ansys Workbench中的Topology Optimization计算，得到优化后模型如图5，其中空心部分(红圈)是建议去除的材料。根据建议去除的材料，从中可以发现软件计算结果去除的部分较多为中间部分，因此可以考虑将连杆2增加凹槽特征，去除部分重量。

重新设计连杆2后模型对比如图6所示，图中均为半剖图，左部分为原来模型，右部分为优化后模型，主要在原本连杆2的上下两个面增加凹槽特征去除部分材料。连杆2在优化前质量为97.67 kg，优化后质量为68.79 kg，对比优化前后前三阶固有频率，具体数据见表2，发现优化后前两阶固有频率有较小的提高，第三阶固有频率下降。

Table 3. Natural frequency comparison before and after optimization

表3. 优化前后固有频率对比

阶数	优化前/Hz	优化后/Hz	变化/Hz
1	128.26	134.07	5.81
2	210.27	211.95	1.68
3	668.57	513.19	−155.38

Figure 5. Topology optimization result

图5. 拓扑优化计算结果

Figure 6. Link 2 redesign schematic

图6. 连杆2重新设计示意图

以上的拓扑优化设计结果可以为加工制造机械臂提供参考思路，但在去除材料的过程中仅仅使用凹槽特征较为单一，实际的拓扑优化过程非常复杂，通常需要不断的迭代才能得到较为合适的结果。

4. 结语

本文通过理论分析，使用DH参数法构建了KUKA机械臂正逆向运动学模型，能够通过关节参数计算机械臂末端位置，以及给定机械臂末端位置反算对应6个关节需要的运动角度。该方法对关节型机器人的运动学模型建立与推导有较强的参考性。在对机械臂连杆2的拓扑优化中，从优化结果分析，重新建模了连杆2的模型，增加了凹槽特征，去除了28.8 kg的质量。该优化设计方法实现了去除工件多余质量，同时保证了工件使用性能不受影响，但还需要不断迭代最优结果，并且进行试验验证。