薄板件在航空航天、能源动力等领域应用广泛，但因刚度低、材料加工性能差、材料去除率高等特点，在铣削加工过程中极易发生颤振^[1]。颤振会影响工件表面的加工质量，加剧刀具的磨损。因此，为保证薄板件的加工质量，铣削稳定性研究势在必行。

目前，判别切削过程稳定性主要通过稳定性叶瓣图(stability lobe diagram, SLD)，一般而言SLD曲线下方为稳定区域，上方为颤振区域。SLD的求解方法主要有频域法和时域法2类^[2]。频域法又分为零阶频域法^[3]和多阶频域法^[4]。相比零阶频域法，多阶频域法在径向切深率较小时预测精度更高，所以更适合薄板件的加工。但是频域法忽略了非线性因素对于稳定性的影响，而时域法考虑了刀齿回转一周过程中非线性切削力的影响，所以预测精度更高。时域法又可分为半离散法^[5-6]和全离散法^[7]。全离散法在半离散法的基础上，对延时微分方程中的周期部分离散化，提高了计算效率。

铣削动力学模型直接影响SLD的预测结果。对于弱刚性薄板件的铣削加工过程来说，需要同时考虑刀具和工件的动态特性^[8]。刀具的动态特性通常是指刀尖点处的频响函数，可通过锤击实验法得到。工件的动态特性通常是指薄板件被切削部位的频响函数。在铣削过程中，由于刀具的进给，被切削部位不断变化，使得各点处的SLD预测结果均不相同，传统的二维SLD无法反映各不同点的稳定性预测结果。Jin等^[9]提出了三维SLD，考虑了工件的不同被切削位置对SLD的影响。对于薄板件粗加工阶段，当材料去除率相对于薄板件壁厚无法忽略时，还应考虑不同切削阶段薄壁件的模态参数^[10]，即考虑工件模态的时变特性。李长城等^[11]对工件切削的不同阶段进行分布建模，分别测量其频响函数并绘制SLD。另外，工件受瞬时切削力的动态变形量也会影响SLD，Li等^[12]提出了一种柔性工件变形预测模型，并分析了工件变形对铣削稳定性的影响。周凯^[13]和Zatarain等^[14]发现在SLD中的主轴高转速段，在稳定域下方会出现附加的不稳定“孤岛”区域，这主要是由于铣刀螺旋角所导致的，但在主轴低转速段并未发现不稳定“孤岛”现象。

根据薄壁板铣削系统的动力学特性，针对薄板件精加工阶段的颤振问题，本文建立了三维多自由度铣削动力学模型，考虑了刀具X、Y方向和工件Z方向的多阶模态，同时考虑了前一次切削遗留的圆弧形刀痕对径向切深率的影响；根据平均切削力法辨识得到的切削力系数和锤击实验法得到的刀具、工件模态参数，利用全离散法得到工件不同点处的SLD，并分析了工件模态刚度、模态频率、模态阻尼比对SLD中不稳定“孤岛”的影响；最后通过多组切削实验验证了预测结果。

1 铣削稳定性预测模型 1.1 刀工接触区域修正

在建立动态铣削力模型时，考虑铣刀在X和Y方向的自由度以及工件在Z方向的自由度(见图 1)。对于球头铣刀来说，由于螺旋角和球头半径的影响，刀齿的径向接触角不仅与时间t有关，还与轴向位置有关，如图 2所示，其中P为刀具切削点。假设球头刀螺旋线导程不变，轴向滞后角与轴向位置z呈线性关系，刀齿j在轴向位置z处的径向接触角φ_jl(t)可表示为

$ {\varphi _{jl}}(t) = \frac{{2{\rm{ \mathsf{ π} }}\mathit{\Omega }}}{{60}}t + \frac{{2{\rm{ \mathsf{ π} }}}}{N} - \frac{{z\tan \beta }}{R}. $

(1)

其中：Ω为主轴转速，N为刀齿数，β为铣刀柄部螺旋角，R为铣刀刀柄处半径，z为轴向位置。

$ z=l d z. $

其中：l为当前轴向位置微元的编号，dz为当前轴向位置微元的厚度。

轴向接触角k(z)为

$ k(z)=\left\{\begin{array}{ll} \arccos \left(\frac{R-z}{R}\right), & 0 \leqslant z \leqslant R ; \\ \frac{\pi}{2}, & z>R . \end{array}\right. $

(2)

判断函数g(z)用于判断刀齿是否处于切削中：

$ g(z)=\left\{\begin{array}{ll} 1, & \varphi_{\mathrm{st}} <\varphi_{i l}(t)<\varphi_{\mathrm{ex}} ; \\ 0, & \text { 其他. } \end{array}\right. $

(3)

其中：φ_st为刀具切入时的径向接触角，φ_ex为刀具切出时的径向接触角，1和0分别表示刀齿处于切削状态和非切削状态。当采用顺铣时，有：

$ \left\{\begin{array}{l} \varphi_{\mathrm{st}}=\pi-\arccos (1-2 \varepsilon(z)); \\ \varphi_{\mathrm{ex}}=\pi. \end{array}\right. $

(4)

其中ε(z)表示径向切深率。

图 3中，对球头刀而言，径向切深率与前一次切削时留下的圆弧形刀痕有关，而圆弧形刀痕与铣刀半径R有关：

$ \varepsilon(z)=\left\{\begin{array}{ll} 1, & 0 \leqslant z <R-\sqrt{R^{2}-\frac{\left(a_{\mathrm{e}}^{0}\right)^{2}}{4}} ; \\ \frac{a_{\mathrm{e}}^{0}}{2 \sqrt{R^{2}-(R-z)^{2}}}, & R-\sqrt{R^{2}-\frac{\left(a_{\mathrm{e}}^{0}\right)^{2}}{4}} \leqslant z<R ; \\ \frac{a_{\mathrm{e}}^{0}, }{2 R}, & z \geqslant R . \end{array}\right. $

(5)

其中a_e⁰表示刀柄处的径向切深。

1.2 三维多自由度铣削动力学模型

三维瞬时动态铣削厚度为^[15]

$ \begin{array}{c} h_{\mathrm{d}, j l}(t)=\left[\left(x_{\mathrm{t}, j l}(t)-x_{\mathrm{t}, j l}(t-T)\right) \sin \varphi_{j l} \sin k(z)+\right. \\ \quad\left(y_{\mathrm{t}, j l}(t)-y_{\mathrm{t}, j l}(t-T)\right) \cos \varphi_{j l} \sin k(z)- \\ \left.\left(z_{\mathrm{w}, j l}(t)-z_{\mathrm{w}, j l}(t-T)\right) \cos k(z)\right] g(z) . \end{array} $

(6)

其中：x_{t, jl}(t)、y_{t, jl}(t)和x_{t, jl}(t－T)、y_{t, jl}(t－T)分别表示铣刀在当前切削周期和前一切削周期的动态位移，z_{w, jl}(t)和z_{w, jl}(t－T)分别表示工件在当前切削周期和前一切削周期的动态位移，T表示延时周期。

刀齿j在轴向位置z的三方向动态微元铣削力为

$ \begin{array}{c} \left[\begin{array}{c} F_{x, j l} \\ F_{y, j l} \\ F_{z, j l} \end{array}\right]=\\ \left[\begin{array}{ccc} -\cos \varphi_{j l} & -\sin \varphi_{j l} \sin k(z) & -\sin \varphi_{j l} \cos k(z) \\ \sin \varphi_{j l} & -\cos \varphi_{j l} \sin k(z) & -\cos \varphi_{j l} \cos k(z) \\ 0 & \cos k(z) & -\sin k(z) \end{array}\right] .\\ \left[\begin{array}{c} K_{\mathrm{tc}} \\ K_{\mathrm{rc}} \\ K_{\mathrm{ac}} \end{array}\right]\left[\begin{array}{c} \sin \varphi_{j l} \sin k(z) \\ \cos \varphi_{j l} \sin k(z) \\ \cos k(z) \end{array}\right]^{\mathrm{T}} \frac{\mathrm{d} z}{\sin k(z)} g(z) \cdot\\ \left[\begin{array}{c} x_{\mathrm{t}, j l}(t)-x_{\mathrm{t}, j l}(t-T) \\ y_{\mathrm{t}, j l}(t)-y_{\mathrm{t}, j l}(t-T) \\ z_{\mathrm{w}, j l}(t)-z_{\mathrm{w}, j l}(t-T) \end{array}\right] . \end{array} $

(7)

其中：K_tc、K_rc和K_ac分别表示当前刀具工件之间的切向、径向和轴向切削力系数。

对式(7)先按轴向位置z积分，再按刀齿j积分即可得到三方向动态铣削力。

考虑刀具在X、Y方向以及工件在Z方向上各阶模态，建立图 4的铣削动力学模型^[13]。

铣刀和工件在X、Y、Z方向上的动态位移可由各阶模态对应的动态位移表示。根据集中质量法，三维多自由度系统的动力学方程可以表示为

$ \begin{array}{c} \mathit{\boldsymbol{M}}\left[ {\begin{array}{*{20}{c}} {{{\ddot x}_{{\rm{t}}, 1}}(t)}\\ {{{\ddot x}_{{\rm{t}}, 2}}(t)}\\ \vdots \\ {{{\ddot x}_{{\rm{t}}, {N_{\rm{m}}}}}(t)}\\ {{{\ddot y}_{{\rm{t}}, 1}}(t)}\\ {{{\ddot y}_{{\rm{t}}, 2}}(t)}\\ \vdots \\ {{{\ddot y}_{{\rm{t}}, {N_{\rm{m}}}}}(t)}\\ {{{\ddot z}_{{\rm{w}} \cdot 1}}(t)}\\ {{{\ddot z}_{{\rm{w}}, 2}}(t)}\\ \vdots \\ {{{\ddot z}_{{\rm{w}}, {N_{\rm{m}}}}}(t)} \end{array}} \right] + \mathit{\boldsymbol{C}}\left[ {\begin{array}{*{20}{c}} {{{\dot x}_{{\rm{t}}, 1}}(t)}\\ {{{\dot x}_{{\rm{t}}, 2}}(t)}\\ \vdots \\ {{{\dot x}_{{\rm{t}}, {N_{\rm{m}}}}}(t)}\\ {{{\dot y}_{{\rm{t}}, 1}}(t)}\\ {{{\dot y}_{{\rm{t}}, 2}}(t)}\\ \vdots \\ {{{\dot y}_{{\rm{t}}, {N_{\rm{m}}}}}(t)}\\ {{{\dot z}_{{\rm{w}}, 1}}(t)}\\ {{{\dot z}_{{\rm{w}}, 2}}(t)}\\ \vdots \\ {{{\dot z}_{{\rm{w}}, {N_{\rm{m}}}}}(t)} \end{array}} \right] + \\ \mathit{\boldsymbol{K}}\left[ {\begin{array}{*{20}{c}} {{x_{{\rm{t}}, 1}}(t)}\\ {{x_{{\rm{t}}, 2}}(t)}\\ \vdots \\ {{x_{{\rm{t}}, {N_{\rm{m}}}}}(t)}\\ {{y_{{\rm{t}}, 1}}(t)}\\ {{y_{{\rm{t}}, 2}}(t)}\\ \vdots \\ {{y_{{\rm{t}}, {N_{\rm{m}}}}}(t)}\\ {{z_{{\rm{w}}, 1}}(t)}\\ {{z_{{\rm{w}}, 2}}(t)}\\ \vdots \\ {{z_{{\rm{w}}, {N_{\rm{m}}}}}(t)} \end{array}} \right] = \mathit{\boldsymbol{G}}\left[ {\begin{array}{*{20}{c}} {{x_{{\rm{t}}, 1}}(t) - {x_{{\rm{t}}, 1}}(t - T)}\\ {{x_{{\rm{t}}, 2}}(t) - {x_{{\rm{t}}, 2}}(t - T)}\\ \vdots \\ {{x_{{\rm{t}}, {N_{\rm{m}}}}}(t) - {x_{{\rm{t}}, {N_{\rm{m}}}}}(t - T)}\\ {{y_{{\rm{t}}, 1}}(t) - {y_{{\rm{t}}, 1}}(t - T)}\\ {{y_{{\rm{t}}, 2}}(t) - {y_{{\rm{t}}, 2}}(t - T)}\\ \vdots \\ {{y_{{\rm{t}}, {N_{\rm{m}}}}}(t) - {y_{{\rm{t}}, {N_{\rm{m}}}}}(t - T)}\\ {{z_{{\rm{w}}, 1}}(t) - {z_{{\rm{w}}, 1}}(t - T)}\\ {{z_{{\rm{w}}, 2}}(t) - {z_{{\rm{w}}, 2}}(t - T)}\\ \vdots \\ {{z_{{\rm{w}}, {N_{\rm{m}}}}}(t) - {z_{{\rm{w}}, {N_{\rm{m}}}}}(t - T)} \end{array}} \right]. \end{array} $

(8)

其中：M表示模态质量矩阵，C表示模态阻尼矩阵，K表示模态刚度矩阵，G表示方向系数矩阵。

将式(8)进行变量替换，转换为延时微分方程，并在不同加工参数下利用全离散法^[13]求解延时微分方程，即可讨论铣削稳定性，从而实现稳定性预测。

2 铣削稳定性预测结果 2.1 切削力系数辨识

本文采用Altintas^[16]提出的平均切削力法辨识切削力系数。实验所用测力仪为Kistler 9255C型三向测力仪，机床为科德KMC600 U型机床，刀具为锥度球头铣刀(球头半径5 mm，锥度4°，3齿，螺旋角35°)，工件为TC4钛合金。实验采用干切削、全齿铣削的方式，结果为：K_tc=2 234.2 MPa，K_rc=623.01 MPa，K_ac=1 200.5 MPa。

2.2 刀具工件模态实验结果

考虑到薄板件铣削系统的动力学特性，仅考虑工件Z方向以及刀具X、Y方向的模态。通过锤击实验法得到刀具和工件的模态如表 1和2所示，均取前3阶模态。所用力锤型号为Kistler9722A500，灵敏度为10 mV/N。传感器为Kistler8778A500型加速度传感器，灵敏度为10.79 mV/g，重量为0.4 g，所以可忽略传感器质量对于刀具、工件模态的影响。利用NI信号采集器(cDAQ-9171型信号采集器，最高采样频率51 200 Hz)采集信号，并用Cutpro软件分析拟合模态参数。由于薄板件各点处动力学响应均不同，因此将工件表面的切削路径离散化为11个测量点(见图 5中的红点)，分别测量其模态参数。考虑到实验板关于X轴对称，表 2仅列出测量点1~6的模态结果。

2.3 铣削稳定性预测结果

根据实验得到的切削力系数和刀具、工件的模态参数，结合三维多自由度铣削动力学模型，利用全离散法可得到离散化各测量点处的SLD，如图 6所示，a_p为轴向切深。图中紫色曲线为预测得到的稳定域分界线，和为实验结果(实验见3.2节)，分别为稳定切削和颤振。由于刀具单次走刀时实验板关于X轴近似呈轴对称，因此仅绘制了测量点1至6的SLD。

因为薄板工件沿切削路径各处的动力学特性不同，所以当球头铣刀沿着图 5所示的Y方向进给过程中，工件的动态响应也不同，导致了工件各处SLD的预测结果随着切削路径的变化而变化。由于实验用薄板工件的刚度较小，在SLD中出现了多个封闭的不稳定“孤岛”区域。并且随着工件模态刚度逐渐增大，这些不稳定“孤岛”区域逐渐变小，甚至某些不稳定“孤岛”区域渐渐变为了稳定区域(对比图 6b和6c)。

2.4 稳定性叶瓣图影响因素分析

为进一步分析不稳定“孤岛”现象产生的原因，本节采用控制变量的方法仿真分析各项参数对SLD的影响。

根据表 2发现工件第3阶模态相比其他各阶模态在不同测量点处差异显著，该阶模态对应的模态参数中模态刚度、模态质量在工件不同位置处的差异尤为明显。而对于某阶模态来说，其模态频率f_m、模态刚度k_m、模态质量m_m满足如下关系式：

$ f_{\mathrm{m}}=\frac{1}{2 \pi} \sqrt{\frac{k_{\mathrm{m}}}{m_{\mathrm{m}}}}. $

(9)

即当模态频率不变时，模态刚度和模态质量是2个相互影响的变量，在分析二者对于SLD的影响时仅需考虑其中之一即可。

首先分析第3阶模态的模态刚度对SLD中不稳定“孤岛”的影响，仿真时工件模态参数如表 3所示。其中，k_{m, 3}分别为5×10⁶、8×10⁶、11×10⁶ N/m，m_{m, 3}通过式(9)计算得到。仿真结果如图 7所示。由图可得当第3阶模态刚度变大时，主轴转速2 275、2 825 r/min附近的不稳定“孤岛”逐渐变小直至消失，这说明了模态刚度会影响不稳定“孤岛”的大小。而主轴转速2 000 r/min附近的不稳定“孤岛”基本不变，推测是由其他模态导致的。主轴转速3 600 r/min附近的不稳定“孤岛”从其形状推测是由2个不稳定“子孤岛”叠加而成，其中一个不稳定“子孤岛”逐渐变小，是由于第3阶模态刚度变大导致的，而另一个不稳定“子孤岛”无明显变化，推测是由其他模态导致的。

在图 7c仿真参数的基础上，进一步修改第1阶模态的模态刚度，仿真参数如表 4所示。其中，k_{m, 1}分别为2×10⁶、2.5×10⁶、3×10⁶ N/m。仿真结果如图 8，当第1阶模态刚度变大时，主轴转速2 000、3 600 r/min附近的不稳定“孤岛”逐渐变小直至消失。结合图 7、8可得各阶模态仅能影响特定的若干个不稳定“孤岛”，对其他不稳定“孤岛”影响很小。

然后分析工件模态频率对于不稳定“孤岛”的影响。仿真参数如表 5所示。其中，f_{m, 3}分别为550、575、600 Hz。仿真结果如图 9所示。可见当模态频率变大时，该阶模态频率对应的不稳定“孤岛”会出现“右移”现象，但“孤岛”的面积变化很小。

最后分析工件模态阻尼比对不稳定“孤岛”的影响。仿真参数如表 6所示。其中，ζ_{m, 3}分别为0.000 5、0.001、0.001 5。仿真结果如图 10所示。可见随着ζ_{m, 3}的增大，不稳定“孤岛”区域逐渐变小直至消失。

3 实验验证 3.1 实验装置

铣削颤振实验现场如图 11所示。实验板材料为钛合金TC4，大小为350 mm×150 mm×9 mm。所用刀具与切削力系数辨识环节的刀具一致。实验系统包含被切实验板、压电陶瓷作动器(芯明天PSt150/10/60VS15型压电陶瓷，行程57，刚度35，标称推力2 300 N)及其驱动器(芯明天XE501-B型驱动器)、加速度传感器(Kistler8778A500型加速度传感器，灵敏度为10.79 mV/g)、力传感器(斯巴拓SBT641力传感器，量程490 N，灵敏度为2 mV/V)、NI信号采集器(cDAQ-9171型)等。

由于实验板厚度较小，在铣削过程中受力易变形，因此需要在切削实验前先对其变形量进行预测并补偿到刀具路径中以提高实验的准确性。具体方法为首先利用Altintas^[16]提出的平均切削力法结合加工参数预测出各组实验的轴向切削力F_z。接着，利用压电陶瓷和丝杠螺母的组合机构在实验板离散的11个测量点处分别提供一个垂向的支撑力，并通过压电陶瓷底部固定的力传感器测量支撑力的大小。然后利用千分表测量施加支撑力前后各测量点处的位移大小，即可得到各点处的等效静刚度。最后结合轴向切削力F_z得到这11个测量点处的变形量预测值，并通过线性插值的方式得到切削路径上各点的变形量预测值补偿到刀具路径中。

3.2 实验结果及分析

各组颤振实验所使用的加工参数如表 7所示。利用粘贴在实验板下方的加速度传感器测得铣削过程中的加速度信号，并通过Fourier变换得到信号的频谱图。当铣削过程为稳定切削状态时，频谱图中仅会出现主轴旋转频率及其倍频；而当铣削过程为不稳定(颤振)切削状态时，频谱图中除了主轴旋转频率及其倍频之外，还存在其他幅值较大的颤振频率。由于实验板各处的SLD预测结果均不同，因此需要对测得的加速度信号进行分段分析，由此来判断各测量点处是否发生了颤振。7组实验各测量点处的实验结果如图 6所示，图中稳定切削状态用表示，颤振状态用表示。下面，选取具有代表性的2组实验(实验5、7)结果进行具体分析。

实验5的稳定性预测结果显示测量点1、2之间和测量点10、11之间(即实验板边缘)为颤振区域，测量点2、10之间为稳定区域。实验测得的加速度信号频谱图如图 12所示。在测量点1、2之间观察到除了主轴频率(spindle frequency, SF)及其倍频之外，还存在489 Hz的颤振频率(chatter frequency, CF)，测量点10、11之间也存在490 Hz的颤振频率。而测量点2、10之间仅存在SF及其倍频。这与预测结果一致。

实验7的预测结果显示各测量点均为颤振区域。实验测得的信号频谱图如图 13所示，由于实验板关于X轴对称，因此仅分析了测量点1~6的加速度信号。从图中可以看出各频谱图中均存在颤振频率，与预测结果一致。以图 13a为例，各颤振频率与其各自相邻的主轴频率的倍频之差均为15 Hz。由于实验板各处动态响应特性均不同，因此各频谱图中的颤振频率也存在变化。

利用相机在相同光照条件下拍摄实验4(稳定切削)与实验7(剧烈颤振)部分铣削表面图像如图 14所示。以测量点4为例，稳定切削和颤振状态下的表面粗糙度值分别为1.263、2.855。通过对比，可以看出颤振现象会影响实验板表面铣削质量，产生明显的颤振振纹。

4 结论

本文建立了球头铣刀铣削薄板件的三维多自由度铣削动力学模型，并考虑了前一次切削遗留的圆弧形刀痕对刀工接触区域判定的影响，通过锤击实验法测得刀具X、Y方向和工件Z方向的模态参数，并利用平均切削力法拟合得到切削力系数，最后利用全离散法绘制了实验薄板件各离散化测量点处的SLD。

预测结果表明，工件各点处动力学特性的不同会影响铣削稳定性。并且在SLD中存在多个封闭的不稳定“孤岛”区域。不稳定“孤岛”与工件各项模态参数有关，具体为工件各阶模态仅能影响特定若干个不稳定“孤岛”，对其他不稳定“孤岛”影响较小；当模态刚度变大时，该阶模态对应的不稳定“孤岛”会变小直至消失；当模态频率增大时，SLD中不稳定“孤岛”会右移；当模态阻尼比增大时，不稳定“孤岛”变小直至消失。为了验证预测结果，设计了铣削颤振实验，测得的加速度信号频谱图与预测结果一致。并且当颤振较为明显时，频谱图中的颤振频率与其各自相邻的主轴频率的倍频之差为定值。