薄壁件动态铣削稳定性叶瓣图高精度预测方法与流程

薄壁件动态铣削稳定性叶瓣图高精度预测方法属于动态铣削领域。

背景技术：

数控加工技术是现代制造技术的基础，它的广泛应用使普通机械被数控机械所代替，这使制造业发生了根本的变化，在数控加工中，以加工系统稳定、被加工表面的高精度与高材料去除率为主要特征的高性能加工广泛应用于飞机的大梁、隔板、火箭的整流罩、武器（导弹）的战斗部壳体等大材料去除量的关键薄壁结构和零件的制造中。在动态切削过程中存在刀具切削刃与工件之间的相对振动，这种相对振动的主要表现形式之一为自激颤振，若这种自激颤振控制不当，极易诱发切削系统失稳，进而损坏加工零件，甚至破坏加工设备。 所以对自激颤振进行研究以及对极限切削条件（稳定性叶瓣图）进行预测进而合理的选择加工工艺参数来指导实际加工显得十分必要。张俊、赵万华、卢秉恒等人发明的专利一种高速铣削稳定性快速判定方法，即： CN101905340A；利用耦合得到的机床整机频响函数特性，绘制出主轴转数与轴向切深的稳定性极限图，最后为切削进行稳定性判定。该发明大大提高了判断速度，适合于在线的切削参数优化。然而，该发明也存在缺陷，即：该发明通过频域法计算切削稳定性，造成时域/频域转化过程中的精度丢失，使得该发明损失计算精度而追求计算速度。实际上，如果不进行时域/频域转化，则不存在时/频域转化精度丢失误差，进而可以提高预测精度。

技术实现要素：

本发明要解决的技术难题是克服现有技术的缺陷，采用数值计算方法预测圆周螺旋立铣刀在薄壁件铣削过程中的稳定性叶瓣图，解决传统发明中存在的利用频域描述预测稳定性叶瓣图时/频域转化过程中精度丢失的缺陷，提高动态铣削系统稳定性叶瓣图预测精度。

本发明采用的技术方案是薄壁件动态铣削系统稳定性叶瓣图高精度的数值预测方法。首先，考虑薄壁件动态铣削过程，将该动态过程的动力学模型简化为一单自由度的二阶时滞微分方程。然后，采用降阶原理，将该单自由度的二阶时滞微分方程降阶为一一阶时滞微分方程组。再次，在机床主轴旋转的一个周期内对该一阶时滞微分方程组进行离散，对于每个离散的小区间，将该一阶时滞微分方程组分解为两个部分（Part I 与 Part II）。分别采用插值法与积分法对一阶时滞微分方程组中的Part I 与 Part II进行简化，得到该方程解的数学表达式，进而得到一个离散区间首尾时间节点处位移值的递推公式。由此，由主轴旋转的一个周期内最后的时间节点的位移与第一个时间节点的位移建立联系，得到传递矩阵。最后，利用该传递矩阵与Floquet原理，即可计算出动态铣削系统稳定性叶瓣图。与传统的发明相比，利用本文发明的薄壁件动态铣削系统稳定性叶瓣图的数值预测方法精度更高，根据该预测的稳定性叶瓣图选择合适的工艺参数对薄壁件进行加工，从而能够稳定、高效、高精度的加工薄壁件。具体步骤如下：

1. 利用现有的商业设备获取薄壁件在加工过程中沿着刚度最弱方向的模态参数，包含自然频率ω、模态质量M、模态刚度K和阻尼比ξ以及沿着铣刀切向与径向静态切削力系数K_t与K_r，x向指的是铣刀的进给方向，y向指的是薄壁件刚度最小的方向，即薄壁方向；

2. 构建动态铣削薄壁件过程的数学模型——可以简化为一单自由度的二阶时滞微分方程，该方程可以用下式表示：

其中，M 为薄壁件铣削区域沿着y方向的模态质量，C为薄壁件铣削区域沿着y方向的阻尼，K为薄壁件铣削区域沿着y方向的模态刚度，a为薄壁件铣削区域铣刀轴向切削深度，K_c(t)为薄壁件铣削过程中沿着y方向的动态切削力系数，a为轴向切削深度，可以有下式表示：

N为刀具所含切削刃个数，φ_j为第j个切削刃t时刻的位置角，K_t和K_r分别为铣削系统沿着铣刀切向与径向静态切削力系数，g(φ_j)函数为一单位跃阶函数，用于确定切削刃是否参与切削，为1时表示此刻第j个切削刃正在参与切削，否则表示第j个切削刃此刻没有参与切削，可以由下式表示：

上式中，φ_st与φ_ex分别表示切削刃的切入切出角，y(t)为薄壁件铣削系统的振动位移；

3. 根据建立的薄壁件动态铣削过程的动力学模型，采用降阶原理，将该动力学模型的二阶时滞微分方程降阶为一一阶时滞微分方程组，该一阶时滞微分方程组可以用下式表示：

具体采用的降阶公式如下：

这时，薄壁件动态铣削数学模型可以由一阶时滞微分方程组表示；

4. 对铣刀旋转一个时间周期进行均匀离散，由于铣刀安装在主轴上，所以铣刀旋转一个周期实际上为主轴旋转一个周期，可以用T表示，具体由下式计算得到：

同时，设离散区间数为m，则每个区间长度可以用h=T/m表示，每个时间节点分别为t₀、t₁、t₂、…、t_m；

5. 在每个小的时间区间内对一阶时滞微分方程组进行求解，例如在第i个时间区间内，该一阶时滞微分方程组的解可用下式表示；

为了书写方便，可以把x(t_i+1)简化为x_i+1，ε∈（t_i t_i+1]为积分参数。于是，上式可以简化为：

；

6. 对一阶时滞微分方程组的解进行分部分（第一部分： Part　I　与　第二部分：　Part II）处理，具体将上式的积分部分分解为两个部分，分别为：

；

7. 对第6步得到的第一部分与第二部分分别采用插值法与积分法进行离散，这里对第一部分采用三阶插值算法，第二部分采用积分算法。对于第一部分：

上式中，ζ∈（0 h], a、b、c和d分别是系数，可以由下式表示：

所以，将上述两式代入第一部分的公式可以得到：

对于第二部分：

由于T是循环周期，所以可以得出：

所以第二部分可以进一步简化：

；

8. 对一阶时滞微分方程组的解计算公式进行整合，对第一部分以及第二部分的公式进行整合并带回到一阶时滞微分方程组的解，得到：

合并同类项得到：

；

9. 对第8步方程进行整合，得到整合方程实际上为一递推公式，对该公式进行变换，得到矩阵递推公式：

；

10. 利用第9步得到的递推矩阵公式，从而可以推得一个周期内Y₀与Y_m之间的联系，可以由下式表示：

定义D为传递矩阵，计算薄壁件在铣削过程中的稳定性极限图时求解传递矩阵D的特征值，根据Floquet原理，若全部特征值绝对值的最大值小于1，则该系统稳定。同时，在改特征矩阵中，未知数为主轴转速n与轴向切削深度a。给定主轴转速n，并给定传递矩阵D特征值绝对值最大值为1，由此可得对应的最大轴向切削深度a。 进而获取了薄壁件在铣削过程中的切削稳定性叶瓣图。

本发明具有以下明显效果：能够有效的提高薄壁件铣削加工叶瓣图的预测精度。

附图说明

图1，某薄壁零件的局部视图；

图2，薄壁零件动态铣削示意图，其中k_y表示薄壁件沿着y轴方向的模态刚度，c_y为薄壁件沿着y轴方向的模态阻尼，φ_j为第j个切削刃t时刻的位置角，j、j-1表示铣刀的第j个与第j-1个切削刃，“第j个切削刃路径”与“第j-1个切削刃路径”表示第j个与第j-1个切削刃切削过程中运动轨迹；

图3，薄壁件铣削过程中的坐标系定义，其中O表达坐标系的原点，x、y、z分别表示该坐标系的x、y与z轴；

图4，主轴旋转一个周期的离散示意图，T表示主轴旋转的周期，m代表对主轴旋转一个周期进行离散的区间数，t₀、t₁、t₂、... ... t_m-2、t_m-1、 t_m分别表示m个离散区间的时间节点，x(t₀)、x(t₁)、x(t₂)、... ... 、x(t_m-2)、x(t_m-1)、x(t_m)表示对应于时间节点处当前切削刃(设为第j个)的位移值，h表示每个离散小区间的时间长度；

图5，对时滞微分方程组的时间滞后项简化示意图，其中m为对主轴旋转一个周期进行离散的区间数，在该示意图中取m=7，x(t₅-T)表示第5个离散区间右端点对应的时间节点(第6个时间节点)处相对于当前切削刃的前一切削刃(设为第j-1个)的时间滞后位移值，x(t_-2)表示相对于当前周期的前一周期倒数第2个时间节点处的当前切削刃的前一切削刃(设为第j-1个)的时间滞后位移值；

图6，利用传统方法预测的薄壁件铣削稳定性叶瓣图，其中由圆圈组成的非连续曲线为实际的薄壁件铣削过程中稳定性叶瓣图，连续曲线为利用传统方法计算的薄壁件铣削过程中的稳定性叶瓣图；

图7，利用本发明提出的方法预测薄壁件铣削过程稳定性极限图, 其中由圆圈组成的非连续曲线为实际的薄壁件铣削过程中稳定性叶瓣图，连续曲线为利用本发明提出的方法计算的薄壁件铣削过程中的稳定性叶瓣图。

具体实施方式

下面结合附图和技术方案，详细说明本发明的具体实施。薄壁件（如图１所示）在数控加工中，常需要优化进给率、定制加工参数等，从而达到加工系统稳定，保证高精度的同时提高铣削效率，这一切的必要条件为高精度的预测薄壁件在加工过程中的稳定性叶瓣图。由于薄壁件铣削动态铣削系统（如图２所示）的动力学模型可以简化为一二阶时间滞后微分方程，不能用解析方法进行预测，为了提高薄壁件铣削系统稳定性叶瓣图的预测精度，本发明提出一种薄壁件动态铣削稳定性叶瓣图高精度预测方法。首先，考虑薄壁件动态铣削过程，将该动态过程的动力学模型简化为一单自由度的二阶时滞微分方程。然后，采用降阶原理，将该单自由度的二阶时滞微分方程降阶为一一阶时滞微分方程组。再次，在机床主轴旋转的一个周期内对该一阶时滞微分方程组进行离散，对于每个离散的小区间，将该一阶时滞微分方程组分解为两个部分（Part I 与 Part II）。分别采用插值法与积分法对一阶时滞微分方程组中的Part I 与 Part II进行简化，得到该方程解的数学表达式，进而得到一个离散区间首尾时间节点处位移值的递推公式。由此，由主轴旋转的一个周期内最后的时间节点的解与第一个时间节点的解建立联系，得到传递矩阵。最后，利用该传递矩阵与Floquet原理，即可计算出动态铣削系统稳定性叶瓣图。该模型预测薄壁件在铣削过程中的稳定性叶瓣图的精度更高；方法的具体步骤如下：

1. 利用现有的商业设备获取薄壁件在加工过程中沿着刚度最弱方向的模态参数，包含自然频率ω、模态质量M、模态刚度K和阻尼比ξ以及沿着铣刀切向与径向静态切削力系数K_t与K_r，x向与y向如图3所示，图3所示的为薄壁件铣削过程中的坐标系定义，其中O表达坐标系的原点，x、y、z分别表示该坐标系的x、y与z轴；

2. 构建动态铣削薄壁件过程的数学模型，如图2所示，图2所示的为薄壁零件动态铣削示意图，其中k_y表示薄壁件沿着y轴方向的模态刚度，c_y为薄壁件沿着y轴方向的模态阻尼，φ_j为第j个切削刃t时刻的位置角，j、j-1表示铣刀的第j个与第j-1个切削刃，“第j个切削刃路径”与“第j-1个切削刃路径”表示第j个与第j-1个切削刃切削过程中运动轨迹；此时，动态铣削薄壁件过程的数学模型可以简化为一单自由度的二阶时滞微分方程，该方程可以用下式表示：

其中，M 为薄壁件铣削区域沿着y方向的模态质量，C为薄壁件铣削区域沿着y方向的阻尼，K为薄壁件铣削区域沿着y方向的模态刚度，a为薄壁件铣削区域铣刀轴向切削深度，K_c(t)为薄壁件铣削过程中沿着y方向的动态切削力系数，a为轴向切削深度，可以有下式表示：

N为刀具所含切削刃个数，φ_j为第j个切削刃t时刻的位置角（如图3所示），K_t和K_r分别为铣削系统沿着铣刀切向与径向静态切削力系数，g(φ_j)函数为一单位跃阶函数，用于确定切削刃是否参与切削，为1时表示此刻第j个切削刃正在参与切削，否则表示第j个切削刃此刻没有参与切削，可以由下式表示：

上式中，φ_st与φ_ex分别表示切削刃的切入切出角，y(t)为薄壁件铣削系统的振动位移；

3. 根据建立的薄壁件动态铣削过程的动力学模型，采用降阶原理，将该动力学模型的二阶时滞微分方程降阶为一一阶时滞微分方程组，该一阶时滞微分方程组可以用下式表示：

具体采用的降阶公式如下：

这时，薄壁件动态铣削数学模型可以由一阶时滞微分方程组表示；

4. 对铣刀旋转一个时间周期进行均匀离散，如图4所示，图4所示的为主轴旋转一个周期的离散示意图，T表示主轴旋转的周期，m代表对主轴旋转一个周期进行离散的区间数，t₀、t₁、t₂、... ... t_m-2、t_m-1、 t_m分别表示m个离散区间的时间节点，x(t₀)、x(t₁)、x(t₂)、... ... 、x(t_m-2)、x(t_m-1)、x(t_m)表示对应于时间节点处当前切削刃(设为第j个)的位移值，h表示每个离散小区间的时间长度；由于铣刀安装在主轴上，所以铣刀旋转一个周期实际上为主轴旋转一个周期，可以用T表示，具体由下式计算得到：

同时，设离散区间数为m，则每个区间长度可以用h=T/m表示，每个时间节点分别为t₀、t₁、t₂、…、t_m；

5. 在每个小的时间区间内对一阶时滞微分方程组进行求解，例如在第i个时间区间内，该一阶时滞微分方程组的解可用下式表示；

为了书写方便，可以把x(t_i+1)简化为x_i+1，ε∈（t_i t_i+1]为积分参数。于是，上式可以简化为：

6. 对一阶时滞微分方程组的解进行分部分（第一部分： Part　I　与　第二部分：　Part II）处理，具体将上式的积分部分分解为两个部分，分别为：

；

7. 对第6步得到的第一部分与第二部分分别采用插值法与积分法进行离散，这里对第一部分采用三阶插值算法，第二部分采用积分算法。 对于第一部分：

上式中，ζ∈（0 h], a、b、c和d分别是系数，可以由下式表示：

所以，将上述两式代入第一部分的公式可以得到：

对于第二部分：

由于T是循环周期，如图5所示，图5所示的为对时滞微分方程组的时间滞后项简化示意图，其中m为对主轴旋转一个周期进行离散的区间数，在该示意图中取m=7，x(t₅-T)表示第5个离散区间右端点对应的时间节点(第6个时间节点)处相对于当前切削刃的前一切削刃(设为第j-1个)的时间滞后位移值，x(t_-2)表示相对于当前周期的前一周期倒数第2个时间节点处（当前切削刃的前一切削刃为第j-1个）的时间滞后位移值；所以可以得出：

所以第二部分可以进一步简化：

；

8. 对一阶时滞微分方程组的解计算公式进行整合，对第一部分以及第二部分的公式进行整合并带回到一阶时滞微分方程组的解，得到：

合并同类项得到：

；

9. 对第8步方程进行整合，得到整合方程实际上为一递推公式，对该公式进行变换，得到矩阵递推公式：

；

10. 利用第9步得到的递推矩阵公式，从而可以推得一个周期内Y₀与Y_m之间的联系，可以由下式表示：

定义D为传递矩阵，计算薄壁件在铣削过程中的稳定性极限图时求解传递矩阵D的特征值，根据Floquet原理，若全部特征值绝对值的最大值小于1，则该系统稳定。同时，在改特征矩阵中，未知数为主轴转速n与轴向切削深度a。给定主轴转速n，并给定传递矩阵D特征值绝对值的最大值为1，由此可得对应的最大轴向切削深度a。 进而获取了薄壁件在铣削过程中的切削稳定性叶瓣图。

实施例：选用图6，图7所示的实施例来说明本发明的特点：

在图6与图7所示的实施例中，选择的加工参数与刀具参数完全相同，铣刀所含切削刃的个数为N=2，铣刀直径为12mm，径向切削深度为12mm，设离散区间数为m=200的情况下的薄壁件稳定性叶瓣图为真实值。

薄壁件的模态参数与铣削过程中的铣削力系数如分别为:

由此可以算出模态刚度的值:

阻尼为:

分别利用传统方法与本发明提出的方法(m=40)进行计算薄壁件在切削过程中的稳定性叶瓣图，在计算过程中主轴转速范围与铣刀轴向切削深度范围分别为:

经过计算得出两种方法预测的稳定性叶瓣图，图6为传统方法预测的稳定性叶瓣图，其中不连续曲线为真实的叶瓣图曲线，连续曲线为利用传统方法预测的叶瓣图曲线，发现二者不能很好的贴合。

图7为利用本发明提出的预测方法预测的在铣削薄壁件过程中的稳定性叶瓣图与真实稳定性叶瓣图的比较，其中不连续曲线为真实的叶瓣图曲线，连续曲线为利用本发明提出的方法预测的叶瓣图曲线，可以发现二者贴合非常紧密。

本实施例说明本发明提出的方法预测的稳定性叶瓣图精度更高，说明了本发明专利的必要性。