一种预测模具型腔数控铣削中颤振的方法，本发明涉及预测模具型腔数控铣削中颤振的方法。本发明的目的是为了解决现有单一刀具路径的铣削稳定性预测方法适用性低，导致铣削颤振预测准确度低，加快刀具失效，影响模具型腔的加工质量的问题。一种预测模具型腔数控铣削中颤振的方法具体过程为：步骤一.建立刀具‑工件系统的相对传递函数；步骤二.将步骤一获得的刀具‑工件系统的相对传递函数引入三维铣削稳定性模型中，得到铣刀铣削颤振频率处的临界轴向切削深度；步骤三、基于步骤二得到的临界轴向切削深度判断模具型腔数控铣削是否发生颤振。本发明用于模具型腔数控铣削领域。
1.一种预测模具型腔数控铣削中颤振的方法，其特征在于：所述方法具体过程为：步骤一.建立刀具-工件系统的相对传递函数；具体过程为：
Gci(jω)＝Xci(jω)/[Fci(jω)] (1)Gwi(jω)＝Xwi(jω)/[Fwi(jω)] (2)式中：Fci(jω)和Fwi(jω)分别为刀具子系统刀尖点和工件子系统工件表面铣削点所受的力，由于刀具和工件之间的作用力大小相等、方向相反，即Fci(jω)＝-Fwi(jω)；Xci(jω)和Xwi(jω)分别为在Fci(jω)和Fwi(jω)作用下产生的位移；i＝1，2，3分别表示遵循右手笛卡尔直角坐标系建立的直角坐标系的X方向、Y方向、Z方向；
由式(1)和(2)，推导出刀具-工件系统的相对传递函数：
步骤二.将步骤一获得的刀具-工件系统的相对传递函数引入三维铣削稳定性模型中，得到铣刀铣削颤振频率处的临界轴向切削深度；
步骤三、基于步骤二得到的临界轴向切削深度判断模具型腔数控铣削是否发生颤振；
所述步骤二中将步骤一获得的刀具-工件系统的相对传递函数引入三维铣削稳定性模型中，得到铣刀铣削颤振频率处的临界轴向切削深度；具体过程为：基于频域法的三维铣削稳定性模型铣刀铣削模具凹槽型腔侧面的过程中，切向力Ft沿铣刀铣削刃铣削速度方向分布，径向力Fr是铣刀进给的径向方向，轴向力Fa沿铣刀轴向作用；
在铣削刃线中，铣削刃被分成M个微分单元，M取值为正整数；
设第j个铣削刃上第k个微分单元的转角量 表达式为：
其中：转角量 为第j个铣削刃上第k个微分单元的转角量，θr(k,t)是铣削刃第k个微分单元在笛卡尔坐标Y轴方向测量的旋转角度，0≤k≤M；Nf表示球头铣刀所含有的刀齿数量，n表示刀具主轴转速，t为时间；j取值为正整数；
其中： 是瞬态切屑厚度，Δa是一个微分单元的长度；
切向力Ft、径向力Fr、轴向力Fa在铣刀铣刃上的微分单元是
dFt＝KtcdAc dFr＝KrcdAc dFa＝KacdAc (5)式中Ktc为切向力的铣削力系数，Krc为径向力的铣削力系数，Kac为轴向力的铣削力系数；
切向力Ft、径向力Fr和轴向力Fa转化为刀具笛卡尔坐标下的三个方向的动态铣削力：式中 为笛卡尔坐标X轴方向的动态铣削力， 为笛卡尔坐标Y轴方向的动态铣削力， 为笛卡尔坐标Z轴方向的动态铣削力，ap为刀具轴向切削深度， 为转角量，Kr为径向力的铣削力系数与切向力的铣削力系数之比，Kr＝Krc/Ktc；Ka为轴向力的铣削力系数与切向力的铣削力系数之比，Ka＝Kac/Ktc；Δxj＝(xc(t)-xc(t-T))-(xw(t)-xw(t-T))，Δyj＝(yc(t)-yc(t-T))-(yw(t)-yw(t-T))，Δzj＝(zc(t)-zc(t-T))-(zw(t)-zw(t-T))，Δxj、Δyj、Δzj均为再生效应振动；
其中xc(t)为铣刀当前刀齿的X方向动态位移，xc(t-T)为铣刀前一个刀齿周期的X方向动态位移，xw(t)为工件在当前刀齿铣削下的X方向动态位移，xw(t-T)为工件在前一个刀齿周期铣削下的X方向动态位移；
yc(t)为铣刀当前刀齿的Y方向动态位移，yc(t-T)为铣刀前一个刀齿周期的Y方向动态位移，yw(t)为工件在当前刀齿铣削下的Y方向动态位移，yw(t-T)为工件在前一个刀齿周期铣削下的Y方向动态位移；
zc(t)为铣刀当前刀齿的Z方向动态位移，zc(t-T)为铣刀前一个刀齿周期的Z方向动态位移，zw(t)为工件在当前刀齿铣削下的Z方向动态位移，zw(t-T)为工件在前一个刀齿周期铣削下的Z方向动态位移；T为周期；
[F(t)]＝apKtc[A(t)]{Δr} (7)其中：F(t)为以 为元素的列向量，{Δr}为以Δxj、Δyj、Δzj为元素的列向量，动态铣削力方向系数矩阵[A(t)]是角频率为ω，周期为T的周期函数；ω＝NfΩ，T＝2π/ω；Ω为刀具主轴转动角速度，单位为rad/s；
将[A(t)]进行傅里叶级数展开，得到平均方向系数矩阵：
式中，Nf表示球头铣刀所含有的刀齿数量， 为刀具刀齿的水平切入角， 为刀具刀齿的水平切出角， 为系数矩阵，为转角量；
由式(8)，将时变动态铣削方程[F(t)]转换为不随时间变化但与刀具-工件接触状态相关的形式：[F(t)]＝apKtc[A(0)]{Δr} (9)其中，
{Δr}＝{rc(t)}-{rc(t-T)}-({rw(t)}-{rw(t-T)})式中，rc(t)为铣刀当前刀齿的X方向、Y方向、Z方向的动态位移列向量，rc(t-T)为铣刀前一个刀齿周期的X方向、Y方向、Z方向的动态位移列向量；
rw(t)为工件在当前刀齿铣削下的X方向、Y方向、Z方向的动态位移列向量，rw(t-T)为工件在前一个刀齿周期铣削下的X方向、Y方向、Z方向的动态位移列向量；
变曲率模具型腔曲线(x(u)、y(u))的曲率表达式为
式中，u为变曲率模具型腔曲线参数方程中的参数变量，ρ(u)为参数u对应点的曲率，变曲率曲线的弯曲方向通过曲线曲率的正负来表示，负数曲率代表的是凹曲线，正数曲率代表的是凸曲线；x(u)、y(u)为笛卡尔坐标下的参数化曲线，x′(u)、y′(u)为在笛卡尔坐标系下的参数化曲线x(u)、y(u)的微分；x″(u)、y″(u)为在笛卡尔坐标系下的参数化曲线x(u)、y(u)的二阶微分；参数u为参数化曲线的自变量；
将刀具铣削的振动位移矢量由时域变换为频域，铣削点的耦合再生位移是
{Δ(jω)}＝{rc(jω)}-{roc(jω)}-({rw(jω)}-{row(jω)})＝(1-ejωT)[Gc(jω)+Gw(jω)]{Fc}ejωT式中，rc(jω)为铣刀当前刀齿的X方向、Y方向、Z方向的频域振动位移矢量，roc(jω)为铣刀当前刀齿的X方向、Y方向、Z方向的时域振动位移矢量，rw(jω)为工件在当前刀齿铣削下的X方向、Y方向、Z方向的频域振动位移矢量，row(jω)为工件在当前刀齿的X方向、Y方向、Z方向的时域振动位移矢量，Fc为铣削力的频域；ωT为产生振动后刀齿周期T之间的相位滞后，得到：{Fc}ejωT＝ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]{Fc}ejωT (13)令ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]行列式值为0，得到ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]的特解：det{[I]-ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]}＝0 (14)式(14)的特征值为：根据(14)和(15)得到ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]的特征方程为：det{[I]+Λ[A(0)][Gc(jω)+Gw(jω)]}＝0 (16)式中，det{·}为方阵函数；
假设只考虑刀具-工件系统的直接传递函数，忽略刀具-工件系统的交叉传递函数，特征方程(16)就简化为一个三次函数：a3Λ3+a2Λ2+a1Λ+1＝0 (17)
求得特征方程的特征值Λ，因为传递函数是复数，a3Λ3+a2Λ2+a1Λ+1＝0的特征值有实部和虚部，因此特征值Λ＝ΛR+iΛI；
式中，ΛR为特征值的实部，ΛI为特征值的虚部，i2＝-1为复数；a1、a2、a3为系数；
式中：K为铣刀铣削圆弧留下的振动波纹的整数，也就是稳定域曲线叶瓣的个数；n为刀具主轴转速。
2.根据权利要求1所述一种预测模具型腔数控铣削中颤振的方法，其特征在于：所述当刀具顺铣时，刀齿的水平切入角和切出角分别为 和 公式为：式中，Rc为刀具半径，ae为径向切削深度；
当刀具逆铣时，刀齿的水平切入角和切出角分别为 和 公式为：
3.根据权利要求2所述一种预测模具型腔数控铣削中颤振的方法，其特征在于：所述刀具子系统的传递函数式中，Gxxc(jω)为刀具子系统在X方向激励、X方向响应的直接传递函数，Gxyc(jω)为刀具子系统在X方向激励、Y方向响应的直接传递函数，Gxzc(jω)为刀具子系统在X方向激励、Z方向响应的直接传递函数，Gyxc(jω)为刀具子系统在Y方向激励、X方向响应的直接传递函数，Gyyc(jω)为刀具子系统在Y方向激励、Y方向响应的直接传递函数，Gyzc(jω)为刀具子系统在Y方向激励、Z方向响应的直接传递函数，Gzxc(jω)为刀具子系统在Z方向激励、X方向响应的直接传递函数，Gzyc(jω)为刀具子系统在Z方向激励、Y方向响应的直接传递函数，Gzzc(jω)为刀具子系统在Z方向激励、Z方向响应的直接传递函数；
式中，Gxxw(jω)为工件子系统在X方向激励、X方向响应的直接传递函数，Gxyw(jω)为工件子系统在X方向激励、Y方向响应的直接传递函数，Gxzw(jω)为工件子系统在X方向激励、Z方向响应的直接传递函数，Gyxw(jω)为工件子系统在Y方向激励、X方向响应的直接传递函数，Gyyw(jω)为工件子系统在Y方向激励、Y方向响应的直接传递函数，Gyzw(jω)为工件子系统在Y方向激励、Z方向响应的直接传递函数，Gzxw(jω)为工件子系统在Z方向激励、X方向响应的直接传递函数，Gzyw(jω)为工件子系统在Z方向激励、Y方向响应的直接传递函数，Gzzw(jω)为工件子系统在Z方向激励、Z方向响应的直接传递函数。
4.根据权利要求3所述一种预测模具型腔数控铣削中颤振的方法，其特征在于：所述步骤三中基于步骤二得到的临界轴向切削深度判断模具型腔数控铣削是否发生颤振；具体过程为：重复执行步骤二得到所有转速下的临界轴向切削深度，将得到的所有临界轴向切削深度aplim定义为数组，利用公式(11)和公式(12)计算刀具刀齿的水平切入角和切出角，再结合刀具主轴转速组成铣削稳定性叶瓣阵列，根据刀具刀齿的水平切入角、切出角以及刀具主轴转速找到铣削稳定性叶瓣阵列中对应的临界稳定的轴向切削深度aplim，aplim与实际轴向切削深度ap相比；
如果ap小于aplim，预测铣削稳定，不发生颤振现象；如果ap大于等于aplim，预测发生颤振现象。
一种预测模具型腔数控铣削中颤振的方法
技术领域
[0001] 本发明涉及预测模具型腔数控铣削中颤振的方法。
背景技术
[0002] 复杂型腔模具广泛应用在汽车、航空航天、船舶和家电等行业，其一般要求有较高的加工精度和表面质量。通常在模具型腔中形状特征多变，具有不规则的尖角、圆角或钝角等大小角度不同的过渡连接和复杂变曲率型腔结构，计算机数控(CNC)铣削加工技术是这类工件的重要加工方式之一。目前计算机数控(CNC)铣削加工技术得到了很好的发展，具有金属切除率高、工件表面质量好和效率高等优点，但对于新产品或部件加工系统仍然需要长时间的试验来获取最佳加工工艺。根据生产阶段的加工条件，仍然存在着对加工结果未知的风险。此外，自适应加工误差控制和补偿系统尚未完全纳入通用数控系统，对于数控加工中心的操作仍依赖于工程师的经验和技术。目前数控铣削加工常用的一些商业软件(例如UG、MasterCAM、Power Mill等)可以进行数控编程，提供加工过程的仿真，但主要限于几何仿真。而模具型腔由于多样化的形状特征，当应用不同的刀具和不同的铣削参数时，其铣削过程刀具与工件接触特性实时变化，铣削力变化大，易发生颤振，导致工件加工质量恶化。颤振一直是数控加工过程中的一个难题，尤其是低刚度铣刀铣削加工高硬度的模具零件，在数控加工企业中很多时候仍然依赖于试切法来确定最佳铣削条件和合适的刀具，这导致企业增加了生产成本，降低了生产效率。
[0003] 尽管有大量学者进行铣削稳定性研究，提出多种铣削稳定性模型，但主要研究集中在单一刀具路径的铣削稳定域生成，实际应用局限于单一刀具轨迹仿真。对于复杂变曲率型腔模具这类铣削条件连续变化的工件而言，单一刀具路径的铣削稳定性预测方法适用性低，导致铣削颤振预测准确度低，加快刀具失效，影响模具型腔的加工质量。
发明内容
[0004] 本发明的目的是为了解决现有单一刀具路径的铣削稳定性预测方法适用性低，导致铣削颤振预测准确度低，加快刀具失效，影响模具型腔的加工质量的问题，而提出一种预测模具型腔数控铣削中颤振的方法。
[0005] 一种预测模具型腔数控铣削中颤振的方法具体过程为：
[0006] 步骤一.建立刀具-工件系统的相对传递函数；
[0007] 分别获取刀具子系统和工件子系统的传递函数
[0008] Gci(jω)＝Xci(jω)/[Fci(jω)] (1)
[0009] Gwi(jω)＝Xwi(jω)/[Fwi(jω)] (2)
[0010] 式中：Fci(jω)和Fwi(jω)分别为刀具子系统刀尖点和工件子系统工件表面铣削点所受的力，由于刀具和工件之间的作用力大小相等、方向相反，即Fci(jω)＝-Fwi(jω)；Xci(jω)和Xwi(jω)分别为在Fci(jω)和Fwi(jω)作用下产生的位移；i＝1，2，3分别表示遵循右手笛卡尔直角坐标系建立的直角坐标系的X方向、Y方向、Z方向；
[0011] 由式(1)和(2)，推导出刀具-工件系统的相对传递函数：
[0012]
[0013] 步骤二.将步骤一获得的刀具-工件系统的相对传递函数引入三维铣削稳定性模型中，得到铣刀铣削颤振频率处的临界轴向切削深度；
[0014] 步骤三、基于步骤二得到的临界轴向切削深度判断模具型腔数控铣削是否发生颤振。
[0015] 本发明的有益效果为：
[0016] 本发明提出了一种用于复杂变曲率模具型腔数控铣削加工中颤振预测的方法，可以达到的效果如下：
[0017] 1.本发明可以进行复杂变曲率模具型腔铣削加工过程颤振预测，针对不同曲率变化的模具型腔，采用不同材质和直径的刀具，选用不同铣削参数(径向切深、主轴转速、轴向切深等)等情况均可进行颤振预测。
[0018] 2.本发明基于相对传递函数综合考虑刀具子系统和工件子系统的动力学特性，同时考虑了横向振动和轴向振动对动态切屑厚度的影响，使颤振预测更加准确。
[0019] 3.本发明在不同曲率形状特征模具型腔的数控几何仿真中，避免了当铣削条件变化时重复生成铣削稳定性叶瓣的必要性，预测过程方便快捷，可以实现铣削颤振快速预测，提高了效率,延长刀具是用寿命，保证了模具型腔的加工质量。
[0020] 解决了现有单一刀具路径的铣削稳定性预测方法适用性低，导致铣削颤振预测准确度低，加快刀具失效，影响模具型腔的加工质量的问题。
[0021] 在普通铣削工程中，判断是否发生颤振，需要借助电脑，通过对公式的编辑运行，得到计算结果，整个过程大约需要耗费30秒，才能得出结果，并且做出相应的调整。而运用铣削稳定性叶瓣阵列的方法在铣削过程中进行预判，比对仅仅需要1秒的时间，缩短了30倍的时间，并且铣削稳定性叶瓣阵列得到之后可以永久使用，省时便利，极大的提高了铣削过程的效率。
附图说明
[0022] 图1为模具型腔铣削不同曲率时刀具-工件接触情况示意图；
[0023] 图2为不同直径刀具铣削时刀具-工件接触情况示意图；
[0024] 图3为不同径向切深下铣削时刀具-工件接触情况示意图，ae1、ae2、ae3为不同的径向切深；
[0025] 图4为模具型腔数控铣削中颤振预测流程示意图；
[0026] 图5为模具型腔铣削加工的刀具路径示意图；
[0027] 图6为实施例获得的传递函数示意图；
[0028] 图7为实施例模具型腔铣削加工的刀具路径示意图；
[0029] 图8为刀具由位置1加工到位置2的切入角和切出角变化示意图；
[0030] 图9为不同主轴转速和切入角时临界轴向切深aplim的变化示意图。
具体实施方式
[0031] 具体实施方式一：本实施方式的一种预测模具型腔数控铣削中颤振的方法具体过程为：
[0032] 步骤一.建立刀具-工件系统的相对传递函数；
[0033] 获取刀具-工件系统的动力学特性，即传递函数矩阵，是进行铣削加工过程稳定性预测的重要前提。以往的研究大多只考虑刀具系统的动力学特性，而忽略了工件系统对整体加工系统动力学特性的影响，本发明综合考虑了刀具子系统和工件子系统的动力学特性，基于相对传递函数建立刀具-工件系统整体动力学模型。
[0034] 分别获取刀具子系统和工件子系统的传递函数
[0035] Gci(jω)＝Xci(jω)/[Fci(jω)] (1)
[0036] Gwi(jω)＝Xwi(jω)/[Fwi(jω)] (2)
[0037] 式中：Fci(jω)和Fwi(jω)分别为刀具子系统刀尖点和工件子系统工件表面铣削点所受的力，由于刀具和工件之间的作用力大小相等、方向相反，即Fci(jω)＝-Fwi(jω)；Xci(jω)和Xwi(jω)分别为在Fci(jω)和Fwi(jω)作用下产生的位移；i＝1，2，3分别表示遵循右手笛卡尔直角坐标系建立的直角坐标系的X方向、Y方向、Z方向；由式(1)和(2)，推导出刀具-工件系统的相对传递函数：
[0038]
[0039] 步骤二.将步骤一获得的刀具-工件系统的相对传递函数引入三维铣削稳定性模型中，
[0040] 得到铣刀铣削颤振频率处的临界轴向切削深度；
[0041] 步骤三、基于步骤二得到的临界轴向切削深度判断模具型腔数控铣削是否发生颤振。
[0042] 具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤二中将步骤一获得的刀具-工件系统的相对传递函数引入三维铣削稳定性模型中，得到铣刀铣削颤振频率处的临界轴向切削深度；具体过程为：
[0043] 基于频域法的三维铣削稳定性模型铣刀铣削模具凹槽型腔侧面的过程中，切向力Ft沿铣刀铣削刃铣削速度方向分布，径向力Fr是铣刀进给的径向方向，轴向力Fa沿铣刀轴向作用；为了考虑一般情况，在铣削刃线中，铣削刃被分成有限数量的小的微分单元。对每个边缘单元的三个方向切屑载荷和相应的微分载荷进行了评估，并用数字积分法分别预测了三个方向的总力。
[0044] 在铣削刃线中，铣削刃被分成M个微分单元，M取值为正整数；对每个边缘微分单元的三个方向切屑载荷和相应的微分载荷进行了评估，并用数字积分法分别预测了三个方向的总力；三个方向即切向力Ft、径向力Fr和轴向力Fa；
[0045] 设第j个铣削刃上第k个微分单元的转角量 表达式为：
[0046]
[0047] 其中：转角量 为第j个铣削刃上第k个微分单元的转角量，θr(k,t)是铣削刃第k个微分单元在笛卡尔坐标Y轴方向测量的旋转角度，0≤k≤M；Nf表示球头铣刀所含有的刀齿数量，n表示刀具的转速，t为时间；j取值为正整数；
[0048] 微分单元区域瞬态铣削层面积 是瞬态切屑厚度，Δa是一个微分单元的长度；
[0049] 切向力Ft、径向力Fr、轴向力Fa在铣刀铣刃上的微分单元是
[0050] dFt＝KtcdAc dFr＝KrcdAc dFa＝KacdAc (5)
[0051] 式中Ktc为切向力的铣削力系数，Krc为径向力的铣削力系数，Kac为轴向力的铣削力系数；
[0052] 切向力Ft、径向力Fr和轴向力Fa转化为刀具笛卡尔坐标下的三个方向的动态铣削力：
[0053]
[0054] 式中 为笛卡尔坐标X轴方向的动态铣削力， 为笛卡尔坐标Y轴方向的动态铣削力， 为笛卡尔坐标Z轴方向的动态铣削力，ap为刀具轴向切削深度， 为转角量，Kr为径向力的铣削力系数与切向力的铣削力系数之比，Kr＝Krc/Ktc；Ka为轴向力的铣削力系数与切向力的铣削力系数之比，Ka＝Kac/Ktc；Δxj＝(xc(t)-xc(t-T))-(xw(t)-xw(t-T))，Δyj＝(yc(t)-yc(t-T))-(yw(t)-yw(t-T))，Δzj＝(zc(t)-zc(t-T))-(zw(t)-zw(t-T))，Δxj、Δyj、Δzj均为再生效应振动；
[0055] 其中xc(t)为铣刀当前刀齿的X方向(笛卡尔坐标)动态位移，xc(t-T)为铣刀前一个刀齿周期的X方向(笛卡尔坐标)动态位移，xw(t)为工件在当前刀齿铣削下的X方向(笛卡尔坐标)动态位移，xw(t-T)为工件在前一个刀齿周期铣削下的X方向(笛卡尔坐标)动态位移；
[0056] yc(t)为铣刀当前刀齿的Y方向(笛卡尔坐标)动态位移，yc(t-T)为铣刀前一个刀齿周期的Y方向(笛卡尔坐标)动态位移，yw(t)为工件在当前刀齿铣削下的Y方向(笛卡尔坐标)动态位移，yw(t-T)为工件在前一个刀齿周期铣削下的Y方向(笛卡尔坐标)动态位移；
[0057] zc(t)为铣刀当前刀齿的Z方向(笛卡尔坐标)动态位移，zc(t-T)为铣刀前一个刀齿周期的Z方向(笛卡尔坐标)动态位移，zw(t)为工件在当前刀齿铣削下的Z方向(笛卡尔坐标)动态位移，zw(t-T)为工件在前一个刀齿周期铣削下的Z方向(笛卡尔坐标)动态位移；
[0058] T为周期；
[0059] [F(t)]＝apKtc[A(t)]{Δr} (7)
[0060] 其中：F(t)为以 为元素的列向量，{Δr}为以Δxj、Δyj、Δzj为元素的列向量，动态铣削力方向系数矩阵[A(t)]是角频率为ω，周期为T的周期函数；ω＝NfΩ，T＝2π/ω；Ω为刀具主轴转动角速度，单位为rad/s；
[0061] 将[A(t)]进行傅里叶级数展开，高阶的傅里叶变换会得到更精确的解，但在实践中迭代特征值搜索并不可取，得到平均方向系数矩阵：
[0062]
[0063] 式中，Nf表示球头铣刀所含有的刀齿数量， 为刀具刀齿的切入角， 为刀具刀齿的切出角， 为系数矩阵， 为转角量；
[0064] 动态铣削过程方程的方向系数是根据刀具和工件的接触条件改变，由式(8)，将时变动态铣削方程[F(t)](公式7)转换为不随时间变化但与刀具-工件接触状态相关的形式：
[0065] [F(t)]＝apKtc[A(0)]{Δr} (9)
[0066] 其中，
[0067] {Δr}＝{rc(t)}-{rc(t-T)}-({rw(t)}-{rw(t-T)})
[0068] 式中，rc(t)为铣刀当前刀齿的X方向、Y方向、Z方向(笛卡尔坐标)的动态位移列向量，rc(t-T)为铣刀前一个刀齿周期的X方向、Y方向、Z方向(笛卡尔坐标)的动态位移列向量；
[0069] rw(t)为工件在当前刀齿铣削下的X方向、Y方向、Z方向的动态位移列向量，rw(t-T)为工件在前一个刀齿周期铣削下的X方向、Y方向、Z方向的动态位移列向量；
[0070] 在单一刀具路径铣削的情况下，在铣削参数确定的情况下，刀具切入角和切出角不变，刀具-工件的接触状态固定。而在复杂变曲率模具型腔铣削加工过程中，随着型腔曲线曲率发生变化时，刀具-工件的接触状态随之变化，刀具切入角和切出角不同，如图1所示。此外，在采用不同刀具半径和不同径向切深时，刀具-工件的接触状态同样变化，刀具切入角和切出角不同，如图2、图3所示。因此，应确定模具型腔变曲率曲线铣削刀具-工件接触状态和曲率、刀具半径、径向切深之间的关系；
[0071] 变曲率模具型腔曲线(x(u)、y(u))的曲率表达式为
[0072]
[0073] 式中，u为变曲率模具型腔曲线参数方程中的参数变量，ρ(u)为参数u对应点的曲率，变曲率曲线的弯曲方向通过曲线曲率的正负来表示，负数曲率代表的是凹曲线，正数曲率代表的是凸曲线；x(u)、y(u)为笛卡尔坐标下的参数化曲线，x′(u)、y′(u)为在笛卡尔坐标系下的参数化曲线x(u)、y(u)的微分；x″(u)、y″(u)为在笛卡尔坐标系下的参数化曲线x(u)、y(u)的二阶微分；参数u为参数化曲线的自变量；
[0074] 当刀具顺铣时，刀齿的水平切入角和切出角分别为 和
[0075] 当刀具逆铣时，刀齿的水平切入角和切出角分别为 和
[0076] 设刀具子系统的频响函数为Gc(jω)；
[0077] 工件子系统的频响函数为Gw(jω)；
[0078] 将刀具铣削的振动位移矢量由时域变换为频域，铣削点的耦合再生位移是[0079] {Δ(jω)}＝{rc(jω)}-{roc(jω)}-({rw(jω)}-{row(jω)})＝(1-ejωT)[Gc(jω)+Gw(jω)]{Fc}ejωT
[0080] 式中，rc(jω)为铣刀当前刀齿的X方向、Y方向、Z方向(笛卡尔坐标)的频域振动位移矢量，roc(jω)为铣刀当前刀齿的X方向、Y方向、Z方向(笛卡尔坐标)的时域振动位移矢量，rw(jω)为工件在当前刀齿铣削下的X方向、Y方向、Z方向(笛卡尔坐标)的频域振动位移矢量，row(jω)为工件在当前刀齿的X方向、Y方向、Z方向(笛卡尔坐标)的时域振动位移矢量，Fc为铣削力的频域；
[0081] ωT为产生振动后刀齿周期T之间的相位滞后，得到：
[0082] {Fc}ejωT＝ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]{Fc}ejωT (13)[0083] 令ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]行列式值为0，得到ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]的特解：
[0084] det{[I]-ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]}＝0 (14)
[0085] 式(14)的特征值为：
[0086]
[0087] 根据(14)和(15)得到ap[A(0)](1-ejωT)[Gc(jω)+Gw(jω)]的特征方程为：
[0088] det{[I]+Λ[A(0)][Gc(jω)+Gw(jω)]}＝0 (16)
[0089] 式中，det{·}为方阵函数；
[0090] 根据给定的铣刀铣削颤振频率(已知的)，切入角切出角(公式11或12)和刀具-工件系统的传递函数(公式3)，求得式(16)的特征值，假设只考虑刀具-工件系统的直接传递函数，忽略刀具-工件系统的交叉传递函数，特征方程(16)就简化为一个三次函数：
[0091] a3Λ3+a2Λ2+a1Λ+1＝0 (17)
[0092] 求得特征方程的特征值Λ，因为传递函数是复数，a3Λ3+a2Λ2+a1Λ+1＝0的特征值有实部和虚部，因此特征值Λ＝ΛR+iΛI；
[0093] 式中，ΛR为特征值的实部，ΛI为特征值的虚部，i2＝-1为复数；a1、a2、a3为系数；
[0094] 求得铣刀铣削颤振频率处的临界轴向切削深度为：
[0095]
[0096] 刀具主轴转速由延迟周期换算得到：
[0097]
[0098] 式中：K为铣刀铣削圆弧留下的振动波纹的整数，也就是稳定域曲线叶瓣的个数；n为刀具主轴转速。
[0099] 其它步骤及参数与具体实施方式一相同。
[0100] 具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述当刀具顺铣时，刀齿的水平切入角和切出角分别为 和 公式为：
[0101]
[0102] 式中，Rc为刀具半径，ae为径向切削深度；
[0103] 当刀具逆铣时，刀齿的水平切入角和切出角分别为 和 公式为：
[0104]
[0105] 其它步骤及参数与具体实施方式一或二相同。
[0106] 具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述刀具子系统的频响函数
[0107] 式中，Gxxc(jω)为刀具子系统在X方向(笛卡尔坐标)激励、X方向(笛卡尔坐标)响应的直接频响函数，Gxyc(jω)为刀具子系统在X方向(笛卡尔坐标)激励、Y方向(笛卡尔坐标)响应的直接频响函数，Gxzc(jω)为刀具子系统在X方向(笛卡尔坐标)激励、Z方向(笛卡尔坐标)响应的直接频响函数，Gyxc(jω)为刀具子系统在Y方向(笛卡尔坐标)激励、X方向(笛卡尔坐标)响应的直接频响函数，Gyyc(jω)为刀具子系统在Y方向(笛卡尔坐标)激励、Y方向(笛卡尔坐标)响应的直接频响函数，Gyzc(jω)为刀具子系统在Y方向(笛卡尔坐标)激励、Z方向(笛卡尔坐标)响应的直接频响函数，Gzxc(jω)为刀具子系统在Z方向(笛卡尔坐标)激励、X方向(笛卡尔坐标)响应的直接频响函数，Gzyc(jω)为刀具子系统在Z方向(笛卡尔坐标)激励、Y方向(笛卡尔坐标)响应的直接频响函数，Gzzc(jω)为刀具子系统在Z方向(笛卡尔坐标)激励、Z方向(笛卡尔坐标)响应的直接频响函数；
[0108] 工件子系统的频响函数
[0109] 式中，Gxxw(jω)为工件子系统在X方向(笛卡尔坐标)激励、X方向(笛卡尔坐标)响应的直接频响函数，Gxyw(jω)为工件子系统在X方向(笛卡尔坐标)激励、Y方向(笛卡尔坐标)响应的直接频响函数，Gxzw(jω)为工件子系统在X方向(笛卡尔坐标)激励、Z方向(笛卡尔坐标)响应的直接频响函数，Gyxw(jω)为工件子系统在Y方向(笛卡尔坐标)激励、X方向(笛卡尔坐标)响应的直接频响函数，Gyyw(jω)为工件子系统在Y方向(笛卡尔坐标)激励、Y方向(笛卡尔坐标)响应的直接频响函数，Gyzw(jω)为工件子系统在Y方向(笛卡尔坐标)激励、Z方向(笛卡尔坐标)响应的直接频响函数，Gzxw(jω)为工件子系统在Z方向(笛卡尔坐标)激励、X方向(笛卡尔坐标)响应的直接频响函数，Gzyw(jω)为工件子系统在Z方向(笛卡尔坐标)激励、Y方向(笛卡尔坐标)响应的直接频响函数，Gzzw(jω)为工件子系统在Z方向(笛卡尔坐标)激励、Z方向(笛卡尔坐标)响应的直接频响函数。
[0110] 其它步骤及参数与具体实施方式一至三之一相同。
[0111] 具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤三中基于步骤二得到的临界轴向切削深度判断模具型腔数控铣削是否发生颤振；具体过程为：
[0112] 从模具型腔几何模型铣削数控代码和频域稳定性叶瓣构造算法出发，推导了颤振稳定性预测算法。该预测模具型腔铣削颤振的方法具体实施过程如图4所示。
[0113] 获取刀具子系统和工件子系统的动力学特性，得到其传递函数。刀具子系统传递函数通过模态实验、RCSA或理论建模等方法获取，工件子系统传递函数通过模态实验、有限元仿真等方法获取；再将两个子系统的传递函数导入公式(3)中，得到刀具-工件系统相对传递函数；步骤一已经得到；
[0114] 根据常用刀具的信息、常见工件信息、径向切削深度范围等确定刀具切入角和切出角的变化范围，该范围尽量覆盖模具型腔铣削的常见情况，再结合主轴转速范围和上一步所获取的刀具-工件系统相对传递函数生成稳定性叶瓣阵列；
[0115] 对于某一给定的模具型腔铣削加工条件，根据主轴转速和铣刀切入切出角度确定铣削稳定性，将步骤二得到的临界稳定的轴向切削深度aplim定义为数组，即在主轴转速范围和切入角切出角范围内，任一的主轴转速和切入角切出角都对应有一个临界稳定的轴向切削深度aplim。步骤二已经得到；
[0116] 针对给定的模具工件型腔的数控铣削加工，按照通常加工经验给定的一个初始数控铣削条件，即选用刀具、选定铣削参数、选定加工策略和规划刀具路径。图5所示为进行模具型腔铣削加工的刀具路径，在此条件下利用公式(11)和公式(12)计算刀具与工件之间的接触情况，获得刀具刀齿切入角和切出角，再结合已知刀具主轴转速信息在铣削稳定性叶瓣阵列中找到对应的临界稳定的轴向切削深度aplim；
[0117] 颤振预测算法是指在NC几何仿真模型中的主轴转速和切入切出角度下的铣削稳定性阵列中的稳定临界轴向切削深度(aplim)与实际轴向切削深度(ap)进行严格比较；
[0118] 重复执行步骤二得到所有转速下的临界轴向切削深度，将得到的所有临界轴向切削深度aplim定义为数组，利用公式(11)和公式(12)计算刀具刀齿切入角和切出角，再结合刀具主轴转速组成铣削稳定性叶瓣阵列，根据刀具刀齿切入角、切出角以及刀具主轴转速找到铣削稳定性叶瓣阵列中对应的临界稳定的轴向切削深度aplim，aplim与实际轴向切削深度ap相比；
[0119] 如果ap小于aplim，预测铣削稳定，不发生颤振现象；如果ap大于等于aplim，预测发生颤振现象。
[0120] 采用以下实施例验证本发明的有益效果：
[0121] 实施例一：
[0122] 本实施例一种预测模具型腔数控铣削中颤振的方法具体是按照以下步骤制备的：
[0123] 将该方法应用于模具型腔铣削加工，按照步骤一获得刀具子系统传递函数、工件子系统传递函数，在此基础上得到刀具-工件系统的相对传递函数，如图6所示。
[0124] 使用CAM软件按规划的刀具路径生成数控代码，图7表示生成的刀具路径，初选的轴向铣削深度是0.7mm，径向铣削深度是1mm，刀具直径是20mm，计算刀具与工件接触特性，图8表示由位置1加工到位置2的时候刀具切入角和切出角的变化情况。
[0125] 随着切入角的减小，稳定性曲线的临界轴向切削深度，如图9所示。对于给定的模具型腔工件和刀具，可以在数控仿真之前生成颤振稳定性阵列 本发明方法根据刀具-工件接触情况计算 的变化，并在颤振稳定性阵列 中提取出
对应的aplim，与实际的ap进行比较。
[0126] 当型腔铣削主轴转速是5500rpm，不同切入角和切出角时临界轴向切深aplim如图9所示，在给定的加工位置时可以得到切入角(1.9rad)和切出角(3.14rad)的临界轴向切深aplim是0.56mm，而初选轴向切深ap是0.7mm，高于aplim，所以该方法预测在主轴转速5500rpm和切入角切出角(1.9rad,3.14rad)下加工变得不稳定，发生颤振。
[0127] 当型腔铣削主轴转速是7500rpm，在图9中可以得到切入角(1.9rad)和切出角(3.14rad)的临界轴向切深aplim是0.89mm，高于初选的轴向切深ap。该方法预测在主轴转速7500rpm和切入角切出角(1.9rad,3.14rad)下加工保持稳定。
