CN116619357A - 一种铣削机器人变姿态稳定性叶瓣图获取方法 - Google Patents

Abstract

本发明提供了一种铣削机器人变姿态稳定性叶瓣图获取方法，包括：将机器人的工作空间网格化，在划分的网格中进行模态试验，辨识机器人不同姿态下的关节刚度和阻尼，建立铣削机器人变姿态动力学模型，获取机器人固有振动属性和刀尖频响函数，基于零阶频域稳定性算法绘制机器人在大范围工作空间下变姿态铣削的稳定性叶瓣图。本发明所提出的方法充分考虑机器人大场景作业时姿态多变对铣削稳定性的影响，结合多体动力学方法与零阶频域法预测铣削稳定域，能够准确、快速地选取稳定铣削时的工艺参数(转速、切深)，从而避免加工过程中的颤振问题。

Description

一种铣削机器人变姿态稳定性叶瓣图获取方法

技术领域

本发明属于工业机器人铣削稳定性预测领域，尤其涉及一种铣削机器人变姿态稳定性叶瓣图获取方法。

背景技术

随着航空航天产业的蓬勃发展，许多工作场景下需要工业机器人进行原位加工。机器人作为铣削加工载体时比传统机床更灵活、更高效，能够在环境复杂的加工现场进行工件的铣削作业，因而在加工智能化、柔性化领域发展前景广阔。

然而，机器人本体的串联结构刚性较差，铣削加工时末端刀具与工件容易发生强相互作用而引起机器人-刀具系统颤振，末端刀具剧烈抖动产生加工噪声甚至造成刀具和工件损伤，严重影响加工精度和加工质量。稳定性叶瓣图是预防机器人铣削加工颤振的重要手段，能够通过选择稳定铣削参数，进而避免颤振现象的发生。传统稳定性叶瓣图主要应用于机床颤振，而有关机器人颤振模型的研究仍不够成熟。

Said Mousavi等在学术期刊《International Journal of AdvancedManufacturing Technology》，2016，88(9-12):3053-3065上发表的文章中建立了基于梁单元的串联机器人多体动力学模型，通过实验识别梁单元的几何、弹性和阻尼参数，将铣削过程中再生颤振作为运动冗余度变量函数，实现了稳定性叶瓣图绘制。

Marcel Cordes等在学术期刊《Robotics and Computer IntegratedManufacturing》2019，55：11-18上发表的文章中建立了具有主轴和工具的铰接式机械手的结构动力学模型，用半离散的时间域和频域方法分析了系统的稳定性。在铝合金和钛合金零件的铣削加工实验中验证了稳定性叶瓣图的准确性。

上述研究建立了机器人的动力学模型并通过时域、频域等方法较为准确地实现了机器人铣削稳定性的预测，但其不足之处在于没有考虑到机器人姿态变化对模态参数与稳定性状态的影响，难以在大空间、复杂加工姿态下建立准确的稳定性预测模型。

发明内容

发明目的：本发明所要解决的技术问题是针对现有技术的不足，提供一种铣削机器人变姿态稳定性叶瓣图获取方法。考虑机器人姿态变化对关节刚度和阻尼的影响，基于多体系统传递矩阵法建立机器人变姿态铣削动力学模型，预测系统固有频率与刀尖点频响函数，结合铣削力模型进行颤振稳定域解析计算，绘制铣削机器人变姿态稳定性叶瓣图。为实现上述目的，本发明采用以下技术方案：

一种铣削机器人变姿态稳定性叶瓣图获取方法，包括以下步骤：

步骤1，基于多体系统传递矩阵法建立机器人动力学模型，同时建立铣削过程切削力模型，形成铣削机器人动力学仿真系统；

步骤2，将机器人铣削工作空间标准化为一个大立方体，设定大立方体体积为P，进一步将空间划分成Z个边长为h的小立方体网格，Z＝P/h³，Z取整数，网格编号为i，i＝1,2,3…Z；

步骤3，设置初始网格边长为h和尺寸变化步长，设置网格划分的次数为r，其中，r＝0,1,2…，r＝0表示未划分网格时的初始状态；

步骤4，控制机器人末端刀具中心点处于第i个待采样网格，通过计算机器人运动学逆解获取待采样网格对应姿态下的关节转角，参考崔光裕等在学术期刊《AppliedMathematical Modelling》2022,105：114-136上发表的文章中机器人关节刚度和阻尼辨识方法，对机器人进行模态测试，辨识待采样网格内机器人关节刚度矩阵K_i和阻尼矩阵C_i；

步骤5，将关节转角、关节刚度矩阵K_i和阻尼矩阵C_i输入铣削机器人动力学仿真系统，计算系统固有频率与刀尖点频响函数；

步骤6，基于零阶频域法建立颤振稳定域预测模型，绘制第i个待采样网格的稳定性叶瓣图，记录绝对稳定区的临界切深；

步骤7，统计并记录所有网格稳定性叶瓣图临界切深的平均值CA_r和最大值CM_r，当r＝0时，有CA_r＝0.01和CM_r＝0.01；

步骤8，将网格尺寸减小X₁(一般取值为50mm)，重新划分网格，重复步骤4～步骤7；

步骤9，当r>0时，计算第r次和第r-1次临界切深平均值变化率x＝|CA_r-CA_r-1|/CA_r以及最大值变化率y＝|CM_r-CM_r-1|/CM_r，当x、y均小于X₂(一般取值为10％)时，输出当前网格划分下的稳定性叶瓣图；否则重复步骤8，直至满足x、y均小于X₂的条件时，输出变姿态稳定性叶瓣图。

步骤1包括：

步骤1-1，建立铣削过程切削力模型；

步骤1-2，测量各元件动力学参数，拼接各元件传递矩阵，推导机器人动力学模型传递方程，代入边界条件求解特征方程，获得固有振动属性；

步骤1-3，根据解耦后的机器人动力学模态方程，将铣削过程切削力模型输入，基于模态叠加原理获取物理空间内振动位移，完成机器人动力学模型建立步骤；

步骤1-4，将铣削过程切削力模型输入机器人动力学模型形成铣削机器人动力学仿真系统，对铣削力和振动位移进行傅里叶变换，计算刀尖点频响函数。

步骤1-1包括：机器人二自由度的瞬时铣削过程切削力模型表示为：

其中，dF_t、dF_r分别表示切向铣削力微元和径向铣削力微元，K_tc、K_rc分别表示切向力系数和径向力系数，a_c为铣削厚度，da为轴向铣削微元，K_te、K_re分别表示切向刃口力系数和径向刃口力系数，ds为切削刃长度微元；

将径向力F_r、切向力F_t投影到工件坐标系，得到瞬时铣削力：

其中，F_x、F_y分别表示工件坐标系下沿x方向的铣削力、沿y方向的铣削力，φ为铣削刀具接触角，φ_ex、φ_st分别表示铣削切出角和切入角，N为铣削刀具刃数，a为轴向切削深度，f为每齿进给量；

对式(2)积分获得x方向平均铣削力y方向平均铣削力/>

通过铣削正交试验获得不同铣削参数下的平均铣削力，结合式(3)辨识K_tc、K_rc。

步骤1-2包括：以工业六自由度串联机器人为例，机器人各连杆通过关节处齿轮连接，末端法兰携带能够铣削的电主轴及刀具；将机器人(各种机器人通用实施例为KUKA-KR210)整体划分为体元件和铰元件，编号为0～14，其中，2、4、6、8、10、12号为机器人本体连杆，建模为刚体元件；14号铣削刀具刚性较差，建模为空间振动Euler-Bernoulli梁元件；各体元件具有独立的随元件移动的连体坐标系；1、3、5、7、9、11、13号为机器人关节处添加的沿连体坐标系x、y、z轴向具有线弹性，绕x、y、z轴具有扭转弹性的6自由度弹性铰元件，铰元件坐标系x、y、z轴方向与机器人基坐标系方向保持一致；0号元件为机器人基座，固定于地面，视为机器人动力学模型的固定边界；14号元件为末端刀具，视为系统的自由边界；整个机器人动力学模型为树形结构，其传递方向由根0到梢14，模型总传递方程为：

Z_14,0＝U_1-14Z_0,1 (4)

其中，为体元件状态矢量，Z_0,1、Z_14,0分别表示与大地相连元件的状态矢量、自由边界相连元件的状态矢量，U_1-14为总传递矩阵，由各个元件传递矩阵依次拼接得到，X、Y、Z分别表示在体元件坐标系下x、y、z方向模态空间下的线位移，Θ_x、Θ_y、Θ_z分别表示在体元件坐标系下x、y、z方向模态空间下的角位移，M_x、M_y、M_z分别表示在体元件坐标系下x、y、z方向模态空间下的内力矩，Q_x、Q_y、Q_z分别表示在体元件坐标系下x、y、z方向模态空间下对应的力；下标b、j分别为体元件的编号和铰元件的编号，固定边界和自由边界的编号为0；

机器人各元件传递矩阵拼接得到机器人动力学模型总传递矩阵：

U_1-14＝U₁₄U₁₃U₁₂U₁₁U₁₀U₉U₈U₇U₆U₅U₄U₃U₂U₁ (5)

其中，U为模型中各组成元件的传递矩阵，U₁₄为第14个元件的传递矩阵；

根据多体系统传递矩阵法链式推导法则，1号元件与固定边界相连，1号元件的线位移、角位移为0；14号元件与自由边界相连，14号元件的内力、内力矩为0，有如下边界条件：

对于b＝2、4、6、8、10、12号刚体元件，C_b为编号为b连杆的质心点，刚体元件传递矩阵U_b在以输入点I_b为原点，O_b为输出点的自身连体系中描述为：

式中，

其中，I₃为3行3列的单位矩阵，O_3×3为3行3列的零矩阵，m_b为第b个体元件质量，为第b个体元件坐标系中输出点O_b相对于输入点I_b的坐标反对称矩阵，/>为第b个体元件坐标系中质心点C_b相对于输出点O_b的坐标反对称矩阵，/>为第b个体元件坐标系中质心点C_b相对于输入点I_b的坐标反对称矩阵，ω_p为系统第p阶固有频率，J_b,I为刚体相对于输入点I_b的惯量矩阵，J_b,x、J_b,y、J_b,z分别为以输入点I_b为原点的元件坐标系中x、y、z方向的转动惯量，J_b,xy、J_b,xz、J_b,yz分别为空间中x、y、z方向转动惯性积，(g_b1,g_b2,g_b3)为输出点坐标，(c_b1,c_b2,c_b3)为质心坐标；

对于b＝14号梁元件，其传递矩阵为：

式中，

u_1,1＝u_10,10＝cosβ_xx_b，u_10,1＝β_xEA_bsinβ_xx_b，u_4,4＝u_7,7＝cosγ_θxx_b， u_6,2＝u_11,9＝λ_yV(λ_yx_b)，/>u_5,3＝u_12,8＝-λ_zV(λ_zx_b)，/> u_8,3＝u_12,5＝-EI_yλ_z ²U(λ_zx_b)，/>u_2,2＝u_6,6＝u_9,9＝u_11,11＝S(λ_yx_b)，u_8,5＝-EI_yλ_zV(λ_zx_b)，u_9,6＝EI_zλ_yV(λ_yx_b)，u_3,3＝u_5,5＝u_8,8＝u_12,12＝S(λ_zx_b)，

其中，中间变量λ_y、γ_θx、β_x、λ_z的计算公式为：

其中，中间函数变量S(u)、T(u)、U(u)、V(u)的计算公式为：

其中u为函数自变量；

其中，传递矩阵U₁₄中参数在自身连体坐标系I₁₄中描述，下标b表示元件序号，为线质量密度，x_b表示梁的长度，ρ为梁的密度，E为弹性模量，G为剪切模量，A_b为横截面积，I_x为横截面对x轴的转动惯量，I_y为横截面对y轴的转动惯量，I_z为横截面对z轴的转动惯量，J_p为横截面对x轴的极惯性矩；

所有刚体元件和梁元件的传递矩阵参数通过三维CAD软件测量获取，且不随网格变化而改变；

对于j＝1、3、5、7、9、11、13号铰元件，其传递矩阵U_j为：

式中，

其中，K″_j表示编号为j的铰元件的线刚度系数矩阵，K′_j表示编号为j的铰元件的扭转刚度系数矩阵，K″_j,x、K″_j,y、K″_j,z分别表示编号为j的铰元件x、y、z方向的线弹性刚度，K′_j,x、K′_j,y、K′_j,z分别表示x、y、z方向的扭转弹性刚度，第i个网格内的关节刚度矩阵K_i包含该网格下所有铰元件的线刚度及扭转刚度信息，需要通过试验辨识获取；

将边界条件代入总传递方程进一步推导系统的特征方程：

det[U_1-14(ω_p)]＝0 (11)

其中，det[u]表示u的行列式函数，u为函数自变量，ω_p为机器人第p阶固有频率；求解式(11)得到机器人各阶固有频率。

步骤1-3包括：通过模态试验进行第i个网格下阻尼矩阵C_i参数辨识，解耦后获取模态阻尼比，根据式(12)计算机器人末端的振动位移响应：

其中，ζ_p为系统第p阶模态阻尼比，q^p(t)、分别为第p阶广义位移、广义速度和广义加速度，V^p为第p阶增广特征矢量，M_p＝<MV^p,V^p>为系统第p阶模态质量，M为系统质量矩阵，f为外力列阵。

步骤1-4包括：将铣削过程切削力力模型作为外部激励输入，基于模态叠加原理将各阶模态响应叠加，计算物理空间振动位移输出，对输入力、输出位移分别进行傅里叶变换，求解刀尖点频响函数G(iω)与颤振频率ω_d，频响函数矩阵为：

其中，G_xx(iω)、G_yy(iω)分别为x方向上的频响函数、y方向上的频响函数，i为虚数单位。

步骤5包括：当系统发生颤振时，在颤振频率ω_d处的特征方程为：

其中，I₂为2行2列的单位矩阵，t为时间变量，A₀为随时间变化的总切削力方向系数矩阵的傅里叶级数展开的第1项；

特征方程的特征值Λ为：

其中，T为切削周期，Λ_R与Λ_I分别为系统特征值的实部和虚部；

特征值满足如下关系式：

b₁Λ²+b₂Λ+1＝0 (16)

式中，b₁、b₂是中间参数，

b₁＝G_xx(iω_d)G_yy(iω_d)(α_xxα_yy-α_xyα_yx)，b₂＝α_xxG_xx(iω_d)+α_yyG_yy(iω_d)，

其中，α_xx、α_yy、α_xy、α_yx为铣削方向系数；

改写式(15)并求解式(16)，求得二自由度铣削系统的临界铣削深度a_lim，对于实际的铣削加工，刀具的轴向切削深度为实数，因此得到：

式中，κ为中间参数，

根据主轴转速n与周期T之间的关系，获得叶瓣图中与式(17)计算的轴向切深对应的转速n为：

式中，k＝1,2,3...10，为铣削稳定性叶瓣图的叶瓣数；

根据实际过程中所需转速范围选择叶瓣数，以主轴转速为横轴，临界轴向铣削深度为纵轴，绘制稳定性叶瓣图。

步骤3中，网格初始尺寸和变化步长可根据机器人实际工作空间确定。

步骤4中，机器人的关节刚度和阻尼随姿态变化发生改变，当末端在同一个网格内运动时引起姿态变化较小，此时认为在该网格内刚度和阻尼近似不变。

步骤4中，进行单一网格内模态试验时，机器人姿态的选取应以实际铣削加工时机器人姿态为参考，进而辨识关节刚度和阻尼。

本发明的有益效果是：

1、本发明提出的通过工作空间网格化建立的铣削机器人变姿态动力学模型，考虑机器人姿态变化引起的关节刚度、阻尼改变，实现机器人不同姿态下的模态特性和动力学响应(位移、速度、加速度)较为快速准确地仿真，能够在保证所需精度的前提下节约试验时间成本；

2、本发明提出的铣削机器人变姿态稳定性叶瓣图获取方法，考虑机器人铣削姿态多变的问题，相较于现有稳定性预测方法，能够更准确、更高效地建立变姿态铣削稳定性预测模型；

3、本发明绘制的变姿态稳定性叶瓣图，针对大范围复杂工作空间内机器人铣削作业，可以在不同姿态下即时提供合适的铣削工艺参数，有效提升机器人铣削加工稳定性，减少颤振现象。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1是本发明方法流程图。

图2是本发明铣削机器人动力学模型拓扑图。

图3是本发明机器人作业空间网格化刚度阻尼辨识过程示意图。

图4是本发明机器人绘制的机器人在姿态Φ₁＝(0,-90,90,0,0,180)°时的稳定性叶瓣图。

图5是本发明机器人绘制的机器人在姿态Φ₂＝(0,-71,90,0,0,180)°时的稳定性叶瓣图。

图6是本发明机器人绘制的机器人在姿态Φ₃＝(-11，-58,95,-96,-95,128)°时的稳定性叶瓣图。

具体实施方式

如图1所示流程，以KUKA-210铣削机器人为例，本发明的一种铣削机器人变姿态稳定性叶瓣图获取方法，包括如下步骤：

1、建立铣削机器人动力学仿真系统；

首先，建立铣削过程切削力模型，完成铣削力系数辨识。

机器人二自由度的瞬时铣削过程切削力模型可以表示为：

式中，dF_t、dF_r分别表示切向铣削力微元和径向铣削力微元，K_tc、K_rc分别表示切向力系数和径向力系数，a_c为铣削厚度，da为轴向铣削微元，K_te、K_re分别表示切向刃口力系数和径向刃口力系数，ds为切削刃长度微元；

将径向力F_r、切向力F_t投影到工件坐标系，得到瞬时铣削力：

式中，F_x、F_y表示工件坐标系下沿x、y方向的铣削力，φ为铣削刀具接触角，φ_ex、φ_st分别表示铣削切出角和切入角，N为铣削刀具刃数，a为轴向切削深度，f为每齿进给量；

对式(2)积分获得x、y方向平均铣削力

通过铣削正交试验获得不同铣削参数下的平均铣削力，结合式(3)辨识K_tc、K_rc。

其次，建立机器人-刀具铣削系统刚柔耦合动力学模型，实现刀尖频响函数计算。

以工业六自由度串联机器人为例，机器人动力学模型结构示意图如图2所示。图2中o为机器人基坐标系原点，x、y、z为基坐标轴，带下标的x、y、z为各元件连体系坐标轴，带下标的I、O、C、K′、K″、U为在各个元件传递矩阵中定义的变量。

机器人各连杆通过关节处齿轮连接，末端法兰携带能够铣削的电主轴及刀具；将机器人整体划分为体元件和铰元件，编号为0～14，其中，2、4、6、8、10、12号为机器人本体连杆，建模为刚体元件；14号铣削刀具刚性较差，建模为空间振动Euler-Bernoulli梁元件；各体元件具有独立的随元件移动的连体坐标系；1、3、5、7、9、11、13号为机器人关节处添加的沿x、y、z轴向具有线弹性，绕x、y、z轴具有扭转弹性的6自由度弹性铰元件，铰元件坐标系x、y、z轴向与机器人基坐标系方向保持一致；0号元件为机器人基座，固定于地面，视为机器人动力学模型的固定边界；14号元件为末端刀具，视为系统的自由边界；整个机器人动力学模型为树形结构，其传递方向由根0到梢14，模型总传递方程为：

Z_14,0＝U_1-14Z_0,1 (4)

其中，为体元件状态矢量，Z_0,1、Z_14,0分别表示与大地相连元件的状态矢量、自由边界相连元件的状态矢量，U_1-14为总传递矩阵，由各个元件传递矩阵依次拼接得到，X、Y、Z分别表示在体元件坐标系下x、y、z方向模态空间下的线位移，Θ_x、Θ_y、Θ_z分别表示在体元件坐标系下x、y、z方向模态空间下的角位移，M_x、M_y、M_z分别表示在体元件坐标系下x、y、z方向模态空间下的内力矩，Q_x、Q_y、Q_z分别表示在体元件坐标系下x、y、z方向模态空间下对应的力；下标b、j分别为体元件的编号和铰元件的编号，固定边界和自由边界的编号为0；

机器人各元件传递矩阵拼接得到机器人动力学模型总传递矩阵：

U_1-14＝U₁₄U₁₃U₁₂U₁₁U₁₀U₉U₈U₇U₆U₅U₄U₃U₂U₁ (5)

其中，U为模型中各组成元件的传递矩阵，U₁₄为第14个元件的传递矩阵；

根据多体系统传递矩阵法链式推导法则，1号元件与固定边界相连，1号元件的线位移、角位移为0；14号元件与自由边界相连，14号元件的内力、内力矩为0，有如下边界条件：

对于b＝2、4、6、8、10、12号刚体元件，C_b为编号为b连杆的质心点，刚体元件传递矩阵U_b在以输入点I_b为原点，O_b为输出点的自身连体系中描述为：

式中，

其中，I₃为3行3列的单位矩阵，O_3×3为3行3列的零矩阵，m_b为第b个体元件质量，为第b个体元件坐标系中输出点O_b相对于输入点I_b的坐标反对称矩阵，/>为第b个体元件坐标系中质心点C_b相对于输出点O_b的坐标反对称矩阵，/>为第b个体元件坐标系中质心点C_b相对于输入点I_b的坐标反对称矩阵，ω_p为系统第p阶固有频率，J_b,I为刚体相对于输入点I_b的惯量矩阵，J_b,x、J_b,y、J_b,z分别为以输入点I_b为原点的元件坐标系中x、y、z方向的转动惯量，J_b,xy、J_b,xz、J_b,yz分别为空间中x、y、z方向转动惯性积，(g_b1,g_b2,g_b3)为输出点坐标，(c_b1,c_b2,c_b3)为质心坐标；

对于b＝14号梁元件，其传递矩阵为：

式中，

u_1,1＝u_10,10＝cosβ_xx_b，u_10,1＝β_xEA_bsinβ_xx_b， u_9,2＝u_11,6＝-EI_zλ_y ²U(λ_yx_b)，u_7,4＝-γ_θx(GJ_p)_isinγ_θxx_b，/> u_6,2＝u_11,9＝λ_yV(λ_yx_b)，/>u_5,3＝u_12,8＝-λ_zV(λ_zx_b)，/> u_8,3＝u_12,5＝-EI_yλ_z ²U(λ_zx_b)，/>u_2,2＝u_6,6＝u_9,9＝u_11,11＝S(λ_yx_b)，u_8,5＝-EI_yλ_zV(λ_zx_b)，u_9,6＝EI_zλ_yV(λ_yx_b)u_3,3＝u_5,5＝u_8,8＝u_12,12＝S(λ_zx_b)

其中，

为方便计算定义的中间变量

其中，

为定义的中间函数变量，u为函数自变量；

其中，传递矩阵U₁₄中参数在自身连体坐标系I₁₄中描述，下标b表示元件序号，为线质量密度，x_b表示梁的长度，ρ为梁的密度，E为弹性模量，G为剪切模量，A_b为横截面积，I_x为横截面对x轴的转动惯量，I_y为横截面对y轴的转动惯量，I_z为横截面对z轴的转动惯量，J_p为横截面对x轴的极惯性矩；

所有刚体元件和梁元件的传递矩阵参数通过三维CAD软件测量获取，且不随网格变化而改变；

对于j＝1、3、5、7、9、11、13号铰元件，其传递矩阵U_j为：

式中，

其中，K″_j表示编号为j的铰元件的线刚度系数矩阵，K′_j表示编号为j的铰元件的扭转刚度系数矩阵，K″_j,x、K″_j,y、K″_j,z分别表示编号为j的铰元件x、y、z方向的线弹性刚度，K′_j,x、K′_j,y、K′_j,z分别表示x、y、z方向的扭转弹性刚度，第i个网格内的关节刚度矩阵K_i包含该网格下所有铰元件的线刚度及扭转刚度信息，需要通过试验辨识获取；

将边界条件代入总传递方程进一步推导系统的特征方程：

det[U_1-14(ω_p)]＝0 (11)

其中，det[u]表示u的行列式函数，u为函数自变量，ω_p为机器人第p阶固有频率；

求解式(11)可得机器人各阶固有频率；

通过模态试验进行第i个网格下阻尼矩阵C_i参数辨识，解耦后获取模态阻尼比，根据式(12)计算机器人末端的振动位移响应：

其中，ζ_p为系统第p阶模态阻尼比，q^p(t)、分别为第p阶广义位移、广义速度和广义加速度，V^p为第p阶增广特征矢量，M_p＝<MV^p,V^p>为系统第p阶模态质量，M为系统质量矩阵，f为外力列阵。

将铣削过程切削力模型作为外部激励输入，基于模态叠加原理将各阶模态响应叠加，计算物理空间振动位移输出，对输入力、输出位移分别进行傅里叶变换，求解刀尖点频响函数G(iω)与颤振频率ω_d，频响函数矩阵为：

其中，G_xx(iω)、G_yy(iω)分别为x方向上的频响函数、y方向上的频响函数，i为虚数单位。

2、进行工作网格划分，绘制变姿态稳定性叶瓣图；

如图3所示为机器人工作空间网格划分示意图，图3中i为网格序号，K_i为第i个网格的关节刚度矩阵，C_i为第i个网格的阻尼矩阵，将机器人铣削工作空间标准化为一个边长为2000mm的大立方体，进一步将整个空间划分成125个边长为400mm的小立方体网格，网格编号为i＝1,2,3…125。设置初始网格边长为h＝400，尺寸变化步长为50，设置网格划分的次数为r，其中，r＝0,1,2…，r＝0表示未划分网格时的初始状态。

控制机器人末端刀具中心点依次处于第1～125个待采样网格，通过计算运动学逆解获取该网格对应机器人姿态下的关节转角，同时对机器人进行模态测试，辨识该网格内机器人关节刚度和阻尼。将转角、刚度、阻尼信息输入动力学仿真系统，计算系统固有频率与刀尖点频响函数。

当系统发生颤振时，在颤振频率ω_d处的特征方程为：

其中，I₂为2行2列的单位矩阵，t为时间变量，A₀为随时间变化的总切削力方向系数矩阵的傅里叶级数展开的第1项；

特征方程的特征值Λ为：

其中，T为切削周期，Λ_R与Λ_I分别为系统特征值的实部和虚部；

特征值满足如下关系式：

b₁Λ²+b₂Λ+1＝0 (16)

式中，b₁、b₂是中间参数，

b₁＝G_xx(iω_d)G_yy(iω_d)(α_xxα_yy-α_xyα_yx)，b₂＝α_xxG_xx(iω_d)+α_yyG_yy(iω_d)，

其中，α_xx、α_yy、α_xy、α_yx为铣削方向系数；

改写式(15)并求解式(16)，求得二自由度铣削系统的临界铣削深度a_lim，对于实际的铣削加工，刀具的轴向切削深度为实数，因此得到：

式中，κ为中间参数，

根据主轴转速n与周期T之间的关系，获得叶瓣图中与式(17)计算的轴向切深对应的转速n为：

式中，k＝1,2,3...10，为铣削稳定性叶瓣图的叶瓣数；

根据实际过程中所需转速范围选择叶瓣数，以主轴转速为横轴，临界轴向铣削深度为纵轴，绘制稳定性叶瓣图。

绘制第1～125个网格的稳定性叶瓣图，统计并记录所有网格稳定性叶瓣图临界切深的平均值CA_rmm和最大值CM_rmm，设置CA_r和CM_r的初值为0.01。

减小网格尺寸至350mm，重新划分网格并重复稳定叶瓣图绘制步骤。计算第1次和第2次临界切深平均值变化率x＝|CA_r-CA_r-1|/CA_r以及最大值变化率y＝|CM_r-CM_r-1|/CM_r，当x、y均小于10％时，输出第2次网格划分下的稳定性叶瓣图；否则继续减小网格尺寸，重复稳定叶瓣图绘制步骤，最终得到变姿态稳定性叶瓣图。

本实例通过上述方法选定网格边长为300mm，通过本发明绘制的KUKA-210铣削机器人姿态Φ₁＝(0,-90,90,0,0,180)°、Φ₂＝(0,-71,90,0,0,180)°、Φ₃＝(-11，-58,95,-96,-95,128)°下的稳定性叶瓣图如图4、图5、图6所示。

具体实现中，本申请提供计算机存储介质以及对应的数据处理单元，其中，该计算机存储介质能够存储计算机程序，所述计算机程序通过数据处理单元执行时可运行本发明提供的一种铣削机器人变姿态稳定性叶瓣图获取方法的发明内容以及各实施例中的部分或全部步骤。所述的存储介质可为磁碟、光盘、只读存储记忆体(read-only memory，ROM)或随机存储记忆体(random access memory，RAM)等。

本领域的技术人员可以清楚地了解到本发明实施例中的技术方案可借助计算机程序以及其对应的通用硬件平台的方式来实现。基于这样的理解，本发明实施例中的技术方案本质上或者说对现有技术做出贡献的部分可以以计算机程序即软件产品的形式体现出来，该计算机程序软件产品可以存储在存储介质中，包括若干指令用以使得一台包含数据处理单元的设备(可以是个人计算机，服务器，单片机。MUU或者网络设备等)执行本发明各个实施例或者实施例的某些部分所述的方法。

本发明提供了一种铣削机器人变姿态稳定性叶瓣图获取方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (7)

1.一种铣削机器人变姿态稳定性叶瓣图获取方法，其特征在于，包括以下步骤：

步骤1，基于多体系统传递矩阵法建立机器人动力学模型，同时建立铣削过程切削力模型，形成铣削机器人动力学仿真系统；

步骤2，将机器人铣削工作空间标准化为一个大立方体，设定大立方体体积为P，进一步将空间划分成Z个边长为h的小立方体网格，Z＝P/h³，Z取整数，网格编号为i，i＝1,2,3…Z；

步骤3，设置初始网格边长为h和尺寸变化步长，设置网格划分的次数为r，其中，r＝0,1,2…，r＝0表示未划分网格时的初始状态；

步骤4，控制机器人末端刀具中心点处于第i个待采样网格，通过计算机器人运动学逆解获取待采样网格对应姿态下的关节转角，对机器人进行模态测试，辨识待采样网格内机器人关节刚度矩阵K_i和阻尼矩阵C_i；

步骤5，将关节转角、关节刚度矩阵K_i和阻尼矩阵C_i输入铣削机器人动力学仿真系统，计算系统固有频率与刀尖点频响函数；

步骤6，基于零阶频域法建立颤振稳定域预测模型，绘制第i个待采样网格的稳定性叶瓣图，记录绝对稳定区的临界切深；

步骤7，统计并记录所有网格稳定性叶瓣图临界切深的平均值CA_r和最大值CM_r，当r＝0时，有CA_r＝0.01和CM_r＝0.01；

步骤8，将网格尺寸减小X₁，重新划分网格，重复步骤4～步骤7；

步骤9，当r>0时，计算第r次和第r-1次临界切深平均值变化率x＝|CA_r-CA_r-1|/CA_r以及最大值变化率y＝|CM_r-CM_r-1|/CM_r，当x、y均小于X₂时，输出当前网格划分下的稳定性叶瓣图；否则重复步骤8，直至满足x、y均小于X₂的条件时，输出变姿态稳定性叶瓣图。

2.根据权利要求1所述的方法，其特征在于，步骤1包括：

步骤1-1，建立铣削过程切削力模型；

步骤1-2，测量各元件动力学参数，拼接各元件传递矩阵，推导机器人动力学模型传递方程，代入边界条件求解特征方程，获得固有振动属性；

步骤1-3，根据解耦后的机器人动力学模态方程，将铣削过程切削力模型输入，基于模态叠加原理获取物理空间内振动位移，完成机器人动力学模型建立步骤；

步骤1-4，将铣削过程切削力模型输入机器人动力学模型形成铣削机器人动力学仿真系统，对铣削力和振动位移进行傅里叶变换，计算刀尖点频响函数。

3.根据权利要求2所述的方法，其特征在于，步骤1-1包括：机器人二自由度的瞬时铣削过程切削力模型表示为：

其中，dF_t、dF_r分别表示切向铣削力微元和径向铣削力微元，K_tc、K_rc分别表示切向力系数和径向力系数，a_c为铣削厚度，da为轴向铣削微元，K_te、K_re分别表示切向刃口力系数和径向刃口力系数，ds为切削刃长度微元；

将径向力F_r、切向力F_t投影到工件坐标系，得到瞬时铣削力：

其中，F_x、F_y分别表示工件坐标系下沿x方向的铣削力、沿y方向的铣削力，φ为铣削刀具接触角，φ_ex、φ_st分别表示铣削切出角和切入角，N为铣削刀具刃数，a为轴向切削深度，f为每齿进给量；

对式(2)积分获得x方向平均铣削力y方向平均铣削力/>

通过铣削正交试验获得不同铣削参数下的平均铣削力，结合式(3)辨识K_tc、K_rc。

4.根据权利要求3所述的方法，其特征在于，步骤1-2包括：机器人各连杆通过关节处齿轮连接，末端法兰携带能够铣削的电主轴及刀具；将机器人整体划分为体元件和铰元件，编号为0～14，其中，2、4、6、8、10、12号为机器人本体连杆，建模为刚体元件；14号铣削刀具刚性较差，建模为空间振动Euler-Bernoulli梁元件；各体元件具有独立的随元件移动的连体坐标系；1、3、5、7、9、11、13号为机器人关节处添加的沿连体坐标系x、y、z轴向具有线弹性，绕x、y、z轴具有扭转弹性的6自由度弹性铰元件，铰元件坐标系x、y、z轴方向与机器人基坐标系方向保持一致；0号元件为机器人基座，固定于地面，视为机器人动力学模型的固定边界；14号元件为末端刀具，视为系统的自由边界；整个机器人动力学模型为树形结构，其传递方向由根0到梢14，模型总传递方程为：

Z_14,0＝U_1-14Z_0,1 (4)

其中，为体元件状态矢量，Z_0,1、Z_14,0分别表示与大地相连元件的状态矢量、自由边界相连元件的状态矢量，U_1-14为总传递矩阵，由各个元件传递矩阵依次拼接得到，X、Y、Z分别表示在体元件坐标系下x、y、z方向模态空间下的线位移，Θ_x、Θ_y、Θ_z分别表示在体元件坐标系下x、y、z方向模态空间下的角位移，M_x、M_y、M_z分别表示在体元件坐标系下x、y、z方向模态空间下的内力矩，Q_x、Q_y、Q_z分别表示在体元件坐标系下x、y、z方向模态空间下对应的力；下标b、j分别为体元件的编号和铰元件的编号，固定边界和自由边界的编号为0；

机器人各元件传递矩阵拼接得到机器人动力学模型总传递矩阵：

U_1-14＝U₁₄U₁₃U₁₂U₁₁U₁₀U₉U₈U₇U₆U₅U₄U₃U₂U₁ (5)

其中，U为模型中各组成元件的传递矩阵，U₁₄为第14个元件的传递矩阵；

根据多体系统传递矩阵法链式推导法则，1号元件与固定边界相连，1号元件的线位移、角位移为0；14号元件与自由边界相连，14号元件的内力、内力矩为0，有如下边界条件：

对于b＝2、4、6、8、10、12号刚体元件，C_b为编号为b连杆的质心点，刚体元件传递矩阵U_b在以输入点I_b为原点，O_b为输出点的自身连体系中描述为：

式中，

其中，I₃为3行3列的单位矩阵，O_3×3为3行3列的零矩阵，m_b为第b个体元件质量，为第b个体元件坐标系中输出点O_b相对于输入点I_b的坐标反对称矩阵，/>为第b个体元件坐标系中质心点C_b相对于输出点O_b的坐标反对称矩阵，/>为第b个体元件坐标系中质心点C_b相对于输入点I_b的坐标反对称矩阵，ω_p为系统第p阶固有频率，J_b,I为刚体相对于输入点I_b的惯量矩阵，J_b,x、J_b,y、J_b,z分别为以输入点I_b为原点的元件坐标系中x、y、z方向的转动惯量，J_b,xy、J_b,xz、J_b,yz分别为空间中x、y、z方向转动惯性积，(g_b1,g_b2,g_b3)为输出点坐标，(c_b1,c_b2,c_b3)为质心坐标；

对于b＝14号梁元件，其传递矩阵为：

式中，

u_2,2＝u_6,6＝u_9,9＝u_11,11＝S(λ_yx_b)，u_8,5＝-EI_yλ_zV(λ_zx_b)，u_9,6＝EI_zλ_yV(λ_yx_b)，

u_3,3＝u_5,5＝u_8,8＝u_12,12＝S(λ_zx_b)，

其中，中间变量λ_y、γ_θx、β_x、λ_z的计算公式为：

其中，中间函数变量S(u)、T(u)、U(u)、V(u)的计算公式为：

其中u为函数自变量；

其中，传递矩阵U₁₄中参数在自身连体坐标系I₁₄中描述，下标b表示元件序号，为线质量密度，x_b表示梁的长度，ρ为梁的密度，E为弹性模量，G为剪切模量，A_b为横截面积，I_x为横截面对x轴的转动惯量，I_y为横截面对y轴的转动惯量，I_z为横截面对z轴的转动惯量，J_p为横截面对x轴的极惯性矩；

对于j＝1、3、5、7、9、11、13号铰元件，其传递矩阵U_j为：

式中，

其中，K″_j表示编号为j的铰元件的线刚度系数矩阵，K′_j表示编号为j的铰元件的扭转刚度系数矩阵，K″_j,x、K″_j,y、K″_j,z分别表示编号为j的铰元件x、y、z方向的线弹性刚度，K′_j,x、K′_j,y、K′_j,z分别表示x、y、z方向的扭转弹性刚度，第i个网格内的关节刚度矩阵K_i包含该网格下所有铰元件的线刚度及扭转刚度信息，需要通过试验辨识获取；

将边界条件代入总传递方程进一步推导系统的特征方程：

det[U_1-14(ω_p)]＝0 (11)

其中，det[u]表示u的行列式函数，u为函数自变量，ω_p为机器人第p阶固有频率；求解式(11)得到机器人各阶固有频率。

5.根据权利要求4所述的方法，其特征在于，步骤1-3包括：通过模态试验进行第i个网格下阻尼矩阵C_i参数辨识，解耦后获取模态阻尼比，根据式(12)计算机器人末端的振动位移响应：

其中，ζ_p为系统第p阶模态阻尼比，q^p(t)、分别为第p阶广义位移、广义速度和广义加速度，V^p为第p阶增广特征矢量，M_p＝<MV^p,V^p>为系统第p阶模态质量，M为系统质量矩阵，f为外力列阵。

6.根据权利要求5所述的方法，其特征在于，步骤1-4包括：将铣削过程切削力力模型作为外部激励输入，基于模态叠加原理将各阶模态响应叠加，计算物理空间振动位移输出，对输入力、输出位移分别进行傅里叶变换，求解刀尖点频响函数G(iω)与颤振频率ω_d，频响函数矩阵为：

其中，G_xx(iω)、G_yy(iω)分别为x方向上的频响函数、y方向上的频响函数，i为虚数单位。

7.根据权利要求6所述的方法，其特征在于，步骤5包括：当系统发生颤振时，在颤振频率ω_d处的特征方程为：

其中，I₂为2行2列的单位矩阵，t为时间变量，A₀为随时间变化的总切削力方向系数矩阵的傅里叶级数展开的第1项；

特征方程的特征值Λ为：

其中，T为切削周期，Λ_R与Λ_I分别为系统特征值的实部和虚部；

特征值满足如下关系式：

b₁Λ²+b₂Λ+1＝0 (16)

式中，b₁、b₂是中间参数，

b₁＝G_xx(iω_d)G_yy(iω_d)(α_xxα_yy-α_xyα_yx)，b₂＝α_xxG_xx(iω_d)+α_yyG_yy(iω_d)，

其中，α_xx、α_yy、α_xy、α_yx为铣削方向系数；

改写式(15)并求解式(16)，求得二自由度铣削系统的临界铣削深度a_lim，对于实际的铣削加工，刀具的轴向切削深度为实数，因此得到：

式中，κ为中间参数，

根据主轴转速n与周期T之间的关系，获得叶瓣图中与式(17)计算的轴向切深对应的转速n为：

式中，k＝1,2,3...10，为铣削稳定性叶瓣图的叶瓣数；

根据实际过程中所需转速范围选择叶瓣数，以主轴转速为横轴，临界轴向铣削深度为纵轴，绘制稳定性叶瓣图。