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

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

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

${\rm{MRR}} = N \cdot \varOmega \cdot {a_{\rm{p}}} \cdot {a_{\rm{e}}} \cdot {f_{\rm{t}}}$

(1)

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

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

1 变工况铣削动力学建模 1.1 考虑过程阻尼的动态铣削力

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

$h\left( {{\varphi _j}} \right) = {\left[\!\!\!\! {\begin{array}{*{20}{c}} {\sin \;{\varphi _j}} \\ {\cos \;{\varphi _j}} \end{array}} \!\!\!\!\right]^{\rm{T}}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {{f_{\rm{t}}} + x\left( t \right) - x\left( {t - T} \right)} \\ {y\left( t \right) - y\left( {t - T} \right)} \end{array}}\!\!\!\! \right]$

(2)

式中：f_t为每齿进给量；T为刀齿切削周期；x(t)、x(t–T)和y(t)、y(t–T)则分别为t、(t–T)时刻刀齿j在x和y向的振动位移； ${\varphi _j} $ 为刀齿瞬时径向接触角，且：

${\varphi _j}\left( t \right) = {{2\text{π} \varOmega t} / {60}} + \left( {j - 1} \right){{2\text{π} } / N}$

(3)

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

$\left\{\!\!\!\! \begin{array}{l} {F_{{\rm{t}}j}}\left( t \right) = {F_{{\rm{ts}}}}\left( t \right) + {F_{{\rm{tp}}}}\left( t \right),\\ {F_{{\rm{r}}j}}\left( t \right) = {F_{{\rm{rs}}}}\left( t \right) + {F_{{\rm{rp}}}}\left( t \right) \end{array} \right.$

(4)

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

$\left\{\!\!\!\! \begin{array}{l} {F_{{\rm{t}}j}} = g\left( {{\varphi _j}} \right){K_{\rm{t}}}h\left( {{\varphi _j}} \right){a_{\rm{p}}},\\ {F_{{\rm{r}}j}} = g\left( {{\varphi _j}} \right){K_{\rm{r}}}h\left( {{\varphi _j}} \right){a_{\rm{p}}} \end{array} \right. $

(5)

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

$g\left( {{\varphi _j}\left( t \right)} \right) = \left\{\!\!\!\! {\begin{array}{l} 1,\;\;{{\varphi _j}\left( t \right) \in \left( {{\varphi _{{\rm{st}}}},{\varphi _{{\rm{ex}}}}} \right)}; \\ 0,\;\;{{\varphi _j}\left( t \right) \notin \left( {{\varphi _{{\rm{st}}}},{\varphi _{{\rm{ex}}}}} \right)} \end{array}} \right. $

(6)

式中， ${\varphi _{{\rm{st}}}} $ 和 ${\varphi _{{\rm{ex}}}} $ 分别为切入、切出角。

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

$ {F}_{\rm{rp}}\approx {C}_{\rm{eq}}{\left[\!\!\!\!\begin{array}{c}\dot{x}\left(t\right)\\ \dot{y}\left(t\right)\end{array}\!\!\!\!\right]}^{\rm{T}}\left[\!\!\!\!\begin{array}{c}\mathrm{sin}\;{\varphi }_{j}\left(t\right)\\ \mathrm{cos}\;{\varphi }_{j}\left(t\right)\end{array}\!\!\!\!\right],\;{C}_{\rm{eq}}=\frac{{K}_{\rm{sp}}{a}_{\rm{p}}{W}^{2}}{4v}$

(7)

${F_{{\rm{tp}}}} = \mu {F_{{\rm{rp}}}}$

(8)

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

${\;\;\;\;\;\;\;\begin{aligned}[b] \left[\!\!\!\! {\begin{array}{*{20}{c}} {{F_x}\left( t \right)}\\ {{F_y}\left( t \right)} \end{array}} \!\!\!\!\right] =& {a_{\rm{p}}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {{h_{xx}}}&{{h_{xy}}}\\ {{h_{yx}}}&{{h_{yy}}} \end{array}}\!\!\!\! \right]\left[\!\!\!\! {\begin{array}{*{20}{c}} {{f_{\rm{t}}} + x\left( t \right) - x\left( {t - T} \right)}\\ {y\left( t \right) - y\left( {t - T} \right)} \end{array}} \!\!\!\!\right] + \\ &{C_{{\rm{eq}}}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {{g_{{\rm{p}}xx}}}&{{g_{{\rm{p}}xy}}}\\ {{g_{{\rm{p}}yx}}}&{{g_{{\rm{p}}yy}}} \end{array}} \!\!\!\!\right]\left[\!\!\!\! {\begin{array}{*{20}{c}} {\dot x\left( t \right)}\\ {\dot y\left( t \right)} \end{array}} \!\!\!\!\right]\\[-10pt] \end{aligned}}$

(9)

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

$\left\{\!\!\!\! {\begin{array}{*{20}{l}} {\left\{ \!\!\!\!{\begin{array}{*{20}{c}} {{h_{xx}}}\\ {{g_{{\rm{p}}xx}}}\!\!\!\!\!\! \end{array}} \right\} = \displaystyle\sum\limits_{j = 1}^N {g\left( {{\varphi _j}} \right)} \left[ {\left\{ \!\!\!\!{\begin{array}{*{20}{c}} { - {K_{\rm{t}}}}\\\!\!\!\! { - \mu }\!\!\!\! \end{array}} \!\!\!\!\right\}\cos \;{\varphi _j} - \left\{ \!\!\!\!{\begin{array}{*{20}{c}} {{K_{\rm{r}}}}\\\!\!\!\! 1 \end{array}}\!\!\!\! \right\}\sin \;{\varphi _j}} \right]\sin \;{\varphi _j},}\\ {\left\{ \!\!\!\!{\begin{array}{*{20}{c}} {{h_{xy}}}\\ {{g_{{\rm{p}}xy}}}\!\!\!\!\!\! \end{array}} \right\} = \displaystyle\sum\limits_{j = 1}^N {g\left( {{\varphi _j}} \right)} \left[ {\left\{ \!\!\!\!{\begin{array}{*{20}{c}} { - {K_{\rm{t}}}}\\\!\!\!\! { - \mu } \end{array}}\!\!\!\! \right\}\cos \;{\varphi _j} - \left\{ \!\!\!\!{\begin{array}{*{20}{c}} {{K_{\rm{r}}}}\\\!\!\!\! 1 \end{array}} \!\!\!\!\right\}\sin \;{\varphi _j}} \right]\cos \;{\varphi _j},}\\ {\left\{ \!\!\!\!{\begin{array}{*{20}{c}} {{h_{{\rm{y}}x}}}\\ {{g_{{\rm{py}}x}}}\!\!\!\!\!\! \end{array}} \right\} = \displaystyle\sum\limits_{j = 1}^N {g\left( {{\varphi _j}} \right)} \left[ {\left\{ \!\!\!\!{\begin{array}{*{20}{c}} {{K_{\rm{t}}}}\\\!\!\!\! \mu \end{array}} \!\!\!\!\right\}\sin \;{\varphi _j} - \left\{ \!\!\!\!{\begin{array}{*{20}{c}} {{K_{\rm{r}}}}\\\!\!\!\! 1 \end{array}} \!\!\!\!\right\}\cos \;{\varphi _j}} \right]\sin \;{\varphi _j},}\\ {\left\{\!\!\!\! {\begin{array}{*{20}{c}} {{h_{yy}}}\\ {{g_{{\rm{p}}yy}}}\!\!\!\!\!\! \end{array}} \right\} = \displaystyle\sum\limits_{j = 1}^N {g\left( {{\varphi _j}} \right)} \left[ {\left\{\!\!\!\! {\begin{array}{*{20}{c}} {{K_{\rm{t}}}}\\\!\!\!\! \mu \end{array}} \!\!\!\!\right\}\sin \;{\varphi _j} - \left\{ \!\!\!\!{\begin{array}{*{20}{c}} {{K_{\rm{r}}}}\\\!\!\!\! 1 \end{array}} \!\!\!\!\right\}\cos \;{\varphi _j}} \right]\cos \;{\varphi _j}} \end{array}} \right.$

(10)

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

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

某凸曲面顺铣工艺过程如图2所示，其中：a_e(t)为时变径向切削宽度；D₁与D₂分别位于待加工表面与已加工表面，且分别对应t时刻刀齿切入与切出点；O_b(t)与O_w(t)则分别为当前时刻上述两点处对应的瞬时曲率中心；X_wY_wZ_w与X_sY_sZ_s分别为工件坐标系与刀具坐标系，且二坐标系在XY平面内的瞬时夹角为θ(t)； ${\varphi _{{\rm{st}}}}$ (t)与 ${\varphi _{{\rm{ex}}}}$ (t)则分别刀齿切入角与切出角，且^[12]：

$\left\{\!\!\!\! \begin{array}{l} {\varphi _{{\rm{st}}}} = \arccos \left[ {{{{a_{\rm{e}}}\left( t \right)} / R} - 1} \right], \\ {\varphi _{{\rm{ex}}}} = \text{π} \end{array} \right.$

(11)

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

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

${f_{\rm{b}}}\left( {{x_{\rm{b}}}\left( t \right),{y_{\rm{b}}}\left( t \right)} \right){\rm{ = }}0$

(12)

${f_{\rm{w}}}\left( {{x_{\rm{w}}}\left( t \right),{y_{\rm{w}}}\left( t \right)} \right){\rm{ = }}0$

(13)

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

$\frac{{x - {x_{D2}}}}{{{x^\prime_{\rm{w}}} \left( t \right)}} = \frac{{y - {y_{D2}}}}{{{y^\prime_{\rm{w}}} \left( t \right)}}$

(14)

$\frac{{\left( {x - {x_{D2}}} \right)}}{{{y^\prime_{\rm{w}}} \left( t \right)}} = \frac{{\left( {y - {y_{D2}}} \right)}}{{{x^\prime_{\rm{w}}} \left( t \right)}}$

(15)

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

$\frac{{\left( {x - {x_{D1}}} \right)}}{{{x^\prime_{\rm{w}}}\left( t \right)}} = \frac{{\left( {y - {y_{D1}}} \right)}}{{{y^\prime_{\rm{w}}} \left( t \right)}}$

(16)

联立式（15）、（16）即可得二者交点D₃坐标（x_D3，y_D3），则：

${a_{\rm{e}}}\left( t \right) = \left\| {D{}_2{D_3}} \right\|{\rm{ = }}\sqrt {{{\left( {{x_{D3}} - {x_{D2}}} \right)}^2} + {{\left( {{y_{D3}} - {y_{D2}}} \right)}^2}} $

(17)

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

1.3 变工况铣削动力学模型

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

${{M}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {\ddot x\left( t \right)} \\ {\ddot y\left( t \right)} \end{array}} \!\!\!\!\right] + {{C}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {\dot x\left( t \right)} \\ {\dot y\left( t \right)} \end{array}} \!\!\!\!\right] + {{K}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {x\left( t \right)} \\ {y\left( t \right)} \end{array}} \!\!\!\!\right]{\rm{ = }}\left[\!\!\!\! {\begin{array}{*{20}{c}} {{F_x}\left( t \right)} \\ {{F_y}\left( t \right)} \end{array}} \!\!\!\!\right]$

(18)

式中：M、C、K分别为切削系统质量矩阵、阻尼矩阵与刚度矩阵，且 ${{M}} \!=\! \left[\!\!\!\! {\begin{array}{*{20}{c}} {{m_{xx}}}\;\;{{m_{yx}}} \\ {{m_{xy}}}\;\;{{m_{yy}}} \end{array}}\!\!\!\! \right]$ ， ${{C}} \!=\! \left[\!\!\!\! {\begin{array}{*{20}{c}} {{c_{xx}}}\;\;{{c_{yx}}} \\ {{c_{xy}}}\;\;{{c_{yy}}} \end{array}}\!\!\!\! \right]$ ， ${{K}}\! =\! \left[ \!\!\!\!{\begin{array}{*{20}{c}} {{k_{xx}}}\;\;{{k_{yx}}} \\ {{k_{xy}}}\;\;{{k_{yy}}} \end{array}}\!\!\!\! \right]$ ；F_x(t)、F_y(t)分别为x、y向铣削力。

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

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

${{\dot X}}\left( t \right) = {{{A}}_0}{{X}}\left( t \right) + {{A}}\left( t \right)\left[ {{{X}}\left( t \right) - {{X}}\left( {t - T} \right)} \right] + {{{F}}_{{f}}}\left( t \right)$

(19)

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

${{X}}\left( t \right) = {\left[ {{{q}}\left( t \right){\rm{ }}{{p}}\left( t \right)} \right]^{\rm{T}}}$

(20)

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

$\left\{\!\!\!\! {\begin{array}{*{20}{l}} {{{q}}\left( t \right) = {{\left[ {x\left( t \right){\rm{ }}y\left( t \right)} \right]}^{\rm{T}}}}, \\ {{{p}}\left( t \right) = {{M\dot q}}\left( t \right) + \left( {{{C}} - {{{C}}_{\rm{P}}}} \right){{q}}\left( t \right)/2} \end{array}} \right.$

(21)

${{{C}}_{\rm{P}}} = {C_{{\rm{eq}}}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {{g_{{\rm{p}}xx}}}&{{g_{{\rm{p}}xy}}} \\ {{g_{{\rm{p}}yx}}}&{{g_{{\rm{p}}yy}}} \end{array}} \!\!\!\!\right]$

(22)

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

${\rm{max}}\left( {\left| \lambda \right|} \right)\left\{\!\!\!\! \begin{array}{l} < 1,\;\;{\text{稳定}} ;\\ = 1,\;\;{\text{临界稳定}} ;\\ > 1,\;\;{\text{不稳定}} \end{array} \right.$

(23)

3 算例与实验验证 3.1 立式加工中心变工况3维铣削稳定性预报 3.1.1 刀尖点非对称动力学行为辨识

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

3.1.2 铣削力系数标定

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

3.1.3 过程阻尼参数确定

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

3.1.4 刀齿切入/切出角

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

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

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

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

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

3.2 实验验证

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

3.2.1 定工况铣削实验（fix_i）

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

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

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

3.2.2 变工况铣削实验（var_i）

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

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

4 结　论

1）在综合考虑铣削过程阻尼、刀尖点非对称动力学行为及时变切宽的基础上，借助时域全离散法针对变工况提出了3维铣削稳定性预报方法体系，为实现面向变工况开展铣削工艺参数优化提供了理论依据与数据支撑。

2）利用所提方法，针对某立式加工中心进行了变工况3维铣削稳定性预报。预报结果表明：过程阻尼及刀尖点非对称动力学行为对低转速区铣削稳定性有较显著的提升作用；而在高转速区，极限稳定轴向切深与径向切宽之间存在较明显的负相关联系。

3）通过开展定工况及变工况铣削实验，对所提方法进行了验证。数据对比分析显示，极限稳定铣削参数（轴向切深、径向切宽）的预报值与实验值吻合程度较高，最大相对误差为5.9%，由此说明所提预报方法具有较高的准确性。