窦 炜，袁胜万，何晓聪

(1.昆明理工大学 机电工程学院，昆明 650500；2.沈机集团昆明机床股份有限公司，昆明 650000)

铣削加工过程中，铣刀各齿周期性的切削工件产生大小和方向不断变化的动态铣削力，引起刀具与工件的受迫振动。振动会使刀具在工件表面留下振纹，影响表面粗糙度，并且在一定不利条件下有可能使铣削力产生不稳定的振荡，进而导致整个铣削系统失稳。为了获得高质量的加工表面，需要将刀具相对工件的振动状态控制在可以接受的范围内。利用动力学模型研究铣削过程的动态特性和稳定性是铣削研究的重要手段。

当铣床主轴带动刀具旋转时，同一刀齿上位于刀具轴向不同位置的刀刃将先后切入或切出工件，其时间间隔或延迟取决于刀刃之间的轴向距离以及刀齿螺旋角。因此，在某一时刻，刀齿有可能只有部分刀刃参与切削，这将导致该刀齿上不同轴向位置的刀刃所受切削力各异，并且不同时刻作用于该刀齿的总切削力也可能不是常数。

针对立铣加工的上述特点，Altintas和Engin等[1-3]较早对不同类型的立铣刀(Helical end mill tool)铣削过程开展了系统性研究，给出了一种不同于经典的直齿铣刀模型的螺旋齿铣刀模型，该模型可以计算不同几何特征的立铣刀铣削力，并进一步归纳总结了立铣加工过程中的 “刀具-工件” 动力学模型及分析的基本方法。Balachandran等[4-5]对不同径向切深条件下的铣削加工过程的动力学问题进行了研究，分析了小径向切深时铣刀与工件的断续接触对系统稳定性的影响。Merdol[6]在采用Budak等的方法分析小径向切深时的铣削加工稳定性得到了相同的结果。但是他们的研究局限于小轴向切入深度，因此刀刃切入和切出“延迟”效应不明显。Insperger等[7-8]对比分析了大轴向切入深度条件下的经典铣刀模型(直齿)和螺旋齿铣刀模型的稳定性表现，指出后者的稳定性叶瓣图中存在封闭的不稳定“孤岛”，并通过试验证明了这一现象。

因此，小轴向切深立铣时若忽略刀刃切入和切出“延迟”效应，即认为铣削过程中参与切削的刀刃段恒等于刀具的轴向切入深度，螺旋齿铣削模型近似退化为直齿铣削模型，铣削分析的结果在一定程度上仍然能够与实际加工情况相吻合。但是，针对诸如加工薄壁构件的大轴向切入深度的铣削研究，则必须考虑螺旋齿铣削时的切削力时变特点。并且，若忽略了刀刃切削部分是随着刀具旋转逐渐增加并逐渐减少的实际过程，铣削力中还将包含失真的高频成分的影响，从而导致后续分析产生错误结果。

随着学者们对铣削的研究不断深入，提出的铣削力模型以及研究的问题越来越复杂[9-15]。然而，在铣削加工动力学分析或求解过程中，螺旋齿铣刀铣削力的计算仍普遍采用将刀具沿轴向离散成微元段逐个计算然后求和的近似方法，这是一个相对耗时的环节。本文推导得到螺旋齿铣刀切削力的一种解析表达式，简化了铣削力的计算过程，从而提高铣削动力学方程的求解效率。

1 立铣加工动力学模型

虽然不同几何特征的立铣刀局部切削力有所区别，但是铣削过程中刀具和工件之间的动力学关系是相同的。如图1所示，铣削加工过程可以抽象为包含“动态切削力”和“刀具-工件系统”两个相互作用的环节所形成的反馈系统。

图1 铣削加工过程方框图Fig.1 Block diagram of the milling process

下面以螺旋齿圆柱铣刀为例简要概括立铣加工的动力学模型。图2所示为二自由度立铣加工力学模型。在此假设工件为刚体，刀具相对工件为柔性体，且由于刀具轴向刚度通常远远高于横切面刚度，因此只考虑铣削力作用于刀具X和Y方向的情况。设刀具有N个刀齿，记第j个刀齿所受切削力铣削力Fj，则根据牛顿定律，铣刀中心的位移向量q满足下式：

图2 二自由度“刀具-工件”系统力学模型Fig.2 Mechanical model of the 2DOF tool-workpiece system

(1)

其中

mx和my，cx和cy，kx和ky分别表示刀具在X方向和Y方向的模态质量，模态阻尼和模态刚度。切削力Fj包含由加工进给引起的“静态力”(Static milling force)和再生效应(Regenerative effect)引起的动态力(Dynamic milling force)。x和y以及δx和δy则分别表示铣刀在铣削力作用下的“静态位移”以及再生振动。

图3 螺旋齿立铣刀及微元段Fig.3 Geometry of a helical tool and an infinitesimal element

以第j个刀齿为研究对象。如图3所示，刀齿最前端Pj,1点Z轴坐标为0，取刀齿上距Pj,1点的轴向距离为z处的一个厚度等于dz微元段刀刃进行分析。刀具转角以Y轴正向起始，顺时针方向为正。记该处刀刃的瞬时角位移为φj(t,z)，其值可由如下函数关系求出：

(2)

式中：变量t表示时间，s。参数ω表示铣床主轴转速，rad/s。β为铣刀螺旋角度，rad；r为刀具半径,m。则该微元段刀齿受到的沿刀刃的切向和法向的铣削分力，ft和fn，可以用瞬时切屑面积和比切削力(Specific cutting force)的乘积表示为:

(3)

(4)

式(3)中，Kt和Kn分别表示刀刃切向和法向线性比切削力系数。h表示瞬时切屑厚度，其值约等于铣削过程中静态和动态切屑厚度沿刀齿法向的投影，即：

(5)

式中:fz表示铣刀进给量。参数τ表示刀齿经过周期 (Tool passing period)。窗函数gj(t,z)的物理意义是：如果φj(t,z)介于刀齿切入角φs和切出角φe形成的区间内，则表示该处刀刃进入切削区参与切削，否则不产生切削力。

随后将ft和fn投影至X-Y坐标系下，整理得到微元段刀刃受到的作用力dFj的表达式为：

(6)

其中矩阵Gj为：

于是，第j齿所受的总铣削力Fj等于所有微元段所受切削力的总和：

(7)

式(7)的积分上限为轴向切入深度b。将式(7)代入式(1)得：

(8)

式(8)即为描述螺旋齿铣刀立铣加工过程的时滞微分控制方程(Delay differential equations)。这类方程的一个主要特点是：方程所描述的动态系统的状态向量q的变化率，同时取决于当前和过去的状态向量值。

式(7)中的窗函数gj(t,z)给计算Fj带来了困难。不论采用符号运算(Symbolic computation)或是数值积分(Quadrature)直接计算式(7)都非常繁琐。因此，在实际运用中，很少通过积分的方法计算式(7)，而是采用离散求和近似的策略。具体而言，取足够大的正整数n，将铣刀沿轴向切入深度b离散为n段，则第k段离散微元的角位移φj,k和所受切削力分别为：

(9)

(10)

于是，当n趋于无穷大时Fj等于：

(11)

取式(11)中的有限项求和(如n=100)替换式(8)中的积分，得到：

(12)

这种处理方法的优点在于：在每一个微元段内，窗函数gj化为常数，而计算式(10)只需初等函数的代数运算。随后，式(11)的微元段连加运算可以非常容易地通过循环算法编程实现。缺点是如果切削力计算精度要求较高，则n值必须取的很大，而n值越大花费的计算时间就越多，因此这一步是限制铣削动力学方程求解效率的重要因素之一。

2 改进的铣削力计算方法

如果式(7)中的定积分被积函数不包含窗函数gj，铣削力计算将非常容易。式(2)表明，任意时刻刀刃微元的轴向位置z和该位置转角存在一一对应的函数关系。因此可以通过变量代换，利用铣削过程中实际进入切削区的刀刃段的角位移作为积分域，从而将被积函数中的窗函数消去。

(13)

图4 螺旋齿立铣刀逆铣加工Fig.4 Up milling with a helical end tool

式中：max和min分别代表最大值和最小值函数。φj表示第j齿前端Pj,1的瞬时角位移，其值由下式计算：

(14)

式中：mod为求余函数。式(13)中的参数φd称为切削滞后角度，它表示第j齿前端Pj,1开始切削工件到刀齿切入末端Pj,2进入切削区域的这段时间内刀具所转过的角位移，如图4中的角∠P2,1OP2,2，其值由下式计算：

(15)

(16)

相应地，进入切削区的刀刃微元的厚度dz等于：

(17)

将式(16)和(17)中定义的φj,k和dz代入式(10)和(11)，对求和项取极限，得到Fj计算式如下：

(18)

经过式(18)的处理，式(7)被积函数中判断刀齿是否进入切削区的窗函数gj被新的积分上下限取代了。式(18)的积分是变限积分，接下来需要确定积分限。

由φd的表达式(15)可知，若刀具半径和螺旋角不变，切削滞后角度将是轴向切入深度b的函数，而铣削过程中的刀具切入角和切出角取决于径向切入深度。因此，根据轴向和径向切入深度的不同取值组合，该积分有以下两种情况。当φs+φd<φe时:

(19)

当φs+φd≥φe时:

(20)

对式(19)和(20)进行积分并整理后，上述两种情况的积分结果可统一表示为：

(21)

矩阵Aj(t)各元素的表达式分别为：

aj,xx(t)=-sin(φL-φU)sin(φU-φL)-

aj,xy(t)=-(φU-φL)-cos(φL+φU)×

aj,yx(t)=φU-φL-cos(φL+φU)×

aj,yy(t)=φU-φL-cos(φL+φU)×

(22)

(23)

将式(23)代入式(1)得：

(24)

3 求解立铣加工动力学方程

目前，只有少量表达式非常简单的迟滞微分方程能给出解析解，形如式(12)或(24)的方程只能通过数值方法求解。铣削动力学方程的求解过程可以概括分解为如下步骤：

1) 选择方程求解时间步长Δt，将时间变量t取值范围离散为时刻ti序列(i=1,2,…)。如将刀齿经过周期τ离散为s段，则Δt=τ/s。

2) 计算ti时刻的切屑厚度h(ti)或H(ti,ti-s):

(25)

3) 计算铣刀所受的铣削力Fj(ti)。

4) 计算铣削力Fj(ti)作用下的铣刀位移q(ti)；

5) 将时刻ti值增加一个步长，并重复第2和第4步，直到最后一个离散时刻。

接下来通过仿真，对比利用解析表达式(23)和离散近似表示式(11)求解铣削动力学方程的效果。两种求解方法均采用Matlab 9.0编程，除了切削力计算部分的代码之外，两种方法其余部分程序内容完全相同并且在同一台电脑上运行(AMD A10-6800K4.4GHz，8GB RAM)。

仿真铣削模型参数来自文献[8]：刀具X和Y方向模态参数相同，频率319.38 Hz，阻尼比1.96%，刚度2.16×107N/m；刀具直径20 mm，齿数为4，螺旋角30°；逆铣加工，径向切深1 mm，进给量0.1 mm/齿，比切削力系数Kt= 804.3×106N/m2，Kn= 331×106N/m2。

图5所示为轴向切深b=20 mm，主轴转速Ω=10 000 r/min时刀具所受Y向切削力仿真结果，图(b)为图(a)的局部放大图形，其中虚线和细实线分别表示采用式(23)和式(11)的结果。图(a)中切削力最大幅值保持在95 N上下，表明模型处于稳定铣削状态。图(b)中的细实线"毛刺"显示了采用式(11)求解方程的离散误差。

图5 Y向铣削力仿真结果(b=20 mm，Ω=10 000 r/min)Fig.5 Simulation results for Y direction force

如式(11)所示，随着轴向微元段数目n趋于无穷大，该误差将收敛为0。取时间步长Δt= τ/200，自变量范围0～0.2 s，采用式(23)求解方程耗时0.672 4 s。以此结果作为基准，图6显示了采用式(11)求解方程时刀具离散微元段数目n对仿真结果的影响。

图6中左边y坐标轴表示相对误差，右边y坐标轴表示两种方法求解消耗时间的比值。从图中可以看出，采用离散近似的方法计算铣削力，耗时与n成正比。n=50的时候，相对误差为6%，相对时间为2.1倍。当计算误差小于0.5%，则n需大于500，耗时约为本文提出方法的18倍以上。

图6 刀具离散微元段数对仿真的影响Fig.6 Effect of the number of elementary tool disks on simulation

图7所示为轴向切深b=12 mm，主轴转速Ω= 10 000 r/min时采用式(23)的仿真结果。图(a)和图(b)所示分别为Y向铣削力和位移。可以看出模型此时处于不稳定切削状态。

图7 不稳定铣削条件的仿真结果(b=12 mm，Ω=10 000 r/minFig.7 Simulation result of unstable milling

和图5情况相比，主轴转速不变，减小轴向切深，系统反而从稳定状态变为不稳定。这是因为图7的铣削条件处于螺旋齿铣刀的不稳定“岛”内[8]，如图8中的阴影区域所示。这也验证了式(22)～(24)的有效性。

图8 铣削稳定性叶瓣图Fig.8 Stability lobe diagram

4 结 论

本文基于线性切削力模型，给出了一种螺旋齿铣刀铣削力的一种解析计算表达式。利用该表达式可以有效简化螺旋齿立铣加工动力学方程的求解或分析过程：

(1) 采用本文所提方法求解铣削动力学方程，不需要对刀具沿轴向离散为微元，因此不存在轴向离散带来的计算误差，也因此不需要在求解方程的每一个时间步进行轴向循环计算切削力，从而提高计算精度和效率。采用本文方法求解螺旋齿立铣动力学方程所需计算量与求解直齿铣削动力学方程相当。

(2) 大部分基于时域解的铣削过程稳定性分析需要针对不同主轴转速和轴向切深求解铣削动力学方程，因此计算非常耗时。半离散法或全离散法虽然不直接采用数值方法求解铣削动力学方程，但是也需要对一个时滞周期内的铣削力进行时域离散近似。因此，本文所提方法可以提高螺旋齿立铣稳定性分析效率。

参 考 文 献

[1] BUDAK E, ALTINTAS Y.Analytical prediction of chatter stability in milling-Part I: General formulation[J].Journal of Dynamic Systems Measurement and Control-transactions of The Asme, 1998, 120(1): 22-30.

[2] ALTINTAS Y, ENGIN S.Generalized modeling of mechanics and dynamics of milling cutters[J].CIRP Annals-Manufacturing Technology, 2001, 50(1): 25-30.

[3] ENGIN S, ALTINTAS Y.Mechanics and dynamics of general milling cutters.: Part I: helical end mills[J].International Journal of Machine Tools and Manufacture, 2001, 41(15): 2195-2212.

[4] BALACHANDRAN B, ZHAO M X.A mechanics based model for study of dynamics of milling operations[J].Meccanica, 2000, 35(2): 89-109.

[5] ZHAO M X, BALACHANDRAN B.Dynamics and stability of milling process[J].International Journal of Solids & Structures, 2001, 38(10): 2233-2248.

[6] MERDOL S D, ALTINTAS Y.Multi frequency solution of chatter stability for low immersion milling[J].Journal of Manufacturing Science & Engineering, 2004, 126(3): 459-466.

[7] INSPERGER T, MUNOA J, ZATARAIN M A, et al.Unstable islands in the stability chart of milling processes due to the helix angle[J].CIRP 2nd International Conference on High Performance Cutting, 2006: 12-13.

[9] SCHMITZ T L, SMITH K S.Machining dynamics: frequency response to improved productivity [M].Springer US, 2008.

[10] SIMS N D, MANN B, HUYANAN S.Analytical prediction of chatter stability for variable pitch and variable helix milling tools[J].Journal of Sound & Vibration, 2008, 317(3/4/5): 664-686.

[11] DOMBOVARI Z, STEPAN G.The effect of harmonic helix angle variation on milling stability[J].Journal of Manufacturing Science & Engineering, 2011, 134(5): 921-922.

[12] 丁烨.铣削动力学——稳定性分析方法与应用 [D]上海：上海交通大学, 2011.

[13] 刘强, 李忠群.数控铣削加工过程仿真与优化-建模、算法与工程应用 [M].北京: 航空工业出版社, 2011.

[14] 宋戈, 李剑峰, 孙杰.基于铣削力精确建模的工件表面让刀误差预测分析[J].机械工程学报, 2013, 49(21): 168-175.

SONG Ge, LI Jianfeng, SUN Jie.Analysis on prediction of surface error based on precision milling cutting force model[J].Journal of Mechanical Engineering, 2013, 49(21): 168-175.

[15] 蒋宇平, 龙新华, 孟光.薄壁结构件铣削加工振动稳定性分析[J].振动与冲击, 2016, 35(2): 45-50.

JIANG Yuping, LONG Xinhua, MENG Guang.Stability analysis for thin-walled milling processes[J].Journal of Vibration and Shock, 2016, 35(2): 45-50.