CN111611725A - 一种基于Cotes数值积分的铣削稳定域预测方法 - Google Patents

Abstract

本发明公开了一种基于Cotes数值积分的铣削稳定域确定方法，其通过建立考虑再生效应的状态空间形式的铣削动力学微分方程：将一个周期内的连续时间t表示为离散的时间节点：计算在区间内的方程表达式：计算关于状态项的等式，得到系统在一个周期内的状态转移矩阵，最终以主轴转速为横坐标，轴向切深为纵坐标进行绘图即可得到稳定性图。选择铣削加工参数时在黑色曲线下方稳定区域，选择对应的主轴转速和轴向切深进行铣削加工，即获得无颤振切削情况。本发明减少了指数矩阵和状态转移矩阵的计算次数，使得计算效率与离散类方法相比有很大提高；并且由于数学近似误差达到O(h ⁶),计算精度与频域法以及离散法相比也有较大提高。

Description

一种基于Cotes数值积分的铣削稳定域预测方法

技术领域

本发明属于数控铣削加工技术，应用于预测数控铣削加工过程中的颤振问题。特别涉及一种基于Cotes数值积分的铣削稳定域预测方法。

背景技术

随着航空、航天、船舶及车辆等现代制造业的飞速发展和科技的不断进步。对于产品性能的要求也不断提升，一些整体尺寸大，结构形状复杂，材料去除率高的结构件，如：飞机大梁、壁板以及发动机整体叶盘等结构件的加工，在满足设计性能要求的同时给制造工艺和加工效率带来了很大难度。如何对这些关键零件进行高质量、高效率的制造，是我国制造业面临的严峻挑战。在众多的材料成型方式中，高速铣削加工仍然是进行整体结构件加工制造的关键技术和应用最广泛的加工方式。由于这些结构件具有结构复杂、弱刚性和材料难加工等特性，在铣削时极易发生切削颤振，严重影响零件的表面质量以及刀具的使用寿命，影响加工制造的生产效率，使高速机床性不能得到充分的发挥，还可能会造成零件的报废。同时在一些航空航天器的设计中，追求轻量化，高效化而使用一些如钛合金等难加工材料更将颤振问题进一步凸显。如何采取有效手段抑制切削过程中的颤振，实现无颤振切削以大幅度提高表面质量和生产效率是现代航空制造业所面临的严峻挑战。

对于铣削颤振预测问题，主要是预报由主轴转速和临界轴向切深构成的铣削过程稳定性图，针对这一类问题，国内外学者提出了许多近似解决方案。如：国外学者发表在《CIRP.Annals-Manufacturing Technology》上的文章“Analytical prediction ofstability lobes in milling”公开了一种计算稳定域的零阶频域方法，能够快速准确预测一般情况下的铣削稳定区域，但无法预测小径向切深情况下稳定区域。又如：2008年国外学者发表在《Journal of Sound and Vibration》杂志上的“On the higher-order semi-discretization for periodic delayed systems”提出了一种半离散法来预测铣削稳定区域，该方法能获得较高的计算精度，但是该方法的计算效率较低。2010年国内学者发表在《International Journal of Machine Tools&Manufacture》上的“A full-discretization method for prediction of milling stability”文章公开了一种全离散方法确定的铣削稳定区域，能获得和半离散法相同的计算精度而计算效率有所提高。再如：中国专利“一种基于正交多项式的铣削稳定性预测方法”(公告号：104680000)，在计算稳定性叶瓣图时，传递矩阵需要离散间隔数m次迭代，计算效率仍然不够高。一般情况下，在实际的工业应用中，对于铣削系统的动力学分析通常涉及较大的参数矩阵，在进行稳定性预测时，将消耗大量的计算时间。对于现有的频域类方法和离散类方法而言，难以同时满足通用性强，计算精度高，计算效率高等要求。

发明内容

基于以上问题，本发明提供一种基于Cotes数值积分的铣削稳定域预测方法，满足高效率、高精度、铣削稳定域预测效果，避免出现切削颤振现象，获得符合要求的零件表面质量并提高加工效率。

本发明的目的是这样实现的。一种基于Cotes数值积分的铣削稳定域确定方法，其步骤如下：

步骤一：建立考虑再生效应的状态空间形式的铣削动力学微分方程：

具体步骤如下：

S01:首先考虑再生效应的n个自由度的铣削动力学微分方程可以表述为：

其中，M、C和K为刀具系统的n个自由度的模态质量矩阵、模态阻尼矩阵和模态刚度矩阵，q(t)是刀具n个自由度的振动位移矢量，K_c(t)为系统所受动态铣削力矩阵，t为连续时间，T为单个刀齿切削周期，a_p是轴向切削深度；

S02:令

方程(1)转换为状态空间形式：

其中：

步骤二：将一个周期内的连续时间t表示为离散的时间节点：

t_i＝t₀+t_f+(i-1)τ； (3)

式(3)中，i＝1,2,...,m,m+1.t₀为开始切削时刻，t_f自由振动时间段，τ为离散间隔长度,m为单个周期离散数；

步骤三：计算在区间[t_i,t_i+1]内的方程(2)表达式：

其具体步骤如下：

S03:将状态方程的a_pB(t)[x(t)-x(t-T)]看作齐次方程

的非齐次项，s为计算过程中的积分因子，则可以将方程(2)表示为如下形式：

S04:根据式(4)可知，当t_i≤t≤t_i+1时，可以得到在区间[t_i,t_i+1]内的表达式如下：

步骤四：计算关于状态项x₁，x₂和x₃的等式：

具体步骤如下：

S05:当刀具不切削的时刻，即在时间段t_f内，此时B(s)为0，方程(5)退化为

S06:为了公式表示的简洁性，下文统一使用x_i代替x(t_i)，x_i-T代替x(t_i-T)，B_i代替B(t_i)，B_i-T代替B(t_i-T)；

S07:在t₁＝t₀+t_f时刻，由式(3)和(5)可以很容易得到以下关于状态项x₁的等式：

S08:在离散点t₂处，状态项x₂可以表示为：

S09:由梯形求积公式，可以得到状态项x₂的近似表达公式：

S10:移项整理后，分离出状态项和时滞项如下式：

S11:类似地，在离散点t₃处，状态项x₃可以得到：

S12:由Simpson求积公式，可以得到x₃的近似表达式为：

S13:同样，分离状态项和时滞项可得到：

步骤五：计算关于状态项x₄的等式

S14:与S11步骤类似，在离散点t₄处，状态项x₄可以得到：

S15:由Newton求积公式，分离状态项和时滞项可以得到：

步骤六：计算第t₅到t_m+1点的求积公式：

其中：

具体步骤如下：

S16:第t₅到t_m+1点的求积公式可表示为：

其中，i＝1,2,…,m-3；

S17:由Cotes求积公式，分离状态项和时滞项可以得到：

步骤七：得到系统在一个周期内的状态转移矩阵，表示为F_IM＝G^-1H，其中，

具体步骤如下：

S18:联立式(6)、(9)、(12)、(14)和(16)得到如下离散映射：

Gy_m+1＝Hy_m+1-T； (17)

S19:系统在一个周期内的状态转移矩阵表示为：

F_IM＝G^-1H； (18)

步骤八：计算状态转移矩阵F_IM的特征值，通过判断特征值的模的大小来判断稳定性，具体公式为：

最终以主轴转速为横坐标，轴向切深为纵坐标进行绘图即可得到稳定性图；

步骤九：在步骤八绘制的图中，选择铣削加工参数时在曲线下方稳定区域，选择对应的主轴转速和轴向切深进行铣削加工，即获得无颤振切削情况，切削后得到符合要求的零件表面；在不稳定区域选择主轴转速和轴向切深进行铣削加工，则会发生切削颤振，切削后得到不合要求的零件表面。

本发明减少了指数矩阵和状态转移矩阵的计算次数，使得计算效率与离散类方法相比有很大提高；并且由于数学近似误差达到O(h⁶),计算精度与频域法以及离散法相比也有较大提高。

附图说明

图1是本发明计算得到的单自由度铣削稳定性区域分布图；

图2是无颤振切削后的零件表面图；

图3是发生颤振时的零件表面图。

具体实施方式

以下结合附图和实施例对本发明做进一步说明。参见图1至图3，一种基于Cotes数值积分的铣削稳定域确定方法，具体步骤如下：

步骤一：建立考虑再生效应的状态空间形式的铣削动力学微分方程：

具体步骤如下：

S01:首先考虑再生效应的n个自由度的铣削动力学微分方程可以表述为：

其中，M、C和K为刀具系统的n个自由度的模态质量矩阵、模态阻尼矩阵和模态刚度矩阵，q(t)是刀具n个自由度的振动位移矢量，K_c(t)为系统所受动态铣削力矩阵，t为连续时间，T为单个刀齿切削周期，a_p是轴向切削深度；

S02:令

方程(1)转换为状态空间形式：

其中：

步骤二：将一个周期内的连续时间t表示为离散的时间节点：

t_i＝t₀+t_f+(i-1)τ； (3)

式(3)中，i＝1,2,...,m,m+1.t₀为开始切削时刻，t_f自由振动时间段，τ为离散间隔长度,m为单个周期离散数；

步骤三：计算在区间[t_i,t_i+1]内的方程(2)表达式：

其具体步骤如下：

S03:将状态方程的a_pB(t)[x(t)-x(t-T)]看作齐次方程

的非齐次项，s为计算过程中的积分因子，则可以将方程(2)表示为如下形式：

S04:根据式(4)可知，当t_i≤t≤t_i+1时，可以得到在区间[t_i,t_i+1]内的表达式如下：

步骤四：计算关于状态项x₁，x₂和x₃的等式：

具体步骤如下：

S05:当刀具不切削的时刻，即在时间段t_f内，此时B(s)为0，方程(5)退化为

S06:为了公式表示的简洁性，下文统一使用x_i代替x(t_i)，x_i-T代替x(t_i-T)，B_i代替B(t_i)，B_i-T代替B(t_i-T)；

S07:在t₁＝t₀+t_f时刻，由式(3)和(5)可以很容易得到以下关于状态项x₁的等式：

S08:在离散点t₂处，状态项x₂可以表示为：

S09:由梯形求积公式，可以得到状态项x₂的近似表达公式：

S10:移项整理后，分离出状态项和时滞项如下式：

S11:类似地，在离散点t₃处，状态项x₃可以得到：

S12:由Simpson求积公式，可以得到x₃的近似表达式为：

S13:同样，分离状态项和时滞项可得到：

步骤五：计算关于状态项x₄的等式

S14:与S11步骤类似，在离散点t₄处，状态项x₄可以得到：

S15:由Newton求积公式，分离状态项和时滞项可以得到：

步骤六：计算第t₅到t_m+1点的求积公式：

其中：

具体步骤如下：

S16:第t₅到t_m+1点的求积公式可表示为：

其中，i＝1,2,…,m-3；

S17:由Cotes求积公式，分离状态项和时滞项可以得到：

步骤七：得到系统在一个周期内的状态转移矩阵，表示为F_IM＝G^-1H，其中，

具体步骤如下：

S18:联立式(6)、(9)、(12)、(14)和(16)得到如下离散映射：

Gy_m+1＝Hy_m+1-T； (17)

S19:系统在一个周期内的状态转移矩阵表示为：

F_IM＝G^-1H； (18)

步骤八：计算状态转移矩阵F_IM的特征值，通过判断特征值的模的大小来判断稳定性，具体公式为：

最终以主轴转速为横坐标，轴向切深为纵坐标进行绘图；

步骤九：设定相应的切削参数，铣刀刀齿数目为2，铣削方式为顺铣，径向浸入比设置为0.1，切向切削力系数和径向切削力系数分别为6×10⁸和2×10⁸，铣刀的一阶固有频率为922×2×π，阻尼比为0.011，模态质量为0.03993。单个刀齿的切削周期离散数为40，将由主轴转速与轴向切深构成的平面划分为400×200网格。

步骤十：将步骤一到步骤八的全部过程使用Matlab软件编写成程序，将步骤九给定的切削参数输入程序进行计算，得到如附图1所示的稳定性图。

步骤十一：在图1中，选择铣削加工参数时在黑色曲线下方稳定区域(空白部分)选择对应的主轴转速和轴向切深进行铣削加工，即可获得无颤振切削情况，切削后得到符合要求的零件表面如附图2所示；在黑色曲线上方不稳定区域(阴影部分)选择主轴转速和轴向切深进行铣削加工，则会发生切削颤振，切削后得到不合要求的零件表面(如图3所示)。

以上所述实施方式仅为本发明的优选实施例，而并非本发明可行实施的穷举。对于本领域一般技术人员而言，在不背离本发明原理和精神的前提下对其所作出的任何显而易见的改动，都应当被认为包含在本发明的权利要求保护范围之内。

Claims (1)

1.一种基于Cotes数值积分的铣削稳定域确定方法，其特征在于，其步骤如下：

步骤一：建立考虑再生效应的状态空间形式的铣削动力学微分方程：

具体步骤如下：

S01:首先考虑再生效应的n个自由度的铣削动力学微分方程可以表述为：

其中，M、C和K为刀具系统的n个自由度的模态质量矩阵、模态阻尼矩阵和模态刚度矩阵，q(t)是刀具n个自由度的振动位移矢量，K_c(t)为系统所受动态铣削力矩阵，t为连续时间，T为单个刀齿切削周期，a_p是轴向切削深度；

S02:令

方程(1)转换为状态空间形式：

其中：

步骤二：将一个周期内的连续时间t表示为离散的时间节点：

t_i＝t₀+t_f+(i-1)τ； (3)

式(3)中，i＝1,2,...,m,m+1.t₀为开始切削时刻，t_f自由振动时间段，τ为离散间隔长度,m为单个周期离散数；

步骤三：计算在区间[t_i,t_i+1]内的方程(2)表达式：

其具体步骤如下：

S03:将状态方程的a_pB(t)[x(t)-x(t-T)]看作齐次方程

的非齐次项，s为计算过程中的积分因子，则可以将方程(2)表示为如下形式：

S04:根据式(4)可知，当t_i≤t≤t_i+1时，可以得到在区间[t_i,t_i+1]内的表达式如下：

步骤四：计算关于状态项x₁，x₂和x₃的等式：

具体步骤如下：

S05:当刀具不切削的时刻，即在时间段t_f内，此时B(s)为0，方程(5)退化为

S06:为了公式表示的简洁性，下文统一使用x_i代替x(t_i)，x_i-T代替x(t_i-T)，B_i代替B(t_i)，B_i-T代替B(t_i-T)；

S07:在t₁＝t₀+t_f时刻，由式(3)和(5)可以很容易得到以下关于状态项x₁的等式：

S08:在离散点t₂处，状态项x₂可以表示为：

S09:由梯形求积公式，可以得到状态项x₂的近似表达公式：

S10:移项整理后，分离出状态项和时滞项如下式：

S11:类似地，在离散点t₃处，状态项x₃可以得到：

S12:由Simpson求积公式，可以得到x₃的近似表达式为：

S13:同样，分离状态项和时滞项可得到：

步骤五：计算关于状态项x₄的等式

S14:与S11步骤类似，在离散点t₄处，状态项x₄可以得到：

S15:由Newton求积公式，分离状态项和时滞项可以得到：

步骤六：计算第t₅到t_m+1点的求积公式：

其中：

具体步骤如下：

S16:第t₅到t_m+1点的求积公式可表示为：

其中，i＝1,2,…,m-3；

S17:由Cotes求积公式，分离状态项和时滞项可以得到：

步骤七：得到系统在一个周期内的状态转移矩阵，表示为F_IM＝G^-1H，

其中，

具体步骤如下：

S18:联立式(6)、(9)、(12)、(14)和(16)得到如下离散映射：

Gy_m+1＝Hy_m+1-T； (17)

S19:系统在一个周期内的状态转移矩阵表示为：

F_IM＝G^-1H； (18)

步骤八：计算状态转移矩阵F_IM的特征值，通过判断特征值的模的大小来判断稳定性，具体公式为：

最终以主轴转速为横坐标，轴向切深为纵坐标进行绘图；

步骤九：在步骤八绘制的图中，选择铣削加工参数时在曲线下方稳定区域，选择对应的主轴转速和轴向切深进行铣削加工，即获得无颤振切削情况，切削后得到符合要求的零件表面；在不稳定区域选择主轴转速和轴向切深进行铣削加工，则会发生切削颤振，切削后得到不合要求的零件表面。

Images