径向浸入率通用型铣削稳定性谱元分析法*

谢建宏，徐菁菁,b，张 宇，李 鸣

(南昌大学 a.信息工程学院；前湖学院，南昌 330031)

0 引言

铣削稳定性是指铣削系统和切削过程相互作用构成的整体系统的稳定性，其实质是描述刀具与工件的相对振动[1]。铣削过程失稳即发生铣削颤振，刀具与工件间产生大幅振动甚至脱离，使得工件加工精度和表面质量超差并显著加剧刀具磨损，同时伴随工作环境极度恶化，甚至造成机床损坏和生产安全事故。保守地选择切削参数虽然能避免发生颤振，却极大地降低了加工效率。因此，对铣削稳定性的准确、高效分析是保证铣削加工质量和效率的关键。考虑切厚再生效应的铣削过程动力学模型较好地描述了颤振机理并得到了大量实验验证[2]。该模型将铣削过程表征为一组周期系数时滞微分方程，其特征方程是超越方程,这给稳定性分析带来了本质性困难[3]。Altintas等开创性地提出了基于Lyapnov判据的单频率法,由此构建了临界切削参数表征的铣削稳定性叶瓣图(Stability Lobe Diagrams, SLDs)以区分稳定切削参数域与不稳定切削参数域[4]。Insperger T等提出了基于Floquet判据的半离散法，具有比单频率法更好的准确性和通用性，成为当前的主流分析方法[5]。此类分析方法的计算瓶颈在于微分方程离散状态向量映射矩阵的获取。针对此问题，国内外研究者提出了许多近似计算方法以提高计算效率从而快速构建SLDs的叶瓣边界[6-9]。基于谱方法和时间有限元法相结合的谱元法给出了一种效率显高于传统半离散法的分析方法[10]，然而大量工程实践发现经典谱元法对于小径向浸入率铣削工况适用性较差，其原因主要在于刀齿切入切出引起的铣削过程强不连续性降低了谱元的收敛性，从而极大地影响了计算效率[11]。

本文提出一种径向浸入率通用型谱元法，将铣削过程时滞周期划分为切削振动和自由振动两个时间段分别进行处理。利用谱元法与微分方程数值解法相结合解决小径向浸入率铣削的强不连续性问题，提高状态向量离散映射矩阵的计算效率，进而保证铣削稳定性分析的准确性、高效性和通用性。

1 铣削稳定性谱元分析方法原理

1.1 铣削过程动力学状态空间模型

经典铣削稳定性分析采用基于再生型动态铣削力模型和铣削工艺系统动力学特性建立的二自由度周期系数时滞微分方程组描述铣削过程动力学，其状态空间模型如式(1)所示：

(1)

式中，τ为时滞周期，在N齿等齿矩铣刀以主轴转速Ω[rpm]铣削的工况下取为一个刀通周期τ=T=60/NΩ，式中状态变量和系数矩阵表达式分别为：

(2)

状态向量u(t)为X,Y方向位移与速度构成的列向量，M,K,C分别表示铣削系统质量、刚度、阻尼矩阵，H(t)表示动态铣削力系数矩阵，如式(3)所示：

(3)

式中，a为轴向切深，Kt为切向切削力系数，Kr为径向与切向切削力系数之比，s表示sin(2φj)，c表示cos(2φj)。瞬时浸入角φj=(2πΩ/60)t- (j- 1)2π/N，j=1,…,N。g表示阶跃函数：

(4)

切入角φst和切出角φex由径向浸入率决定，g(φj)反映了刀齿是否正参与切削。可见，式(1)中，A是常系数矩阵，B(t)是以T为周期的周期系数矩阵。

Floquet稳定性判据依据铣削过程动力学模型在一个周期上的离散状态映射矩阵的特征乘子在Z平面的分布进行铣削过程判稳：若特征乘子都位于单位圆内则铣削过程稳定，否则发生颤振。

1.2 铣削稳定性分析的经典谱元法

经典谱元法将[-τ,τ]时间线划分为2M个时域谱元，每个谱元长h=τ/M，如图1所示。

图1 经典谱元法时间线划分

(5)

(6)

(7)

式(6)回代式(1)，采用加权余量法即可导出一个周期上离散状态向量的映射关系，进而计算离散状态向量映射矩阵特征乘子并根据Floquet判据判别铣削稳定性。由于采用了LGL节点高阶谱插值，使得谱元法相较采用等距节点插值和均值近似的传统半离散法具有更高的计算效率。然而对于小径向浸入率铣削，由于刀齿切入和切出导致H(t)具有强不连续性，使得加权余量法计算过程中LGL数值积分的收敛速度和精度都明显下降，影响了经典谱元法对小径向浸入率铣削的适用性。

2 径向浸入率通用型铣削稳定性谱元分析法

消除刀齿切入切出不连续性影响的一种思路是将图1按铣削过程划分为切削振动时段tc和自由振动时段tf分别求解式(1)，如图2所示。基于此思路提出一种径向浸入率通用型铣削稳定性谱元分析方法。

图2 径向浸入率通用型谱元法时间线划分

式(1)解可描述为积分形式：

(8)

式中，u(t0)为初始时刻t0的状态向量值。不失一般性，设t0为刀具切出时刻，在自由振动区间[t0,t0+tf]不存在切厚再生效应，式(8)中B(η)项为0，退化为无时滞常系数微分方程齐次解形式，得式(9)：

u(t)=eA(t-t0)u(t0)

(9)

解得刀具切入时刻状态向量为：

u(t0+tf)=eA(tf)u(t0)

(10)

在切削振动区间[t0+tf,t0+τ]，存在切厚再生效应使B(η)项不为0，在此区间及上一周期区间 [t0+tf-τ,t0]上分别以时间段h=tc/M等距划分谱单元并将时间坐标t∈[(k-1)h,kh]规范化为[-1,1]区间局部坐标ξ，采用式(6)插值近似，导出kth谱单元上方程余量为：

(11)

(12)

式中，

(13)

(14)

权函数ψi(ξ)取Legendre正交多项式函数族：

(15)

采用LGL数值积分计算式(13)、式(14)得：

(16)

(17)

式中，I表示4×4单位阵，lj′(ξ)表示重心Lagrange插值函数对ξ的一阶导。wj为LGL数值积分系数，如式(18)，其中Ln表示n阶Legendre多项式。

(18)

可见，n+1点LGL数值积分具有2n+1阶代数精度，相较Newton-Cotes积分极大提高了计算精度，且避免了传统等距节点高阶插值的Runge现象导致系统病态无法实现谱收敛的问题。

综合考虑自由振动时段和切削振动时段，将式(10)和式(12)相结合，即可导出铣削过程一个时滞周期上离散状态向量的映射关系：

NUτ=PU0

(19)

(20)

(21)

(22)

则式(1)在一个时滞周期上的离散状态向量映射矩阵Φ为：

Φ=N-1P

(23)

求解矩阵Φ的特征值，根据Floquet稳定性判据即可分析铣削过程稳定性，并据此绘制铣削稳定性叶瓣图。

3 仿真验证

3.1 仿真设置

为验证本文所提算法的准确性、高效性和通用性，针对表1所示相同工艺参数、结构参数、材料参数的二自由度铣削过程动力学模型分别使用径向浸入率通用型谱元法和经典谱元法进行铣削稳定性分析，不失一般性地假设铣削工艺系统动态特性在X,Y方向对称且互不耦合，即式(2)中M,M-1K,M-1C分别退化为对角阵mI,ωn2I, 2ζωnI, 其中m,ωn,ζ分别表示铣削工艺系统的模态质量、固有频率和阻尼比。所有算例均在HP-ZBOOK15G3移动工作站(Core i7 6700hq处理器, 32 GB内存, QUADRO M1000M GPU)上完成。

3.2 仿真结果与分析

图3展示了径向浸入率取1、 0.2、0.05时，分别采用经典谱元法和径向浸入率通用型谱元法进行铣削稳定性分析获取的铣削稳定性叶瓣图，表2给出了两种分析方法仿真铣削稳定性叶瓣图的计算时间对比。可见，在最大分析误差不超过1%的前提下，径向浸入率通用型谱元法在大浸入率工况下具有与经典谱元法相同的计算效率，但在小浸入率工况下计算效率则明显高于经典谱元法，这是由于其较好地处理了铣削力方向系数的不连续性，使得采用较低的谱元阶次即可实现SLDs曲线收敛到准确值。图4和图5分别展示了径向浸入率为0.1的顺铣过程稳定性叶瓣图的经典谱元法和径向浸入率通用型谱元法的收敛情况。对比结果表明对于小径向浸入率铣削，经典谱元法收敛效率明显降低，当谱元阶次取60时叶瓣图基本收敛，而径向浸入率通用型谱元法在谱元阶次取25时已经收敛到准确值，且谱元插值计算仅需在切削振动段进行，进一步缩小了计算区间，极大提高了计算效率，使得算法具有更强的工况通用性。

表1 铣削稳定性仿真算例参数

表2 两种铣削稳定性分析方法的计算效率对比

(a) 径向浸入率= 1经典谱元法 (b) 径向浸入率= 1径向浸入率通用型谱元法

(c) 径向浸入率= 0.2经典谱元法 (d) 径向浸入率= 0.2径向浸入率通用型谱元法

(e) 径向浸入率= 0.05经典谱元法 (f) 径向浸入率= 0.05径向浸入率通用型谱元法图3 经典谱元法和径向浸入率通用型谱元法铣削稳定性叶瓣图仿真结果对比

图4 经典谱元法稳定性叶瓣图收敛情况

图5 径向浸入率通用型谱元法稳定性叶瓣图收敛情况

4 结论

本文阐述了一种用于铣削过程稳定性分析的径向浸入率通用型谱元法，并针对经典铣削工况进行了稳定性叶瓣图仿真，结果表明该方法通过将LGL节点高阶谱元与微分方程数值解法相结合，较好解决了小径向浸入率铣削强不连续性导致的谱元收敛效率降低的问题。在分析误差不超过1%的前提下，相较经典谱元法显著提高了计算效率，实现了大/小径向浸入率铣削工况稳定性分析的强通用性，使无颤振铣削参数在线优化成为可能。