变工况3维铣削稳定性预报方法研究

2021-07-29郑华林

工程科学与技术 2021年4期

胡腾，王虎，阳红，郑华林，刘雁

(1.西南石油大学机电工程学院，四川成都 610500；2.中国工程物理研究院机械制造工艺研究所，四川绵阳 621999；3.四川普什宁江机床有限公司，四川都江堰 611830)

目前，颤振已成为高效铣削领域亟待解决的关键共性问题[1]。准确预报铣削稳定性并将其作为约束条件，可为铣削工艺参数优化提供坚实的理论基础和数据支撑，有利于实现高效、无颤振铣削加工，亦有利于延长刀具及机床寿命，具有较重要的科学意义与工程价值。

近年来，国内外学者针对铣削稳定性预报方法开展了较深入的研究，取得了不少成果：Ramirez等[2]提出了基于Lyapunov-Krasovskii分析的铣削稳定性预报方法，揭示了稳定性与切削参数之间的复杂联系；Sekar等[3]联合自适应步长R-K方法与误差控制算法，提出了铣削稳定域高精度仿真方法，该方法亦可推广至多时滞动力学特征值问题的求解；Li等[4]提出了一种直接利用加工系统模态参数或频响函数快速绘制稳定性叶瓣图的方法，较显著地提升了切削稳定性预报效率；Ahmadi等[5]分别基于多频率法和半离散法对铣削稳定性进行了预报，并借助等效粘性模型探明了过程阻尼对系统低转速区切削颤振的抑制作用；Li等[6]建立了考虑模态耦合和过程阻尼的铣削动力学模型，并利用二阶半离散法对铣削稳定性进行了预报；Ezugwu等[7]提出了基于最小二乘近似的全离散法，并以此分别对再生型极限轴向切深和极限径向切宽进行了分析与预报。王民等[8]从再生型颤振机理出发，基于零阶频域法构建了3维稳定性叶瓣图，系统探讨了主轴转速、轴向切深与径向切宽之间复杂的非线性映射关系；宋清华等[9]通过考虑刀具偏心建立了计及变径向切宽的铣削动力学模型，并基于半离散法预报了该切削系统的稳定域；李欣等[10]建立了考虑过程阻尼的铣削动力学模型，并基于隐式龙格库塔法绘制了稳定性叶瓣图，进而在低转速区对极限轴向切深进行了较准确预报。

现阶段，铣削稳定性预报研究大多仅针对切削宽度ae恒定的定工况，这并不利于面向最大金属去除率开展铣削工艺参数优化设计。一方面，根据金属去除率（metal removal rate，MRR）定义[11]：

从切削动力学的角度来看， Ω、ap、ae三者间存在着相互制约的复杂关系，若只针对固定切宽ae进行切削稳定性预报，等同于在设计空间中缩小了ap与ae的优化可行域，致使难以获取变工况条件下MRR全局最大值；另一方面，文献[7-9]虽然针对切宽可变条件下的铣削稳定性进行了较系统研究，但在构建铣削动力学模型时鲜有完整考虑切削过程阻尼及刀具动力学行为非对称性等因素，并不利于提升铣削稳定性预报精度。

为解决上述问题，以时域全离散方法为理论基础，综合考虑铣削过程阻尼与铣刀径向动力学行为非对称性，面向可变切削深度及切削宽度的变工况提出一种3维铣削稳定性预报方法。以再生效应所致动态切屑厚度为切入点，系统分析两自由度铣削几何，并引入铣削过程阻尼，建立动态铣削力模型；其次，以曲面铣削时径向切削宽度的时变特征为工程背景，推导切宽及刀齿切入/切出角数学表达，进而结合动态铣削力模型与刀尖点径向非对称动力学特性，构建变工况铣削动力学模型；在此基础上，借助时域全离散法对变工况铣削稳定性进行预报，并绘制3维铣削稳定性叶瓣曲面；最后，以某型立式加工中心为平台开展变工况铣削实验，对所提稳定性预报方法进行验证。

1 变工况铣削动力学建模

1.1 考虑过程阻尼的动态铣削力

不失一般性地，构建如图1所示两自由度铣削动力学模型。图1中，Ω为主轴转速，铣刀刀齿j的瞬时径向接触角为φj。由再生效应[6]可知，瞬时动态切削厚度可表示为：

图1 两自由度铣削动力学模型Fig.1 2-DOF milling dynamic model

式中：ft为每齿进给量；T为刀齿切削周期；x(t)、x(t-T)和y(t)、y(t-T)则分别为t、(t-T)时刻刀齿j在x和y向的振动位移；φj为刀齿瞬时径向接触角，且：

另一方面，作用在刀齿j上的切向切削力Ftj与径向切削力Frj可表示为[5]：

式中，下标s与p分别表示剪切力与过程阻尼所致犁耕力。若轴向切深为ap，则铣削剪切力可表示为：

式中：Kt、Kr分别为切向、径向切削力系数；g为判定函数，且：

式中，φst和φex分别为切入、切出角。

Ahmadi等[12]指出，可将过程阻尼等效为黏性阻尼Ceq，进而径向和切向犁耕力为：

式（7）～（8）中，Ksp、μ、W与v分别表示压入系数、库仑摩擦系数、刀具后刀面磨损带宽度与切削速度。将切向切削力Ftj与径向切削力Frj在刀具坐标系下作进一步正交分解，可最终获得铣刀在x、y方向上所受铣削力为：

式中，hkl、gpkl（k，l=x、y）为方向系数矩阵元素，且：

1.2 基于时变径向切宽的刀齿切入/切出角

通常情况下，铣削加工均采用径向切削宽度恒定的工艺过程。然而，针对如叶片、机匣等具有复杂几何型面的零件进行铣削加工时，径向切宽则具有时变特性，致使刀齿切入/切出角在切削进给期间并非恒定不变，进而切削系统动力学行为亦随之变化。

某凸曲面顺铣工艺过程如图2所示，其中：ae(t)为时变径向切削宽度；D1与D2分别位于待加工表面与已加工表面，且分别对应t时刻刀齿切入与切出点；Ob(t)与Ow(t)则分别为当前时刻上述两点处对应的瞬时曲率中心；XwYwZw与XsYsZs分别为工件坐标系与刀具坐标系，且二坐标系在XY平面内的瞬时夹角为θ(t)；φst(t)与φex(t)则分别刀齿切入角与切出角，且[12]：

图2 凸面铣削过程时变切宽Fig.2 Time-varying cutting width during convex milling process

式中，R为刀具半径。不难看出，曲面切削时，刀齿切入/切出角取决于R与ae(t)。

不妨设待加工表面与已加工表面在XY平面内分别表示为fb与fw，且：

在t时刻，D2在工件坐标系下的坐标（xD2，yD2）可借助式（13）及NC代码获取，则该点切线方程与法线方程可分别写作：

另一方面，D1在工件坐标系下的坐标（xD1，yD1）可联立式（12）与刀具端面包络圆方程求得；进而，可得过D1且与式（14）平行的直线方程：

联立式（15）、（16）即可得二者交点D3坐标（xD3，yD3），则：

至此，便求得时变径向切削宽度ae(t)的一般表达式，进而借助式（11）即可计算得任意时刻刀齿切入/切出角。

1.3 变工况铣削动力学模型

综合第1.1、1.2节所述内容，考虑变径向切宽、过程阻尼和刀尖点动力学行为非对称性的变工况两自由度铣削动力学模型可写为:

2 时域全离散的变工况铣削稳定性预报

利用柯西变换，可得变工况铣削动力学模型（18）的状态空间方程为[13]：

式中，A0和A(t)为周期系数矩阵，Ff(t)为周期系数向量，CP为等效黏性阻尼矩阵，且：

其中，q(t)与p(t)可写为：

基于时域全离散法[14]求解铣削系统状态空间方程（19），即可得到该方程的传递矩阵及其特征值 λ，进而可对系统稳定性进行判定：

3 算例与实验验证

3.1 立式加工中心变工况3维铣削稳定性预报

3.1.1 刀尖点非对称动力学行为辨识

基于某立式加工中心搭建锤击模态实验平台（图3）。所选刀具为两齿硬质合金铣刀，直径φ=20 mm。利用所测刀尖点各向原点与交叉加速度频响函数，可辨识得如表1所示的刀尖点非对称动力学特性参数。

图3 锤击模态实验Fig.3 Impact modal test

表1 刀尖点非对称动力学特性参数Tab. 1 Asymmetrical dynamic parameters of tool center point

3.1.2 铣削力系数标定

利用文献[15]所述实验方法，借助Kistler切削力测试系统对铣削力系数进行标定，如图4所示。所选工件材料为AL7075-T6，辨识得平均铣削力系数为：Kt=804.2 MPa，Kr=177.3 MPa。

图4 铣削力系数标定实验Fig.4 Calibrations of milling force coefficients

3.1.3 过程阻尼参数确定

依据式（7）可知，过程阻尼取决于刀具后刀面磨损带宽度与压入系数、库仑摩擦系数以及切削速度。故一方面，利用超景深显微镜对铣刀后刀面磨损带宽度进行测量（图5），可测得后刀面磨损带宽度W为47 μm；另一方面，由文献[16]可查得压入系数Ksp为1.5×105N/mm3及库仑摩擦系数μ为0.3；最终，将上述W、Ksp、及μ代入式（7）即可求得等效过程阻尼，进一步依据式（22）可得等效黏性阻尼矩阵CP。

图5 刀具后刀面磨损带宽度Fig.5 VB measurements on tool flank face

3.1.4 刀齿切入/切出角

基于该立式加工中心设计如图6所示的变工况铣削工艺过程。其中，刀具进给速度f=200 mm/min，且保持恒定；时变径向切宽可借助其通式（17）推导求得，且ae(0)=5 mm；进而结合式（11）即可得任意时刻刀齿切入/切出角。

图6 变工况铣削工艺过程Fig.6 Milling process on varying cutting condition

3.1.5 变工况3维铣削稳定性预报

综合第3.1.1、3.1.4节及式（19），即可得图3所示立式加工中心铣削动力学状态空间方程。进而，借助时域全离散法，在如表2所示铣削参数创成的变工况设计空间内进行稳定性预报，可绘制如图7所示变工况3维铣削稳定性叶瓣曲面。

表2 变工况铣削参数Tab. 2 Milling parameters on varying cutting conditions

由图7可看出：在3维叶瓣曲面上，部分区域（虚线框）极限轴向切深最高可达20 mm（例如A点），但该点对应极限径向切宽仅为5 mm；在2×103～8×103r/min转速范围内，当径向切宽达到20 mm时，对应极限轴向切深最大值为12.3 mm（B点）；由于充分考虑了过程阻尼，铣削稳定性在主轴低转速区域（双点划线框）得到显著提升。如图7中C点的极限径向切宽ae及极限轴向切深ap均达到了20 mm，但该点对应的主轴转速 Ω仅为1.5×103r/min。

图7 3维铣削稳定性叶瓣曲面Fig.7 3D milling stability lobe surface

由此，主轴转速、轴向切深、径向切宽三者间隐含着复杂的非线性映射关系：一方面，对于某一固定主轴转速 Ω而言，极限稳定径向切宽和极限稳定轴向切深之间呈现较显著的负相关联系，这与实际工程选择切削工艺参数时ap与ae出现相互制约的情况是一致的；另一方面，对于某一固定轴向切深ap而言，可有多种主轴转速 Ω与径向切宽ae与之匹配，从而实现无颤振铣削；此外，过程阻尼对低转速区域的铣削稳定性有提升作用，可为低速重载切削提供一定理论与数据支撑。

3.2 实验验证

为验证上述变工况3维铣削稳定性预报结果，再次利用图4所示实验平台，针对铝合金AL7075-T6分别实施定工况（fixi）和变工况（vari）铣削实验。期间，借助Kistler测力仪对铣削力信号进行采集，并以其时频特征与加工表面质量为依据对极限稳定铣削参数进行判定。

3.2.1 定工况铣削实验（fixi）

利用ae=20.0 mm的参考平面截取图7所示变工况铣削稳定性3维叶瓣曲面，可得如图8（a）所示定工况铣削稳定性叶瓣曲线（实线）。由图8（a）可以看出，铣削过程阻尼及刀尖点非对称动力学行为对极限稳定轴向切深ap，尤其在主轴低转速区，有较显著的影响。为此，在不同转速 Ω（fixi）（i=1，2， ···, 5）下开展0～8 mm变切深铣削实验，以验证上述预报结果。以Ω（fix1）=2 600 r/min为例，该转速下极限轴向切深预报值ap_lim_pre=5.4 mm；图8（b）、（c）分别为该转速下定工况铣削实验所测切削力时频数据与零件加工表面质量观测。

根据图8（b）时域数据可知，随着轴向切深逐渐增大，铣削力幅值亦逐渐增大，并在时域30 s附近出现突变，而该铣削力频域内亦出现了1阶模态附近的颤振频率成分；由图8（c）可以观测到，当ap增大至5.3 mm时工件表面开始出现振纹，此时切深即为转速 Ω（fix1）=2 600 r/min实测极限轴向切深ap_lim_exp。

图8 定工况铣削实验（fix1）Fig.8 Milling experimentation on fixed cutting condition(fix1)

类似地，可完成其他转速对应定工况铣削实验，最终可得如表3所示极限轴向切深测试数据与预报结果对比。其中，最大相对误差为4.0%，表明所提方法能较准确预报定工况铣削稳定性。

表3 定工况极限稳定铣削参数预报值与实测值Tab. 3 Predicted and experimental limit stable milling parameters on fixed cutting conditions

3.2.2 变工况铣削实验（vari）

若以任一主轴转速（Ω ）为参考平面截取图7所示3维叶瓣曲面，可得如图9（a）所示变工况极限稳定铣削参数关系曲线。为验证该预报结果，在对应转速下开展变工况铣削实验。

以 Ω（var2）=2 800 r/min为例，该转速下极限轴向切深ap_lim_pre与径向切宽ae_lim_pre预报值分别为6.4与17.7 mm；图9（b）、（c）则分别为该转速下变工况铣削实验所测铣削力时频特征与零件加工表面质量观测，综合二者分析可知，该转速下极限稳定铣削参数实测值分别为ap_lim_exp=6.8 mm、ae_lim_exp=17.2 mm。表4为不同转速下实验与预报对比结果，其中最大相对误差为5.9%，结果表明该方法在各种工况下均有较好的准确性。