在模具、电子、航空、航天等领域，常常要求高效、高精度地从毛坯中铣削去除大量材料.实践表明，颤振是影响高效、高精度铣削的重要因素之一^[1].颤振使得刀具与工件的相对位置发生改变，降低铣削效率和加工表面的精度，还会加剧刀具的磨损甚至破坏，缩短机床和刀具的使用寿命.因此，铣削过程稳定性分析和铣削颤振的控制理论研究已经成为提高铣削加工效率的重要手段之一，具有重要的研究价值.

铣削稳定性分析是合理选取主轴转速和切深等加工参数以达到优化工艺参数，避免发生颤振，提高生产效率的重要途径.铣削过程中多齿刀具旋转而引起的切削力周期变化和切削过程间断等因素，使得铣削稳定性分析比连续正交切削稳定性分析困难得多.

目前国内外学者针对这类铣削稳定性问题，已经提出许多近似计算方法预测由临界加工参数构成的铣削过程稳定性边界曲线即稳定性叶瓣图(stability lobe diagram，SLD)^[2].文献[3]提出了用于判断延时系统稳定性的半离散法(semi - discretization method，SDM).文献[4]针对单自由度断续切削的情况，提出了时间有限元分析法(temporal finite element analysis，TFEA).文献[5]基于铣削过程动态响应的直接积分格式提出了适用于大(小)径向切深、大(小)轴向切深等多种工况半解析的全离散法(full-discretization method，FDM).文献[6]证明，全离散法与零阶半离散法具有同阶局部截断误差.

以上的铣削稳定性预测方法，不论是半离散法、时间有限元法还是全离散法都是将铣削延时微分方程的时间周期平均分割.本文提出的切比雪夫分割方法是将时间周期T不均匀分割成m个时间元，刀具任意时刻的状态项和时间周期项用时间元端点的状态项线性表示，延时项采用二阶Newton差商的方法用时间元端点的状态项插值表示.首先求出铣削系统在每个时间元上的传递矩阵，再求铣削系统在一个时间周期内的传递矩阵，再用Floquet理论^[7]来判断铣削系统的稳定性.通过时域内仿真结果可以看出：在相同精度要求下，切比雪夫不平均分割时间周期方法的计算效率比平均分割时间周期方法的计算效率高.

逼近法近似求函数解析表达式是在函数定义域区间段上取n+1个互异的节点，给定这些节点的函数值，采用Lagrange或者Newton插值方法^[8]构造出一个唯一的n次插值多项式来逼近原函数.采用均匀取点法构造出来的n+1个节点的插值型积分公式一般具有n次代数精度.如果适当地选取节点来构造求积公式，就可以得到代数精度更高的求积公式.切比雪夫分割法^[9]就是这样一种取点分割方法，它的取点分割过程如下.

如图 1所示，将半圆弧按照角θ分割，圆弧上的分割点在x轴上的投影点即是该区间段[-1,1]上的分割点.

图1(Fig. 1)

图1 切比雪夫分割图(N0=4) Fig. 1 Image of Chebyshev segmentation

在区间[a,b]上，x轴上取点的位置表达式为

只考虑再生颤振的铣削延时微分方程为

切比雪夫分割方法将时间周期T不均匀分割为m个时间元：

在第一个时间元内t_k≤t≤t_k+1(k=0，1，…，m-1)，式(4)的响应可以写成直接积分格式：

X(t_k)表示在t_k时刻刀具的状态向量^[6]，则在t_k+1时刻刀具的状态向量X(t_k+1)为

此时式(10)中的积分时间区间为[0,Δt_k].对周期系数矩阵B(t_k+1-ξ)，X(t_k+1-ξ)在时间元t_k≤t≤t_k+1上做线性逼近^{[10, 11, 12]}；对于延时项X(t_k+1-ξ-T)可以用时间元[t_k-T,t_k+1-T]两端点时刻铣刀的状态向量X(t_k+1-T)和X(t_k-T)采用二阶Newton差商的方法线性表示，见式(11)；将式(11)带入式(10)得式(12).

式(12)中的变量F_k+1，F_k-1，F_0,k分别为

设[I-F_k+1]^-1存在，则铣削系统每个时间元上的传递矩阵Ψ_k可以表示为

铣削系统在一个时间周期T上的传递矩阵Ψ可以用这个周期内时间元的传递矩阵序列Ψ_k (k=0，1，…，m-1)表示，即

最后，根据Floquet理论，铣削系统的稳定性可由系统传递矩阵Ψ的特征值决定：若传递矩阵Ψ的所有特征值的模长均小于1，则系统稳定；否则系统不稳定.

两自由度铣削系统如下^[12]：两个刀齿的槽式直铣刀；铣刀齿数N=2；质量m_tx=0.03993kg，m_ty=0.03993kg；阻尼比ζ=0.011；系统的固有频率f_n=ω_n/(2π)=922Hz；铣刀转速Ω变化范围为0.5×10⁴~2.5×10⁴r/min；切深w的变化范围为0~0.01m；K_t=6×10⁸N/m²；K_n=2×10⁸N/m²；铣削系统特性的常数矩阵A：

再生效应动态铣削力决定的周期矩阵B：

K_t是线性化切向切削力系数，K_n是线性化法向切削力系数，a/D表示铣刀的径向侵入比.

将转速-轴向切深平面平均划分为200×100网格点，即转速步长100r/min，轴向切深步长为0.1mm.分别计算径向侵入比等于0.1和0.2两种工况下的铣削稳定性极限.根据前面所述的理论基础和算法原理，编写计算程序，利用计算机绘制稳定性叶瓣图,如表 1所示.

表1(Table 1)

表1 两自由度铣削模型的算法比较 Table 1 Algorithm comparison of two degrees of freedom milling model 切比雪夫分割法下的全离散(m=30) 平均分割法下的全离散法(m=40) a/D=0.1 计算时间/s 242.819 501.195 a/D=0.2 计算时间/s 251.978 509.060

表1 两自由度铣削模型的算法比较 Table 1 Algorithm comparison of two degrees of freedom milling model

由表 1可以看出在用叶瓣图来预测铣削系统的稳定性时，切比雪夫分割法取的点数比平均分割法取的点数少的情况下，叶瓣图仍能达到相同的精度.由两种分割方法的计算时间的对比可知，获得相同精度叶瓣图时，切比雪夫分割法的计算效率比平均分割法计算效率高.

图 2表示a/D=0.1，m =40时切比雪夫分割法和平均分割法得到的稳定性叶瓣对比图，这两种分割方法得到的铣削系统的稳定性范围有所不同.在图 2中取4个工作点：点1，点2，点3，点4，这4个工作点对应的铣削参数和两种分割方法预测的4个工作点的稳定性如表 2所示.对表 2中所列工作点进行时域仿真，仿真结果如图 3所示.工作点的仿真结果如图 3所示.由图 3可知，工作点1，点2，点3对应的铣削过程是稳定的，点4对应的铣削过程是不稳定的.图 3表明4个工作点的时域仿真结果与表 2中切比雪夫分割法预测的4个工作点的稳定性完全吻合，而与平均分割法预测的4个工作点的稳定性相差较多.因此，取相同分割点时，切比雪夫分割法对铣削系统的稳定性预测结果比平均分割法对铣削系统的稳定性预测结果更准确.

图2(Fig. 2)

图2 两种分割方法下的稳定性叶瓣图对比(顺铣N=2，a/D=0.1 ) Fig. 2 Comparison of two different segmentation’s stability lobes diagram

表2(Table 2)

表2 两种取点分割法得到的部分工作点的稳定性 Table 2 Stability of sectional working points for two way to take segmentation 工作点 Ω w 切比雪 夫分割法 平均分 割法 r·min-1 mm 点1 15000 1.5 稳定 不稳定 点2 20000 1 稳定 稳定 点3 20000 2 稳定 不稳定 点4 20000 3 不稳定 不稳定

表2 两种取点分割法得到的部分工作点的稳定性 Table 2 Stability of sectional working points for two way to take segmentation

图3(Fig. 3)

图3 部分工作点的时域仿真结果图 Fig. 3 Time domain simulation’s results for part of working points (a)—点1； (b)—点2; (c)—点3； (d)—点4.

本文介绍了切比雪夫分割方法和Newton二阶差商法对离散项进行插值积分的方法.分析两自由度铣削系统示例，由绘制的系统稳定性叶瓣图可知：在达到相同预测精度的条件下，切比雪夫不平均分割法需要取的点数比平均分割法所需的点数少，并且切比雪夫不平均分割法的计算效率比平均分割法计算效率高；在分割点数相同的 情况下，切比雪夫不平均分割法预测的系统稳定性比平均分割法预测的系统稳定性更加准确.