考虑模态耦合的球头铣刀颤振稳定域建模方法*

2019-06-24代月帮李宏坤魏兆成

振动、测试与诊断 2019年3期

代月帮，李宏坤，魏兆成

(大连理工大学机械工程学院大连，116024)

1 问题的引出

铣削过程中不合适的切削参数导致的颤振严重地影响到了加工效率、精度、质量以及稳定性，制约了铣削技术的快速发展。在多种颤振形成的机理模型中，再生型颤振被认为是切削过程中产生颤振最直接、最根本的原因。如图1所示，由于机床结构振动，上一次切削与本次切削振动位移之间的相位差异导致刀具切削厚度不均而引起自激振动，即为再生型颤振[1]。

图1 再生型颤振的理论模型Fig.1 Mechanism model of regenerative chatter

对铣削动力学方程进行时域离散化处理，获得振动位移的时序表达式，绘制出轴向临界切削深度随主轴转速变化的叶瓣图是目前避免颤振最有效的方法。国内外学者研究了多种颤振时域数值求解算法，如：基于龙格库塔的时域求解法[2]；全离散时域求解法[3]；完全离散时域数值求解法[4]。但是，不论是龙格库塔法、全离散法还是改进的完全离散法，本质上都属于差分类方法，其优点是易于实现，但差分引起的误差很难从根本上消除[5]，同时它们所建立的模型忽略了刀具系统的模态耦合因素。

针对上述不足，笔者利用精细积分法构造出考虑模态耦合的球头铣刀颤振稳定域叶瓣图，并对其进行试验验证，克服传统模型中存在的弊端。

2 球头铣刀-工件动力学方程

2.1 动力学方程的建立

考虑进给方向x和法向y方向的模态耦合，建立如式(1)所示的方程

(1)

其中：

2.2 动态切削力求解

2.2.1 切削刃几何模型

式(2)为第j刀刃上第i个微元坐标表达式[6-7]

(2)

其中：R为球头铣刀半径；β为切削刃螺旋角；t为切削过程中的时间；k为切削微元的轴向接触角；ψji(k)为切削微元的径向滞后角；φ10(t)为刀具转动角度；n为刀具转速；φji(t)为切削微元的径向接触角；Nf为切削刃数目。

2.2.2 动态切削力模型

切向力dFt,ji(φji(t),k)、径向力dFr,ji(φji(t),k)和轴向力dFa,ji(φji(t),k)表示为

(3)

其中：h(φji(t),k)为瞬时动态切削厚度；Ktc为切向力系数；Krc为径向力系数；Kac为轴向力系数；db为瞬时切削宽度，db=Rdk；R为铣刀球头半径。

2.2.3 动态切削厚度

切削微元的动态切削厚度如式(4)所示

(4)

其中：x(t)-x(t-T)，y(t)-y(t-T)为t和(t-T)时刻在x和y方向的振动矢量；T为时滞周期，T=60/(Nfn)。

2.2.4 动态切削力

通过坐标变换，获得切削微元x，y方向瞬时动态切削力为

(5)

其中：M为坐标转换矩阵。

确定切削刀刃最小和最大轴向角，即可得到如式(6)所示的动态切削力

(6)

其中

3 切触区域建模

3.1 切触区域边界构成

图2所示的球头铣刀-工件切触区域边界主要由如图3所示的a，b，c，d这4条线组成[7]。

如图4所示，以铣刀球头顶点为原点，构建三维笛卡尔坐标系x-y-z，其中：z轴为铣刀轴线，远离球头铣刀的顶点为正方向；x轴为进给方向，正方向指向待加工表面。在该坐标系下，铣刀球头轮廓的方程为x2+y2+(z-R)2=R2，其中R为铣刀球头半径。

图2 切触区域Fig.2 The schematic for contact area of ball-end mill-workpiece

图3 切触区域边界组成Fig.3 The schematic for borderline of contact area

图4 三维坐标系建立示意图Fig.4 The schematic for 3D coordinate system

3.2 求解切触区域边界的投影方程

1)a号线投影方程

在x-y平面内，a号线投影方程为x2+y2=R2-(R-L_jg)2，L_jg为切削深度。

2)b号线投影方程

在x-y平面内，b号线投影为1条与y轴重合的直线，方程为x=0。

3)c号线投影方程

如图5所示，在y-z坐标下，获取球头截面方程y2+(z-R)2=R2，将其沿y轴负方向平移刀轨间距L_xl，即可得到与本次刀位点相对应的上一刀轨的铣刀球头截面方程，该方程为(y+L_xl)2+(z-R)2=R2。y2+(z-R)2=R2与(y+L_xl)2+(z-R)2=R2的交点M点为c号线在y-z平面的投影，因此c号线投影方程为y=-L_xl/2。

图5 M点示意图Fig.5 The schematic for point M

4)d号线投影方程

d号线上的其它点的投影坐标可以认为是不同轴向切削深度所对应的a号线与e号线的投影方程的交点。因此，通过插值理论，可以获得d号线多项式方程。

如图6所示，a，b，c，d投影共同围成的区域即为工件与球头铣刀切触区域在x-y平面的投影。

图6 切触区域在x -y平面的投影Fig.6 The projection of the contact region in the x-y plane

4 考虑模态耦合的动力学方程时域数值求解

4.1 基于精细积分法的动力学方程全离散

由式(6)可将式(1)表示为

(7)

将式(7)表示为如下所示的哈密顿系统

(8)

其中：

将A(t)v(t)-A(t)v(t-T)用f(t)来表示，则对于非齐次方程(8)，由常微分理论可知，一般解[8]为

(9)

将时滞周期T均分为m份，即T=mτ，在[tp,tp+1]中，将f(t)表示为如下形式

f(t)=r0+r1(t-tp)

(10)

其中：

r0=f(tp)=A(tp)v(tp)-A(tp)v(tp-mτ)；

将r0，r1进一步表示为

r0=Apvp-Apvp-m

(11)

(12)

由式(9)和式(10)可将v(tp+1)表示为

v(tp+1)=T1[v(tp)+A0-1(r0+A0-1r1)]-

A0-1(r0+A0-1r1+r1τ)

(13)

其中：T1=eA0ε。

将式(11)、式(12)带入到式(13)中，可得

(I-NN/τAp+1)vp+1=(T1+MMAp-NN/τAp)·

vp-NN/τAp+1vp+1-m+(NN/τAp-

MMAp)vp-m

(14)

其中：MM=T1A0-1-A0-1；NN=T1A0-2-A0-2-A0-1τ。

若(I-NN/τAp+1)可逆，则式(14)可表示为

vp+1=(I-NN/τAp+1)-1(T1+MMAp-

NN/τAp)vp-(I-NN/τAp+1)-1NN/

τAp+1vp+1-m+(I-NN/τAp+1)-1(NN/

τAp-MMAp)vp-m

(15)

4.2 瞬时切削位置确定方法

将铣刀全部切削刃在tp时刻投影到x-y坐标系下。由切削刃投影方程与接触区域边界方程的关系便能够得到在此刻刀刃所对应的最大和最小的轴向角。通过该方法能够获得整个周期内所有离散时刻的每一个刀刃所对应的最大和最小的轴向角[7]。

5 叶瓣图构建方法

构建出矩阵Cp，满足如下映射

vp+1=Cpvp

其中：

矩阵Cp中的PK为4×4矩阵，等于式(15)中的(Ι-NN/τAp+1)-1(T+MMAp-NN/τAp)；RK1为4×2矩阵，为式(15)中- (I-NN/τAp+1)-1NN/τAp+1的前2列；RK2为4×2矩阵，等于式 (15)中(I-NN/τAp+1)-1(NN/τAp-MMAp)的前2列。

通过使用一系列离散Cp(p=0,1,…,m-1)，获得整个周期内的过渡矩阵Φ，亦即

vp=Φv0

其中：Φ=Cm-1Cm-2…C1C0。

由Floquet理论可知，矩阵Φ特征值模的最大值小于1、等于1和大于1时，分别表示切削处于稳定状态、临界状态和不稳定状态。

改变主轴转速，获得不同主轴转速所对应的临界切削深度，最终构建出稳定域叶瓣图。