基于加工振动的立铣刀关键几何参数设计方法

摘要

本发明提供的基于加工振动的立铣刀关键几何参数设计方法，包括：对加工系统进行模态试验，得到关键模态参数；对加工系统进行动力学建模，建立多时滞二阶微分动力学方程；建立并得到变换后的状态空间方程；利用GRK法判定加工系统的稳定性并获得加工参数空间中的稳定性图谱Lobe图；改变设计参数即铣刀齿间距和螺旋角的值以获得不同的Lobe图；以获得最大加工效率为目标，通过比较不同设计参数条件下的Lobe图，得到优化后的立铣刀齿间距和螺旋角。本发明与传统等齿距标准铣刀加工相比，采用GRK法得到加工系统的动力学特性，获得优化后的立铣刀关键几何设计参数即齿间距和螺旋角，极大地提高了加工效率。

说明书

技术领域

本发明涉及机械加工技术领域，具体地，涉及一种基于加工振动的立铣刀关键几何参数设计方法。 

背景技术

在机械加工领域，铣削是最为常见和常用的加工方式之一。随着科技的进步，人们对于铣削加工质量的要求也越来越高，而加工振动是影响加工质量的关键因素之一。减少振动，尤其是消除剧烈的振动即颤振的影响，是保证加工质量的必然要求。颤振的成因包括再生效应、模态耦合、磨擦作用、力热耦合等，其中再生效应是引起铣削颤振的最主要原因。 

再生效应的避免有许多方法，其中采用不等齿距铣刀破坏再生颤振的产生机制是行之有效的手段之一。现有的某些商业铣刀也采用了变齿距设计来提高加工稳定性，但是商业铣刀在设计前并不明确其具体的加工环境与加工对象，因此很难将不等齿距铣刀的避振性能发挥到最大。因此，站在用户的角度，根据具体的加工环境定制具有最佳避振性能的铣刀具有十分重要的意义和前景。 

发明内容

本发明针对现有技术中存在的上述不足，提供了一种基于加工振动的立铣刀关键几何参数设计方法，该方法利用GRK法(广义Runge-Kutta法，Generalized Runge-Kutta Method)，依据具体加工环境与条件，考虑加工振动的影响，设计立铣刀关键几何参数(齿距及螺旋角)，在保证无再生颤振高质量加工的前提下，使加工参数(主轴转速、轴向切削深度)空间内的稳定切削区域最大化，进而获得最大的加工效率，产生良好的经济效益。 

本发明是通过以下技术方案实现的。 

一种基于加工振动的立铣刀关键几何参数设计方法，包括如下步骤： 

步骤1：以材质、长度、直径相同的标准等距立铣刀作为替代刀具进行模态试验，获得刀具-机床加工系统的模态参数； 

步骤2：对加工系统进行动力学建模，得到刀具-机床加工系统的动力学方程； 

步骤3：对动力学方程进行状态空间变换，得到状态空间方程； 

步骤4：应用GRK法对刀具-机床加工系统的状态空间方程进行稳定性分析，获得刀具-机床加工系统在加工参数空间的稳定性图谱，即Lobe图； 

步骤5：改变替代刀具不同齿距参数和螺旋角参数，获得不同齿距参数和螺旋角参数条件下的刀具-机床加工系统Lobe图； 

步骤6：以获得最大加工效率为目标，通过比较不同刀具参数条件下的Lobe图，得到优化后的立铣刀关键几何参数，即齿距和螺旋角。 

优选地，所述步骤1，具体为： 

步骤1.1，以具体加工环境为依托，在机床上装夹与所设计立铣刀材质、长度、直径参数相同的标准等距立铣刀，进行刀具-机床加工系统的模态试验； 

步骤1.2，采用脉冲锤敲击实验，测得刀具-机床加工系统所需的模态参数，即：模态质量矩阵M，模态阻尼矩阵C，模态刚度矩阵K。 

优选地，所述步骤2，具体为： 

步骤2.1，在二自由度系统下，不考虑模态耦合作用，所述刀具-机床加工系统的动力学方程，如公式(1)所示： 

$M\left(\begin{array}{c}\stackrel{..}{x}\\ \stackrel{..}{y}\end{array}\right)+C\left(\begin{array}{c}\stackrel{.}{x}\\ \stackrel{.}{y}\end{array}\right)+K\left(\begin{array}{c}x\\ y\end{array}\right)=F---\left(1\right)$

其中，质量矩阵$M=\left(\begin{array}{cc}{m}_x& 0\\ 0& m_y\end{array}\right),$阻尼矩阵$C=\left(\begin{array}{cc}{c}_x& 0\\ 0& c_y\end{array}\right),$刚度矩阵$K=\left(\begin{array}{cc}{k}_x& 0\\ 0& k_y\end{array}\right),$切削力矩阵m、c、k、F分别表示质量、阻尼、刚度、切削力；x，y分别表示两个方向； 

步骤2.2，沿轴向将替代刀具离散为Z个厚度很小的微圆盘，对于每个微元，替代刀具在第j个刀齿处所受的切向力和径向力如公式(2)所示： 

其中，dF_jt(t，z)和dF_jn(t，z)分别表示刀齿j在高度为z的圆盘处于时刻t所受到的切向微元力和径向微元力，表示圆弧角，表示开关函数，h_j(t，z)表示切厚，K_tc、 K_te、K_nc和K_ne分别表示切向切削力系数、切向刃口力系数、径向切削力系数和径向刃口力系数； 

圆弧角的表达式如公式(3)所示： 

其中，Ω，β，和R分别表示主轴转速(rpm)，螺旋角(rad)和半径(m)，表示刀齿j与紧邻的上一个刀齿(j+1)之间的圆心角；对于等齿距铣刀，其中N为刀齿数目； 

开关函数用于判断对应的微元是否正在切削，其表达式如公式(4)所示： 

其中，分别表示刀齿切入、切出角，表示除以2π的余数； 

所述切屑瞬时切厚表达式为： 

其中，f_j表示刀齿j在x方向的进给率； 

步骤2.3，定义a/D为刀具的径向切深比，其中a为径向切深，D为刀具直径；对于替代刀具的逆铣，对于替代刀具的顺铣 

步骤2.4，为了获得在刀具正交坐标系下的动力学方程，将切向力和径向力向x、y方向投影，结果如公式(6)所示： 

步骤2.5，沿轴向对公式(6)积分，得到替代刀具刀齿j在正交坐标系下的切削力，如公式(7)所示： 

步骤2.6，整个替代刀具所受切削力为各个刀齿所受切削力在正交方向的叠加，如公式(8)所示： 

$\left(\begin{array}{c}{F}_x=\underset{j=1}{\overset{N}{\Sigma }}{F}_{\mathrm{jx}}\\ F_y=\underset{j=1}{\overset{N}{\Sigma }}{F}_{\mathrm{jy}}\end{array}\right)---\left(8\right)$

所以，两自由度端铣系统(加工系统)的动力学方程如下： 

$M\stackrel{..}{q}\left(t\right)+C\stackrel{.}{q}\left(t\right)+\mathrm{Kq}\left(t\right)=\underset{j=1}{\overset{N}{\Sigma }}{K}_j\left(t\right)[q\left(t\right)-q(t-T_j)]+F_0---\left(9\right)$

其中，q(t)＝[x(t)，y(t)]^T， 

其中，M为模态质量矩阵，K为模态刚度矩阵，K_j(t)为切削系数矩阵，T_j为齿间距， 为加速度状态向量，为速度状态向量，q(t)为位移状态向量，F₀为与动态切厚无关的切削力分量，t为时间，x(t)为x方向振动位移，y(t)为y方向振动位移，a_p为轴向切深，表示圆弧角，表示开关函数，f_j表示刀齿j在x方向的进给率，K_tc、K_te、K_nc和K_ne分别表示切向切削力系数、切向刃口力系数、径向切削力系数和径向刃口力系数，dz为轴向微元厚度，N为刀具沿轴向离散微元数目。 

F₀与动态切厚无关，F₀对线性动力学稳定性图谱的影响忽略不计。 

优选地，所述步骤3，具体为： 

令则经过状态空间变换之后，系统动力学方程变为如下状态空间表达式： 

$\stackrel{.}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\underset{j=1}{\overset{N}{\Sigma }}{B}_j\left(t\right)[x\left(t\right)-x(t-T_j)]---\left(10\right)$

其中，$A=\left(\begin{array}{cc}0& I\\ {-M}^{-1}K& -M^{-1}C\end{array}\right),B_j\left(t\right)=\left(\begin{array}{cc}0& 0\\ M^{-1}{K}_j\left(t\right)& 0\end{array}\right);$

其中，I为单位矩阵，为x方向振动速度，为y方向振动速度。 

优选地，所述步骤4，具体为： 

公式(10)所示一阶微分方程的解析解为： 

$x\left(t\right)=e^{A(t-t_0)}x\left(t_0\right)+\underset{j=1}{\overset{N}{\Sigma }}{\int }_{t_0}^t\left\{e^{A(t-\xi )}{B}_j\left(\xi \right)\right[x\left(\xi \right)-x(\xi -T_j)\left]\right\}\mathrm{d\xi }---\left(11\right)$

公式(11)为基本矩阵运算表达式，其中，t₀表示计算起始时间点；为了推导简洁，记 K_j(t，ξ，x(ξ))＝e^A(t-ξ)B_j(ξ)[x(ξ)-x(ξ-T_j)]，x(t_i)简写为x_i；应用GRK法求解多时滞微分方程(10)的具体步骤如下： 

步骤4.1，假设计算起始点为(t_2i，x_2i)，i＝0，1，…，i_m，其中i_m是一个与时间离散数m有关的取整函数： 

$i_m=\mathrm{ceil}\left(\frac{m}{2}\right)-1---\left(12\right)$

其中，ceil的作用为取比自变量大的右端紧邻整数值(如果自变量本身为整数点则取自变量本身)，举例说明ceil(2.1)＝3； 

步骤4.2，假设x_2i+1已知，则利用Simpson公式可得： 

$x_{2i+2}=g(t_{2i+2},t_{2i})+\frac{h}{3}\underset{j=1}{\overset{N}{\Sigma }}\{K_j(t_{2i+2},t_{2i},x_{2i})+4K_j(t_{2i+2},t_{2i+1},x_{2i+1})+K_j(t_{2i+2},t_{2i+2},x_{2i+2})\}---\left(13\right)$

其中，h表示离散步长，h＝T/m，T为时间周期； 

步骤4.3，计算x_2i+1，引入中间点x_2i+1/2，然后利用积分形式的四阶Runge-Kutta法可得： 

$\left(\begin{array}{c}{x}_{2i+1}=g(t_{2i+1},t_{2i})+\frac{h}{6}\underset{j=1}{\overset{N}{\Sigma }}\{K_j(t_{2i+1},t_{2i},x_{2i})+2K_j(t_{2i+1},t_{2i+1/2},x_{2i+1/2})\\ +2K_j(t_{2i+1},t_{2i+1/2}{,x}_{2i+1/2})+K_j(t_{2i+1},t_{2i+1},x_{2i+1})\}\end{array}\right)---\left(14\right)$

中间点x_2i+1/2用三点Lagrange插值函数计算， 

$x_{2i+1/2}\approx \frac{3}{8}{x}_{2i}+\frac{3}{4}{x}_{2i+1}-\frac{1}{8}{x}_{2i+2}---\left(15\right)$

因此，公式(14)中的中间点可以表示为： 

$K_j(t_{2i+1},t_{2i+1/2},x_{2i+1/2})=K_j(t_{2i+1},t_{2i+1/2},\frac{3}{8}{x}_{2i}+\frac{3}{4}{x}_{2i+1}-\frac{1}{8}{x}_{2i+2})---\left(16\right)$

x(t_2i-T_j)简写为x_2i-Tj。根据公式(13)-(16)，可以推得相应的迭代公式。偶数项公式可由(13)求得： 

$F_{2i}{x}_{2i}+F_{2i+1}{x}_{2i+1}+F_{2i+2}{x}_{2i+2}=\underset{j=1}{\overset{N}{\Sigma }}(F_{2i-T_j}{x}_{2i-T_j}+F_{2i+1-T_j}{x}_{2i+1-T_j}+F_{2i+2-T_j}{x}_{2i+2-T_j})---\left(17\right)$

其中，$F_{2i}=-e^{\mathrm{Ah}}-\frac{h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}}{B}_{j,2i}\right)-\frac{3}{8}·\frac{4h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}/2}{B}_{j,2i+1/2}\right),$

$F_{2i+1}=I-\frac{3}{4}·\frac{4h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}/2}{B}_{j,2i+1/2}\right)-\frac{h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(B_{j,2i+1}\right),F_{2i+2}=\frac{1}{8}·\frac{4h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}/2}{B}_{j,2i+1/2}\right),$

$F_{2i-T_j}=-\frac{h}{6}{e}^{\mathrm{Ah}}{B}_{j,2i}-\frac{3}{8}·\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{j,2i+1/2},F_{2i+1-T_j}=-\frac{3}{4}·\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{j,2i+1/2}-\frac{h}{6}{B}_{j,2i+1},$

$F_{2i+2-T_j}=\frac{1}{8}·\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{j,2i+1/2};$

步骤4.4，奇数项公式可由(14)求得： 

$G_{2i}{x}_{2i}+G_{2i+1}{x}_{2i+1}+G_{2i+2}{x}_{2i+2}=\underset{j=1}{\overset{N}{\Sigma }}(G_{2i-T_j}{x}_{2i-T_j}+G_{2i+1-T_j}{x}_{2i+1-T_j}+G_{2i+2-T_j}{x}_{2i+2-T_j})---\left(18\right)$

其中，$G_{2i}={-e}^{2\mathrm{Ah}}-\frac{h}{3}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{2\mathrm{Ah}}{B}_{j,2i}\right),G_{2i+1}=-\frac{4h}{3}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}}{B}_{j,2i+1}\right),G_{2i+2}=I-\frac{h}{3}\underset{j=1}{\overset{N}{\Sigma }}{B}_{j,2i+2},$

$G_{2i-T_j}=-\frac{h}{3}{e}^{2\mathrm{Ah}}{B}_{j,2i},G_{2i+1-T_j}=-\frac{4h}{3}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}}{B}_{j,2i+1}\right),G_{2i+2-T_j}=-\frac{h}{3}{B}_{j,2i+2}.$

基于公式(17)-(18)，动力学方程两个周期之间的离散映射关系为： 

$\left(\begin{array}{c}{P}_1{[{x_0}^T,{x_1}^T,{x_2}^T,{x_3}^T,{x_4}^T,···,{x_{{2i}_m+1}}^T,{x_{{2i}_m+2}}^T]}^{T=}\\ Q{[{x_{0-T}}^T,{x_{1-T}}^T,{x_{2-T}}^T,{x_{3-T}}^T,{x_{4-T}}^T,···,{x_{{2i}_m+1-T}}^T,{x_{{2i}_m+2-T}}^T]}^T\\ +P_2{[{x_0}^T,{x_1}^T,{x_2}^T,{x_3}^T,{x_4}^T,···,{x_{{2i}_m+1}}^T,{x_{{2i}_m+2}}^T]}^T\end{array}\right)---\left(19\right)$

其中， 

P₁、P₂和Q中的分块矩阵的位置由公式(17)-(19)决定；具体的，P₁中的分块矩阵的位置是固定的，如公式(20)所示，P₂和Q中的分块矩阵的行位置也是固定的，要与P₁对应，而P₂和Q中的分块矩阵的列位置不是固定的，与对应时滞T_i的大小有关；式(20)-(22)中，F_x表示奇数项对应的分块矩阵，下标x表示对应的离散时刻；G_x表示偶数项对应的分块矩阵，下标x表示对应的离散时刻； 

状态转移矩阵可表示为： 

其中，表示矩阵(P₁-P₂)的Moor-Penrose广义逆； 

根据Floquet理论，如果状态转移矩阵Φ的所有特征值的模都小于1，则系统是稳定的，反之，如果状态转移矩阵Φ的任何一个特征值的模大于1，则系统是不稳定的；因此，根据Floquet理论画出系统在切削参数空间的稳定边界，即稳定性图谱Lobe图。 

优选地，所述步骤5，具体为： 

改变替代刀具的螺旋角β和齿间距T_j，重复步骤4，获得新的参数条件下的稳定性图谱；设计齿间距时，将铣刀动平衡的约束条件考虑在内，且满足T为时间周期。 

优选地，所述步骤6，具体为： 

根据不同齿间距和螺旋角条件下得到的稳定性图谱，根据具体加工时所采用的转速范围，以获得最大加工效率为优化目标，获得优化后的立铣刀关键几何参数即齿间距和螺旋角。 

与现有技术相比，本发明具有如下的有益效果： 

本发明提出的基于加工振动的立铣刀关键几何参数设计方法，是用户根据具体加工背景定制设计立铣刀关键几何参数即齿距和螺旋角的新方法-GRK法，与标准等距铣刀加工相比，加工效率得到了极大的提高，与现有商业不等齿距铣刀相比，也使铣刀的稳定切削范围得到很大提升。 

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显： 

图1为不等齿距立铣刀端铣加工系统示意图，以顺铣为例；图中，x和y分别表示正交坐标系的两个方向，k和c分别表示刚度系数和阻尼系数，F_n和F_t分别表示径向和切向切削力，表示刀齿的圆心角，f_z表示工件的进给方向。 

图2为考虑螺旋角的立铣刀沿轴向离散示意图；图中，表示齿间距，β表示螺旋角，dz表示微圆盘的厚度。 

图3为不等齿距铣刀加工实例示意图；图中，序号1、2、3、4分别表示四个刀齿的标号，T₁、T₂、T₃和T₄分别表示刀齿1与2之间的延迟时间、刀齿2与3之间的延迟时间、刀齿3与4之间的延迟时间、刀齿4与1之间的延迟时间。 

图4为不等齿距铣刀加工实例对应的Lobe图；图中，实线表示用GRK法绘出的 Lobe图，带星号的实线表示用文献1中频域法绘出的Lobe图，两种方法均没有考虑螺旋角的影响。 

图5为不同齿距条件下对应的Lobe图；图中，螺旋角为0°，粗虚线代表的齿间角为90°-90°-90°-90°，细虚线代表的齿间角为80°-100°-80°-100°，粗实线代表的齿间角为70°-110°-70°-110°，细实线代表的齿间角为60°-120°-60°-120°。 

图6为不同径向切深比条件下对应的Lobe图；齿间角为70°-110°-70°-110°，螺旋角为0°，粗虚线表示的径向切深比为a/D＝0.5，细实线表示的径向切深比为a/D＝0.3，粗实线表示的径向切深比为a/D＝0.1。 

图7为高径向切深比、不同螺旋角条件下对应的Lobe图；图中，齿间角为70°-110°-70°-110°，径向切深比为a/D＝0.5，细实线代表的螺旋角为0°，粗实线代表的螺旋角为30°，细虚线代表的螺旋角为45°。 

图8为低径向切深比、不同螺旋角条件下对应的Lobe图；图中，齿间角为70°-110°-70°-110°，径向切深比为a/D＝0.05，细虚线代表的螺旋角为0°，细实线代表的螺旋角为60°，粗虚线代表的螺旋角为30°，粗实线代表的螺旋角为45°。 

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。 

请同时参阅图1至图8。 

以长度、直径、材质等参数相同的标准等距立铣刀为模态试验对象，测得刀具-机床系统模态参数； 

对加工系统进行动力学建模，建立多时滞二阶微分方程。 

对动力学方程进行状态空间变换，建立状态空间方程； 

用GRK法求取加工系统在加工参数空间(主轴转速-轴向切深空间)的无颤振加工区域，即稳定性图谱Lobe图； 

改变立铣刀的齿距及螺旋角，获得不同设计参数下的Lobe图； 

以获得最大加工效率为优化目标，通过比较不同设计参数条件下的Lobe图，得到优化后的立铣刀关键几何参数，即齿距与螺旋角。 

具体地，本实施例提供了一种基于加工振动的立铣刀关键几何参数设计方法，包括如下步骤： 

步骤1：以材质、长度、直径相同或相似的标准等距立铣刀为替代刀具，获得刀具-机床加工系统的模态参数； 

步骤2：建立加工系统的动力学方程； 

步骤3：对动力学方程进行状态空间变换； 

步骤4：应用GRK法对加工系统进行稳定性分析，获得加工参数空间的稳定性图谱，即Lobe图； 

步骤5：获得不同齿距参数和螺旋角参数条件下的Lobe图； 

步骤6：以获得最大加工效率为目标，通过比较不同设计参数条件下的Lobe图，得到优化后的立铣刀关键几何参数即齿距和螺旋角。 

优选地，所述步骤1，具体为： 

以具体加工环境为依托，在机床上装夹与所设计立铣刀材质、长度、直径等参数相同或相近的标准等距铣刀，进行刀具-机床系统的模态试验，可以采用脉冲锤敲击实验，测得设计所需的模态参数，即：模态质量矩阵M，模态阻尼矩阵C，模态刚度矩阵K。 

优选地，所述步骤2，具体为： 

按照图1所示加工系统示意图，以二自由度系统为例，不考虑模态耦合作用，所述系统动力学方程，如公式(1)所示： 

$M\left(\begin{array}{c}\stackrel{..}{x}\\ \stackrel{..}{y}\end{array}\right)+C\left(\begin{array}{c}\stackrel{.}{x}\\ \stackrel{.}{y}\end{array}\right)+K\left(\begin{array}{c}x\\ y\end{array}\right)=F---\left(1\right)$

其中，质量矩阵$M=\left(\begin{array}{cc}{m}_x& 0\\ 0& m_y\end{array}\right),$阻尼矩阵$C=\left(\begin{array}{cc}{c}_x& 0\\ 0& c_y\end{array}\right),$刚度矩阵$K=\left(\begin{array}{cc}{k}_x& 0\\ 0& k_y\end{array}\right),$切削力矩阵m，c，k，F分别表示质量、阻尼、刚度、切削力；x，y分别表示两个方向。 

如图2所示，沿轴向将刀具离散为Z个厚度很小的微圆盘。对于每个微元，铣刀在第j个刀齿处所受的切向力和径向力如公式(2)所示： 

其中，dF_jt(t，z)和dF_jn(t，z)分别表示刀齿j在高度为z的圆盘处于时刻t所受到的切向微 元力和径向微元力，表示圆弧角，表示开关函数，h_j(t，z)表示切厚，K_tc，K_te，K_nc和K_ne分别表示切向切削力系数，切向刃口力系数，径向切削力系数和径向刃口力系数。 

圆弧角的表达式如公式(3)所示： 

其中，Ω，β，和R分别表示主轴转速(rpm)，螺旋角(rad)，和半径(m)，表示刀齿j与紧邻的上一个刀齿(j+1)之间的圆心角。对于等齿距铣刀，其中N为刀齿数目。 

开关函数用于判断对应的微元是否正在切削，其表达式如公式(4)所示： 

其中，分别表示刀齿切入、切出角，表示除以2π的余数。 

定义a/D为刀具的径向切深比，其中a为径向切深，D为刀具直径。对于逆铣， 对于顺铣

切屑瞬时切厚表达式为： 

其中，f_j表示刀齿j在x方向的进给率。 

为了获得在刀具正交坐标系下的动力学方程，将切向力和径向力向x、y方向投影，结果如公式(6)所示： 

沿轴向对公式(6)积分，得到刀齿j在正交坐标系下的切削力，如公式(7)所示： 

整个刀具所受切削力为各个刀齿所受切削力在正交方向的叠加，如公式(8)所示： 

$\left(\begin{array}{c}{F}_x=\underset{j=1}{\overset{N}{\Sigma }}{F}_{\mathrm{jx}}\\ F_y=\underset{j=1}{\overset{N}{\Sigma }}{F}_{\mathrm{jy}}\end{array}\right)---\left(8\right)$

所以，两自由度端铣系统动力学方程如下： 

$M\stackrel{..}{q}\left(t\right)+C\stackrel{.}{q}\left(t\right)+\mathrm{Kq}\left(t\right)=\underset{j=1}{\overset{N}{\Sigma }}{K}_j\left(t\right)[q\left(t\right)-q(t-T_j)]+F_0---\left(9\right)$

其中，q(t)＝[x(t)，y(t)]^T， 

由于F₀与动态切厚无关，其对线性动力学稳定性图谱的影响可以忽略不计，所以在下述分析中，F₀没有被包含在内。 

优选地，所述步骤3，具体为： 

令则经过状态空间变换之后，系统动力学方程变为如下状态空间表达式： 

$\stackrel{.}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\underset{j=1}{\overset{N}{\Sigma }}{B}_j\left(t\right)[x\left(t\right)-x(t-T_j)]---\left(10\right)$

其中，$A=\left(\begin{array}{cc}0& I\\ {-M}^{-1}K& -M^{-1}C\end{array}\right),B_j\left(t\right)=\left(\begin{array}{cc}0& 0\\ M^{-1}{K}_j\left(t\right)& 0\end{array}\right).$

优选地，所述步骤4，具体为： 

公式(10)所示一阶微分方程的解析解为： 

$x\left(t\right)=e^{A(t-t_0)}x\left(t_0\right)+\underset{j=1}{\overset{N}{\Sigma }}{\int }_{t_0}^t\left\{e^{A(t-\xi )}{B}_j\left(\xi \right)\right[x\left(\xi \right)-x(\xi -T_j)\left]\right\}\mathrm{d\xi }---\left(11\right)$

其中，t₀表示计算起始时间点。为了推导简洁，记K_j(t，ξ，x(ξ))＝e^A(t-ξDB_j(ξ)[x(ξ)-x(ξ-T_j)]，x(t_i)简写为x_i。那么，应用GRK法求解多时滞微分方程(10)的具体步骤如下。 

假设计算起始点为(t_2i，x_2i)，i＝0，1，…，i_m，其中i_m是一个与时间离散数m有关的取整函数。 

$i_m=\mathrm{ceil}\left(\frac{m}{2}\right)-1---\left(12\right)$

其中，ceil的作用为取比自变量大的右端紧邻整数值(如果自变量本身为整数点则取自变量本身)，举例说明ceil(2.1)＝3。 

假设x_2i+1已知，则利用Simpson公式可得： 

$x_{2i+2}=g(t_{2i+2},t_{2i})+\frac{h}{3}\underset{j=1}{\overset{N}{\Sigma }}\{K_j(t_{2i+2},t_{2i},x_{2i})+4K_j(t_{2i+2},t_{2i+1},x_{2i+1})+K_j(t_{2i+2},t_{2i+2},x_{2i+2})\}---\left(13\right)$

其中，h表示离散步长，h＝T/m。 

为了计算x_2i+1，引入中间点x_2i+1/2，然后利用积分形式的四阶Runge-Kutta法可得： 

$\left(\begin{array}{c}{x}_{2i+1}=g(t_{2i+1},t_{2i})+\frac{h}{6}\underset{j=1}{\overset{N}{\Sigma }}\{K_j(t_{2i+1},t_{2i},x_{2i})+2K_j(t_{2i+1},t_{2i+1/2},x_{2i+1/2})\\ +2K_j(t_{2i+1},t_{2i+1/2}{,x}_{2i+1/2})+K_j(t_{2i+1},t_{2i+1},x_{2i+1})\}\end{array}\right)---\left(14\right)$

中间点x_2i+1/2用三点Lagrange插值函数计算， 

$x_{2i+1/2}\approx \frac{3}{8}{x}_{2i}+\frac{3}{4}{x}_{2i+1}-\frac{1}{8}{x}_{2i+2}---\left(15\right)$

因此，公式(14)中的中间点可以表示为： 

$K_j(t_{2i+1},t_{2i+1/2},x_{2i+1/2})=K_j(t_{2i+1},t_{2i+1/2},\frac{3}{8}{x}_{2i}+\frac{3}{4}{x}_{2i+1}-\frac{1}{8}{x}_{2i+2})---\left(16\right)$

x(t_2i-T_j)简写为根据公式(13)-(16)，可以推得相应的迭代公式。偶数项公式可由(13)求得： 

$F_{2i}{x}_{2i}+F_{2i+1}{x}_{2i+1}+F_{2i+2}{x}_{2i+2}=\underset{j=1}{\overset{N}{\Sigma }}(F_{2i-T_j}{x}_{2i-T_j}+F_{2i+1-T_j}{x}_{2i+1-T_j}+F_{2i+2-T_j}{x}_{2i+2-T_j})---\left(17\right)$

其中，$F_{2i}=-e^{\mathrm{Ah}}-\frac{h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}}{B}_{j,2i}\right)-\frac{3}{8}·\frac{4h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}/2}{B}_{j,2i+1/2}\right),$

$F_{2i+1}=I-\frac{3}{4}·\frac{4h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}/2}{B}_{j,2i+1/2}\right)-\frac{h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(B_{j,2i+1}\right),F_{2i+2}=\frac{1}{8}·\frac{4h}{6}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}/2}{B}_{j,2i+1/2}\right),$

$F_{2i-T_j}=-\frac{h}{6}{e}^{\mathrm{Ah}}{B}_{j,2i}-\frac{3}{8}·\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{j,2i+1/2},F_{2i+1-T_j}=-\frac{3}{4}·\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{j,2i+1/2}-\frac{h}{6}{B}_{j,2i+1},$

$F_{2i+2-T_j}=\frac{1}{8}·\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{j,2i+1/2}.$

奇数项公式可由(14)求得： 

$G_{2i}{x}_{2i}+G_{2i+1}{x}_{2i+1}+G_{2i+2}{x}_{2i+2}=\underset{j=1}{\overset{N}{\Sigma }}(G_{2i-T_j}{x}_{2i-T_j}+G_{2i+1-T_j}{x}_{2i+1-T_j}+G_{2i+2-T_j}{x}_{2i+2-T_j})---\left(18\right)$

其中，$G_{2i}={-e}^{2\mathrm{Ah}}-\frac{h}{3}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{2\mathrm{Ah}}{B}_{j,2i}\right),G_{2i+1}=-\frac{4h}{3}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}}{B}_{j,2i+1}\right),G_{2i+2}=I-\frac{h}{3}\underset{j=1}{\overset{N}{\Sigma }}{B}_{j,2i+2},$$G_{2i-T_j}=-\frac{h}{3}{e}^{2\mathrm{Ah}}{B}_{j,2i},G_{2i+1-T_j}=-\frac{4h}{3}\underset{j=1}{\overset{N}{\Sigma }}\left(e^{\mathrm{Ah}}{B}_{j,2i+1}\right),G_{2i+2-T_j}=-\frac{h}{3}{B}_{j,2i+2}.$

基于公式(17)-(18)，动力学方程两个周期之间的离散映射关系为： 

$\left(\begin{array}{c}{P}_1{[{x_0}^T,{x_1}^T,{x_2}^T,{x_3}^T,{x_4}^T,···,{x_{{2i}_m+1}}^T,{x_{{2i}_m+2}}^T]}^{T=}\\ Q{[{x_{0-T}}^T,{x_{1-T}}^T,{x_{2-T}}^T,{x_{3-T}}^T,{x_{4-T}}^T,···,{x_{{2i}_m+1-T}}^T,{x_{{2i}_m+2-T}}^T]}^T\\ +P_2{[{x_0}^T,{x_1}^T,{x_2}^T,{x_3}^T,{x_4}^T,···,{x_{{2i}_m+1}}^T,{x_{{2i}_m+2}}^T]}^T\end{array}\right)---\left(19\right)$

其中， 

需要指出的是P₁，P₂，和Q中的分块矩阵的位置由公式(17)-(19)决定。具体的，P₁中的分块矩阵的位置是固定的，如公式(20)所示，P₂，和Q中的分块矩阵的行位置也是固定的，要与P₁对应，而P₂，和Q中的分块矩阵的列位置不是固定的，与对应时滞T_j的大小有关。 

状态转移矩阵可表示为： 

其中，表示矩阵(P₁-P₂)的Moor-Penrose广义逆。 

根据Floquet理论，如果状态转移矩阵Φ的所有特征值的模都小于1，则系统是稳定的，反之，如果状态转移矩阵Φ的任何一个特征值的模大于1，则系统是不稳定的。因此，可以根据Floquet理论画出系统在切削参数空间的稳定边界，即稳定性图谱Lobe图。 

优选地，所述步骤5，具体为： 

改变以上步骤中的螺旋角β和齿间距T_j，重复步骤4，获得新的参数条件下的稳定 性图谱。需要注意的是，设计齿间距时要将铣刀动平衡的约束条件考虑在内，且满足 $T=\underset{j=1}{\overset{N}{\Sigma }}{T}_j.$

优选地，所述步骤6具体为：根据不同齿间距和螺旋角条件下得到的稳定性图谱，根据具体加工时所采用的转速范围，以获得最大加工效率为优化目标，获得优化后的立铣刀关键几何参数即齿间距和螺旋角。 

下面结合具体加工实例说明本发明的具体实施方案，实例参数引自文献1Altintas et a1.，1999，Analytical Stability Prediction and Design of Variable Pitch Cutters，ASME Transactions on Journal of Mechanical Science and Engineering。如图3所示，铣刀直径D＝1.905×10^-2m，螺旋角β＝30/180*π，齿数N＝4，齿间距分别为70°-110°-70°-110°(因此T₁＝7/36T,T₂＝11/36T,T₃＝7/36T,T₄＝11/36T)，径向切深比a/D＝0.5，模态参数如表1所示： 

表1不等距立铣刀实例模态试验参数 

模态振型>f_n(Hz)>ξ(％)>m(kg)>x>563.6>0.055801>1.4986>y>516.21>0.025004>1.199>

其中，f_n表示固有频率，ζ表示相对阻尼比，m表示模态质量。 

将已知参数代入发明内容中的步骤1-步骤4，在步骤2中，总的时间离散步数为m＝72，每一个齿间距的离散步数分别为m₁＝14，m₂＝22，m₁₃＝14，m₄＝22。 

在切削参数空间，横坐标为主轴转速，纵坐标为轴向切削深度，不等齿距立铣刀加工系统的稳定性图谱如图4所示。两种不同方法的结果对比证明了本发明采用的GRK法的有效性与精确性。 

改变立铣刀的关键几何参数齿间距和螺旋角会得到不同的加工参数空间的稳定性图谱，如图5所示。实例所采用的铣刀直径为D＝1.905×10^-2m，齿数为N＝4，考虑到排屑的畅通性，齿间角不小于55°。为了图5的清晰性，图中只绘出了四种齿间距条件下对应的Lobe图，齿间角分别为90°-90°-90°-90°、80°-100°-80°-100°、70°-110°-70°-110°、60°-120°-60°-120°，螺旋角从控制变量的角度出发设为0°。 

选定要加工工件的转速，如选取5000rpm-6000rpm作为加工转速，则最大稳定切削轴向切深对应的加工效率最高，因此根据图5中的Lobe图，选取齿间角为70°-110°-70°-110°的不等齿距立铣刀能够获得最大的加工效率。从图中可以看出，齿 间角对Lobe图的影响主要是使Lobe图沿横轴平移，对Lobe图沿纵轴变化的影响不大。 

需要指出的是： 

(1)初步选定70°-110°-70°-110°作为齿间角后，可以对此齿间角进行小范围调整，对比不同的Lobe图，找到更为高效稳定加工的齿间角。 

(2)图5所示的Lobe图是在径向切深比a/D＝0.5条件下得出的，如果径向切深比a/D发生变化，Lobe图主要是沿纵轴发生变化，沿横轴的变化不大，因此基本不影响指定转速范围内的齿间角选择。以70°-110°-70°-110°不等齿距立铣刀为例，其在不同径向切深比条件下的稳定性图谱如图6所示。 

选定齿间角之后，就再选定螺旋角，还是利用GRK法画出不同螺旋角下的Lobe图，通过分析Lobe图随螺旋角的变化趋势选定加工效率最高的螺旋角。从切削性能考虑，螺旋角的选择范围一般在30°-60°之间，常见的螺旋角为30°、45°、60°三种。研究证明，螺旋角对Lobe图的影响很小，尤其是在径向切深比较大时其影响几乎可以忽略，仍然以70°-110°-70°-110°不等齿距立铣刀为例，其在不同螺旋角条件下的稳定性图谱如图7所示，此时a/D＝0.5。因此此时的螺旋角可以在30°-60°范围内任意指定。 

需要指出的是，如果指定的加工转速范围内存在倍周期分岔，那么适当增加螺旋角会有利于减小倍周期分岔，增加稳定切削轴向切深。以70°-110°-70°-110°不等齿距立铣刀为例，选用小径向切深比a/D＝0.05，画出加工系统在螺旋角分别为0°、30°、45°、60°条件下的稳定性图谱，结果如图8所示。螺旋角为0°时，在主轴转速为6400rpm-6700rpm范围内出现了倍周期分岔，随着螺旋角的增加，倍周期分岔也越来越小，螺旋角为60°时，倍周期分岔完全消失，但是出现了孤岛效应。所以在此主轴转速范围内，可以选取45°的螺旋角，稳定加工轴向切深由9mm提高至13mm，加工效率得到极大的提高。 

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。