一种机器人铣削响应动态仿真方法

技术领域

本发明涉及机器人铣削领域，具体涉及一种机器人铣削响应动态仿真方法。

背景技术

机器人铣削技术作为一种前景广阔的智能制造技术和先进制造技术，具有加工空间大、相对体积小、位形灵活、成本较低等特点，并且在加工大尺寸复杂形状零件时表现出突出的优势，引起了行业内人员的广泛关注。与此同时，机器人铣削技术也存在不可忽视的缺陷，它的稳定性较差，对加工过程中的颤振现象非常敏感；机器人铣削颤振会降低零件表面加工质量，降低加工效率，甚至导致零件成批报废。为实现铣削过程颤振抑制，预测机器人铣削稳定性很有必要，这种预测可以通过机器人铣削响应仿真实现。

机器人铣削响应仿真算法涉及机器人结构动力学模型、铣削力模型以及两种模型之间的耦合，算法设计难度较大。现有技术中，例如在机器人结构动力学方面，存在有基于经典的刚柔耦合结构动力学模型提出的机器人刚柔耦合动力学的残差形式，旨在快速实现机器人振动仿真；再例如在铣削力建模方面，存在有以经典的Altinta铣削力模型为基础提出的Dystamill铣削力仿真程序，旨在实现铣削力的高效仿真。

在机器人动力学建模方面，现有的机器人结构动力学或性能较差，或结构臃肿；在系统控制领域，机器人模型通常为刚体动力学模型，在振动研究领域，机器人通常被建模为模态参数结构动力学模型；动力学领域缺乏综合考虑机器人刚性属性与结构柔性，使用机器人刚柔耦合动力学模型描述机器人振动状态的动力学模型。

现有机器人刚柔耦合动力学多为内容臃肿的微分矩阵方程，系数项较多，影响因素复杂，求解效率低下；在铣削力建模方面，现有整体立铣刀铣削力模型内容经典，形式简洁，但算法结构复杂，应用难度较大。Dystamill仿真程序使用计算机图形学理论实现整体立铣刀与工件的建模与离散，程序的编写与维护相对困难；使用C++语言编写底层程序并执行仿真，仿真结果需要转移至Matlab平台进行可视化处理与数据分析，算法的应用方式相对繁琐；算法通过几何属性计算未变形切屑厚度，并保存加工后工件的实时形状，消耗的算力资源与存储资源相对较多。此外，上述两种模型之间的耦合算法仍有待研究。

发明内容

针对现有技术的不足，本发明的主要目的是提供一种机器人铣削响应动态仿真方法，来实现机器人动力学模型与铣削力模型之间的耦合，并能够提供高质量的机器人铣削响应仿真数据。具有耦合效率高、仿真效果好的优点。

为了实现上述主要目的，本发明提供了一种机器人铣削响应动态仿真方法，其包括以下步骤：

S1：建立机器人刚柔耦合动力学模型，并对机器人刚柔耦合动力学模型进行解耦，得到机器人结构动力学模型；

S2：将螺旋刃立铣刀沿刀具轴向离散为一系列直刃铣刀片的叠加；基于Altinta铣削力模型理论建立直刃铣刀片的铣削力模型，描述铣削力与铣削过程中切下的未变形切屑厚度之间的线性关系；

S3：将仿真时间离散为大量的时间微元；通过机器人结构动力学预测当前时间微元结束时刻机器人末端的实际位置，根据机器人末端在当前时间微元结束时刻的理论位置来获取各刀片在当前时间微元的未变形切削厚度；通过步骤S2计算得到各刀片的铣削力，叠加得到该时间微元的总铣削力；

S4：将步骤S3得到的总铣削力作为机器人的外力，计算机器人在当前时刻的末端振动，结合末端理论位置得到末端修正位置；经过振动修正的末端修正位置与之前的末端理论位置之间存在误差，使用刚刚得到的末端修正位置作为新的预测值并再次计算新的末端修正位置，直至相邻次末端修正位置的计算结果足够接近时，取当前的末端修正位置为该时间微元结束时刻的实际末端位置；

S5：将当前时间微元结束时刻的数据赋值给下一时间微元的开始时刻，在下一时间微元中继续循环执行步骤S4，直至仿真结束并得到机器人铣削响应力信号。

根据本发明的一种具体实施方式，步骤S1展开如下：

S11：建立机器人刚体动力学模型为：

其中，广义坐标q是表示各关节输出角度的6维向量，惯性矩阵M用于描述系统惯性力对运动的影响，克里斯托弗尔项C表示离心力与科里奥利力对系统运动的作用；粘滞摩擦与滑动摩擦F

S12：将关节轴等效具有一定动刚度的扭转弹簧模型，并完善机器人刚体动力学，得到的机器人刚柔耦合动力学形式为：

其中，广义坐标q

S13：将机器人系统视为刚性系统与柔性系统的叠加，并将机器人的动力学解耦为刚性动力学与结构动力学，得到的结构动力学模型为：

其中，M为惯性矩阵，广义坐标q是表示各关节输出角度的6维向量，C为克里斯托弗尔项、D为关节轴阻尼矩阵、K为刚度矩阵。

根据本发明的一种具体实施方式，步骤S3中，离散后某一刀片上的某一刀齿所受到的铣削力在铣削加工平面坐标系中表示为：

其中，k与j为刀片和刀齿的编号；dh、dz和ds分别为未变形切屑厚度、离散后的刀片厚度和局部切削刃长度；系数K

铣削加工平面坐标系用于描述铣刀状态与切削情况，该铣削加工平面坐标系的位姿与工件固连，其原点与机器人惯性坐标系原点坐落在同一竖直平面内，该平面与惯性坐标系X轴重合。

根据本发明的一种具体实施方式，在步骤S3中所获得的铣削力位于切向-法向-轴向坐标系，可以通过如下方式向机器人惯性坐标系转换：

其中，F

根据本发明的一种具体实施方式，在步骤S1中需要先确定机器人的几何参数与动力学参数，其包括D-H参数、杆件惯性参数、关节柔性参数以及机器人初始位形。

根据本发明的一种具体实施方式，在步骤S2中需要先确定螺旋刃立铣刀的几何参数和切削力系数；其中螺旋刃立铣刀的几何参数包括刀具尺寸、刀具齿数、刀具进给速度以及刀具螺旋角。

本发明具备以下有益效果：

本发明可实现机器人刚柔耦合结构动力学与铣削力模型之间静态的数据交换和动态的状态更新过程；仿真过程动态进行，结构动力学模型与铣削力模型同时参与仿真，保证了仿真结果的有效性；另外，本发明还具有数据交换与状态更新速度快、仿真效率高的优点。

为了更清楚地说明本发明的目的、技术方案和优点，下面结合附图和具体实施方式对本发明作进一步的详细说明。

附图说明

图1是机器人关节模型构成和运动传递过程的示意图；

图2是仿真过程中常用的惯性参考系和铣削加工平面参考系的示意图；

图3是机器人铣削响应动态仿真过程的流程图；

图4是机器人铣削响应动态仿真结果与现有技术的对照图。

具体实施方式

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不限于下面公开的具体实施例的限制。

本发明实施例提供了一种机器人铣削响应动态仿真方法，其包括以下步骤：

S1：建立机器人刚柔耦合动力学模型，并对机器人刚柔耦合动力学模型进行解耦，得到机器人结构动力学模型；

其中，需要先确定机器人的几何参数与动力学参数，其包括D-H参数、杆件惯性参数、关节柔性参数以及机器人初始位形。

(1)机器人D-H参数；

如表1所示，在模型中，机器人底座并未被视为结构的一部分；

为保证数据单位的统一，在动态仿真算法中将长度与角度单位统一为m和rad。

表1：机器人D-H参数辨识结论

(2)机器人惯性参数。

使用制造商提供的几何模型与质量参数，结合商用CAD软件，识别本发明实施例机器人中各构件惯性参数，结果如表2所示。

表中每个连杆的惯性参数包括质量、质心位置与中心转动张量，参考系指上述惯性参数的观察坐标系，也是几何模型中的关节坐标系。

表2：机器人惯性参数辨识结论

(3)机器人工作位形；

本发明实施例中机器人工作位形位于[–10.09,9.36,17.81,21.30,–28.85,–18.86]附近，单位为deg，需要转化为rad后再参与仿真。

(4)柔性参数矩阵；

通过虚拟实验模态分析可以辨识系统柔性参数矩阵。常用的柔性参数矩阵为对角矩阵，考虑到机器人系统动态特性的非对称性，以及各关节变形之间的耦合特性，参考系统振型，给出一种非对称形式的柔性参数矩阵形式，辨识得到的阻尼矩阵D和刚度矩阵K如表达式(Ⅰ)和表达式(Ⅱ)所示；使用这两种柔性参数矩阵和惯性矩阵，执行机器人铣削仿真的结果收敛。

进一步地，步骤S1展开如下：

S11：建立机器人刚体动力学模型为：

其中，广义坐标q是表示各关节输出角度的6维向量，惯性矩阵M用于描述系统惯性力对运动的影响，克里斯托弗尔项C表示离心力与科里奥利力对系统运动的作用；粘滞摩擦与滑动摩擦F

S12：将关节轴等效具有一定动刚度的扭转弹簧模型，并完善机器人刚体动力学，得到的机器人刚柔耦合动力学形式为：

其中，广义坐标q

S13：将机器人系统视为刚性系统与柔性系统的叠加，并将机器人的动力学解耦为刚性动力学与结构动力学，得到的结构动力学模型为：

其中，M为惯性矩阵，广义坐标q是表示各关节输出角度的6维向量，C为克里斯托弗尔项、D为关节轴阻尼矩阵、K为刚度矩阵。

S2：将螺旋刃立铣刀沿刀具轴向离散为一系列直刃铣刀片的叠加；基于Altinta铣削力模型理论建立直刃铣刀片的铣削力模型，描述铣削力与铣削过程中切下的未变形切屑厚度之间的线性关系；

其中，需要先确定螺旋刃立铣刀的几何参数和切削力系数；螺旋刃立铣刀的几何参数包括刀具尺寸、刀具齿数、刀具进给速度以及刀具螺旋角。在本发明实施例中，定义刀具半径为10mm，刀具螺旋角为30°，工作状态为半浸没逆铣，进给速度为2.4mm/s，仿真总时间为1s。具体的，使用三菱硬质合金螺旋刃整体立铣刀铣削6082铝制工件，通过逆向求解方法识别铣削力系数，得到的铣削力系数和边缘铣削系数分别为1197.7、3857.5、-518.0和-1.1、-20.4、2.4，上述参数可以保证仿真结果收敛。

在确定上述参数的基础上，结合图1建立机器人结构动力学模型，以及建立铣削力模型。

S3：将仿真时间离散为大量的时间微元；通过机器人结构动力学预测当前时间微元结束时刻机器人末端的实际位置，根据机器人末端在当前时间微元结束时刻的理论位置获得各铣刀片上各刀齿的位置；对当前时刻的刀齿位置与前一刀齿切削轨迹进行比较，将两者之间的距离记为未变形切屑厚度；通过步骤S2计算得到各刀片的铣削力，叠加得到该时间微元的总铣削力；

在步骤S3中，离散后某一刀片上的某一刀齿所受到的铣削力在铣削加工平面坐标系中表示为：

其中，k与j为刀片和刀齿的编号；dh、dz和ds分别为未变形切屑厚度、离散后的刀片厚度和局部切削刃长度；系数K

铣削加工平面坐标系用于描述铣刀状态与切削情况，该铣削加工平面坐标系的位姿与工件固连，其原点与机器人惯性坐标系原点坐落在同一竖直平面内，该平面与惯性坐标系X轴重合。

在步骤S3中所获得的铣削力通过如下方式向机器人惯性坐标系转换：

其中，F

S4：将步骤S3得到的总铣削力作为机器人的外力，计算机器人在当前时刻的末端振动，结合末端理论位置计算末端修正位置；经过振动修正的末端修正位置与之前的末端理论位置之间存在误差，使用刚刚得到的末端修正位置作为新的预测值并再次计算新的末端修正位置，直至相邻次末端修正位置的计算结果间距小于1×10

S5：将当前时间微元结束时刻的数据赋值给下一时间微元的开始时刻，在下一时间微元中继续循环执行步骤S4，直至仿真结束并得到机器人铣削响应力信号。

在本发明实施例中，根据图2定义机器人参考系之后，按照图3所示的流程开展机器人仿真，操作参数为主轴转速取2000rpm至8000rpm，间隔100rpm，产生61种主轴转速；轴向切削深度为0.5mm至2.2mm，间隔0.1mm，产生18种轴向切削深度。使用这些操作参数进行仿真，共得到61×18＝1098组仿真结果。

在本发明实施例仿真结束后，使用能量熵稳定性判据处理仿真结果，得到的系统稳定性分布情况如图4中等高线图所示，红色线条为临界稳定曲线。图4中的红色三角、黄色方框和绿色圆圈分别代表实验中表现为颤振、临界稳定和稳定性的加工参数，其分布与仿真结果匹配情况好；青色和蓝色曲线分别为使用经典的全离散方法得到的稳定性边界线及其优化后的曲线，两条曲线与实验结果之间的匹配程度不如本算法仿真结果。此外，仿真结果中的熵值凸起与临界稳定成分现象与实验结果相对应。上述结论均说明，本发明实施例的机器人铣削响应动态仿真方法具有较高的仿真质量，可以实现机器人铣削稳定性的高质量预测。

虽然本发明以具体实施例揭露如上，但这些具体实施例并非用以限定本发明实施的范围。任何本领域的普通技术人员，在不脱离本发明的发明范围内，当可作些许的变化/修改，即凡是依照本发明所做的同等变化/修改，应为本发明的保护范围所涵盖。