CN109740264A - 一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法 - Google Patents

Abstract

本发明提供一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法，将时滞周期等分成若干离散小区间，然后在每个小区间上采用三次牛顿插值多项式和分段三次埃尔米特插值多项式分别插值时滞微分方程中的状态项和时滞项，采用节点上函数值相同及节点上相同的导数值来逼近所需项，减少了计算方法的拟合误差。本发明可以有效提高计算精度和计算效率，从而更加准确高效地预测铣削稳定域，为切削技术人员选择合理地切削参数来提高零件的表面质量提供理论依据。

Description

一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法

技术领域

本发明属于先进制造技术领域，尤其涉及一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法。

背景技术

铣削作为一种高效的机械加工，具有加工精度高、材料去除率大和加工成本低等优点，因而广泛应用于航空、航天、船舶和汽车关键零部件加工等领域。颤振常常存在于铣削加工过程中，是造成加工过程不稳定的主要因素。此外，铣削过程中的颤振不但会影响加工零件的表面质量，而且会降低数控机床的寿命。通过铣削过程颤振稳定性预测，可选取无颤振铣削参数组合，以达到避免颤振的发生，提高加工效率，降低加工成本的目的。因此，高精度、高效率的稳定性预测方法对充分发挥高端数控设备的高性能加工特性具有重要意义。

目前，文献“Insperger T,Stépán G.Updated semi-discretization method forperiodic delay-differential equations with discrete delay[J].InternationalJournal for Numerical Methods in Engineering,2004,61(1):117-141.”中介绍的半离散法(SDM)，该方法具有较高的计算精度，由于仅对时滞项进行离散，致使该方法的计算效率比较低。文献“Ding Y,Zhu LM,Zhang X J,et al.A full-discretization method forprediction of milling stability[J].International Journal of Machine Tools andManufacture,2010,50(5):502-509.”提出了全离散法(FDM)，该方法离散状态项、时滞项和系数矩阵，相比SDM，该方法计算效率得到相对较大提高。然后文献“Quo Q,SunY,JiangY.Onthe accurate calculation of milling stability limits using third-order full-discretization method[J].International Journal of Machine Tools andManufacture,2012,62:61-66.”提出了一种用三阶全离散法精确计算铣削稳定性极限的方法，该方法较一阶、二阶的全离散法在计算精度上有所提高，但计算方法的效率还是很低。文献“Ozoegwu C G,Omenyi S N,Ofochebe S M.Hyper-third order full-discretization methods in milling stability prediction[J].InternationalJournal of Machine Tools and Manufacture,2015,92:1-9.”推广到了更高阶(四阶、五阶)全离散法，阶数越高，消耗计算时间越多，而且时滞项还是采用一阶线性插值，导致计算精度提高相对较少。

目前申请号为：“201510068454.7”的专利，在获取稳定性叶瓣图时，传递矩阵需要离散间隔数m次迭代，降低了计算的效率；申请号为：“201510067259.2”的发明专利，只采用节点上函数值相同来拟合时滞项，计算精度不高。

发明内容

发明目的：为解决上述计算精度和计算效率不能同时提高，并且计算量大的问题，本发明提供了一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法。

技术方案：本发明提供一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法，该方法具体包括如下步骤；

步骤1：将考虑再生颤振铣削动力学的时滞微分方程变换为空间状态形式，得到该时滞微分方程的空间状态方程；

步骤2；将刀齿切削周期等分成若干个小间隔，在任意一个小间隔的时间长度内对上述时滞微分方程的空间状态方程进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程；

步骤3；利用插值多项式算法逼近上述的状态项、时滞项、周期系数项，从而得到状态项、时滞项、周期系数项的方程表达式；并将状态项、时滞项、周期系数项的方程表达式代入步骤2中的状态项、时滞项、周期系数项三者之间的关系方程中，并利用迭代法计算得到相邻周期刀齿响应之间的映射关系；根据所述的映射关系构建相邻周期刀齿响应的传递矩阵Ф；

步骤4：计算步骤3中的相邻周期刀齿响应传递矩阵Ф的特征值，判定铣削系统的稳定性，具体的判定如公式1所示：

其中λ()的含义为求特征值。

进一步的，在单自由度系统中所述步骤1的具体方法为：

步骤1.1：构建考虑再生颤振的单自由度铣削动力学的时滞微分方程，如公式2所示：

其中ξ、m_t、ω_n分别为刀具的阻尼比、质量和固有频率；h(t)表示周期为T的切削力系数；q(t)为刀具的坐标，为q(t)的一阶导数，为q(t)的二阶导数；w为轴向切深；

步骤1.2：令X(t)＝[q(t) p(t)]^T，其中将公式2转换成空间状态形式如公式3所示：

其中A₀为常数矩阵 为X(t)的一阶导数；B(t)表示周期为T的系数矩阵，

进一步的，所述步骤2的具体方法为：将刀齿切削周期T等分成m个小间隔，则具有y＝m+1个节点，步长将公式2在任意一个时间区间[t_K，t_K+1]内进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程，如公式4所示，其中t_K表示第K个节点的时刻；所述K＝1,2,…,y；

其中X(δ)为状态项、X(δ-T)为时滞项、B(δ)为周期系数项，X_k+1表示时刻t_K+1的响应；X_k表示时刻t_K的响应；δ为时间变量，0≤δ≤τ或者t_K≤δ≤t_K+1。

进一步的，所述步骤3中确定状态项、时滞项、周期系数项的方程表达式的具体方法为：

在当前周期的时间区域[t_K-2，t_K+1]中，利用牛顿插值多项式算法将该时间区域中的四个时刻t_K-2、t_K-1、t_K、t_K+1以及这四个时刻的响应X_K-2、X_K-1、X_K、X_K+1逼近状态项，K>2，得到状态项X(δ)的方程表达式为：

X(δ)＝a₁X_K+1+b₁X_K+c₁X_K-1+d₁X_K-2 (5)

其中，a₁、b₁、b₁、d₁为系数，具体的表示如公式6所示；

在上一个周期的时间区域[t_K，-m，t_K+2，-m]中，利用分段三次埃米尔特插值多项式算法将三个时刻t_K，-m、t_K+1，-m、t_K+2，-m以及这三个时刻的响应X_K，-m、X_K+1，-m、X_K+2，-m逼近时滞项，得到时滞项X(δ-T)的方程表达式如公式7所示，所述t_K，-m、t_K+1，-m、t_K+2，-m分别表示上一个周期中节点K、K+1、K+2的时刻：

X(δ-T)＝a₂X_K，-m+b₂X_K+1，-m+c₂X_K+2，-m (7)

其中a₂、b₂、c₂为系数，具体表达式如公式8所示；

用线性插值算法逼近周期系数项，得到周期系数项B(δ)的方程表达式为：

其中B_K+1表示在当前周期中第K+1个节点的周期系数；B_K表示当前周期中第K个节点的周期系数。

进一步的，步骤3中的构建相邻周期刀齿响应的传递矩阵Ф具体方法如下：

步骤A：将公式5、7、9代入状态项、时滞项、周期系数项三者之间的关系方程中，整理得到公式10；

(H_K，1—I)X_K+1+(H₀+H_k，0)X_K+H_k，-1X_K-1+H_k，-2X_K-2

＝H_k，mX_K，-m+H_k，m-1X_K+1，-m+H_k，m-2X_K+2，-m (10)

其中，I为单位矩阵，H_k，1为当前周期中第K+1个节点响应的切削力系数；H_k，0为当前周期中第K个节点响应的切削力系数，H_k，-1为当前周期中第K-1个节点响应的切削力系数；H_k，-2为当前周期中第K-2个节点响应的切削力系数；H_k，m-2为上一个周期中第K-2个节点响应的切削力系数；H_k，m-1为上一个周期中第K-1个节点响应的切削力系数；H_k，m为上一个周期中第K个节点响应的切削力系数；

其中，

步骤B：将上述公式10进行m次的迭代计算，得到相邻周期刀齿响应之间的映射关系，用矩阵表示如公式11所示：

其中X_m为当前周期的中时刻t_m的响应；X_m，-m为上一个周期中时刻t_m，-m的响应；D₁和D₂表达式如下所示：

步骤C：构建相邻周期刀齿响应的传递矩阵Ф，如公式12所示：

Φ＝(D₁)^-1D₂ (12)。

进一步的，上述用分段三次埃米尔特插值多项式算法逼近时滞项，得到时滞项X(δ-T)的方程表达式的具体方法为：

步骤S1:在上一个周期的每个时间区域[t_K，-m，t_K+1，-m]内，构造分段三次埃尔米特插值多项式，得到公式13：

其中，为响应X_K，-m的一阶导数，为响应X_K+1，-m的一阶导数；系数ε₀(δ)、ε₁(δ)、β₀(δ)、β₁具体的表述如公式14所示：

当δ＝t_K，-m时，

当δ＝t_K+1，-m时，

步骤S2：由于令τ＝t_K+1，-m—t_K，-m；并将上述公式15和公式16代入公式13中；得到X(δ-T)的表达式如公式17所示：

进一步的，在两自由度系统中所述步骤1的具体方法为：

构建两自由度再生颤振铣削动力学的时滞微分方程如公式18所示：

其中q₁(t)和q₂(t)分别表示在两自由度系统中x，y方向上刀具的坐标，h_xx(t)与h_yy(t)表示在两自由度系统中x，y方向上周期为T的切削力系数；h_xy(t)表示x方向与y方向之间的耦合作用，h_yx(t)表示y方向与x方向之间的耦合作用；ξ、m_t、ω_n分别为刀具的阻尼比、质量和固有频率；

令x₁(t)＝[q₁(t) q₂(t)]，将公式18变换为空间状态形式，如公式19所示

其中，分别表示为q₁(t)，q₂(t)的一阶导数，为x₁(t)的一阶导数，A₁为常数矩阵，B₁(t)表示周期为T的系数矩阵，具体A₁和B₁(t)的表示如公式20、21所示：A₁为常数矩阵，B₁(t)表示周期为T的系数矩阵，具体A₁和B₁(t)的表示如公式20、21所示：

其中w为轴向切深。

有益效果：本发明采用分段三次埃尔米特插值多项式来插值空间状态动力学方程中的时滞项，采用节点上函数值相同及节点上相同的导数值来逼近所需项，减少了计算方法的拟合误差，从而提高了预测方法的精度；同时在求解稳定性叶瓣图的过程中，构建传递矩阵单乘，减少传递矩阵计算过程中的迭代次数，显著地提高计算效率。

附图说明

图1为本发明单自由度时径向浸入比为0.1时的稳定性叶瓣图；

图2为本发明单自由度时径向浸入比为1时的稳定性叶瓣图；

图3为本发明两自由度时径向浸入比为0.1时的稳定性叶瓣图；

图4为本发明两自由度时径向浸入比为1时的稳定性叶瓣图；

图5为本发明轴向切深w为0.2mm，主轴转速Ω为5000rmp的收敛图；

图6为本发明轴向切深w为0.5mm，主轴转速Ω为5000rmp的收敛图；

图7为本发明的整体流程图。

具体实施方式

构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

如图7所示，本发明提供一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法，该方法具体包括如下步骤；

步骤1：构建考虑单自由度再生颤振铣削动力学的时滞微分方程，如公式22所示：

其中，ζ、ω_n、m_t、w和q(t)分别是刀具的阻尼比、固有频率、质量、轴向切深和坐标，为q(t)的一阶导数，为q(t)的二阶导数，h(t)表示周期为T的切削力系数，h(t)的表达式如公式23所示。本实施例中周期T与时滞量T_s相等，T_s＝60/NΩ，N为铣刀齿数，本实施例中N的值为2，Ω为铣刀主轴转速，铣刀的一阶固有频率为922*2π，相对阻尼为0.11，质量为0.03993。

其中，s_t和s_n分别是切向和法向的切削力系数，本实施例中分别取6*10⁸和2*10⁸；φ_j(t)是第j刀齿的角位移，表达式为φ_j(t)＝(2πΩ/60)t+2π(j-1)/N；函数g(φ_j(t))表达式为：

其中φ_st和φ_ex分别为铣刀切入角和切出角。当逆铣时，φ_st＝0，φ_ex＝arcos(1-2a/D)；当顺铣时，φ_st＝arcos(2a/D-1)和φ_ex＝π，a/D是径向浸入比。

将公式22变换为空间状态形式，并令X(t)＝[q(t) p(t)]^T， 得到公式23：

其中A₀为常数矩阵B(t)表示周期为T的系数矩阵，

步骤2：将刀齿切削周期T等分成m个小间隔，本实施例中m取值为40；则具有y＝m+1个节点，步长将公式23在任意一个时间区间[t_K，t_K+1]内进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程，如公式24所示，其中t_K表示第K个节点的时刻；所述K＝1,2,…,y；

其中X(δ)为状态项、X(δ-T)为时滞项、B(δ)为周期系数项，X_k+1表示时刻t_K+1的响应；X_k表示时刻t_K的响应；δ为时间变量，0≤δ≤τ或者t_K≤δ≤t_K+1。

步骤3：在当前周期的时间区域[t_K-2，t_K+1]中，利用牛顿插值多项式算法将该时间区域中的四个时刻t_K-2、t_K-1、t_K、t_K-2以及这四个时刻的响应X_K-2、X_K-1、X_K、X_K+1拟合状态项，K>2，得到状态项的X(δ)的方程表达式为：

X(δ)＝a₁X_K+1+b₁X_K+c₁X_K-1+d₁X_K-2 (25)

其中，a₁、b₁、b₁、d₁为系数，具体的表示如公式26所示；

在上一个周期的时间区域[t_K，-m，t_K+2，-m]中，利用分段三次埃米尔特插值多项式算法将三个时刻t_K，-m、t_K+1，-m、t_K+2，-m以及这三个时刻的响应X_K，-m、X_K+1，-m、X_K+2，-m逼近时滞项，得到时滞项X(δ-T)的方程表达式如公式27所示，所述t_K，-m、t_K+1，-m、t_K+2，-m分别表示上一个周期中节点K、K+1、K+2的时间：

X(δ-T)＝a₂X_K，-m+b₂X_K+1，-m+c₂X_K+2，-m (27)

其中a₂、b₂、c₂为系数，具体表达式如公式28所示；

用线性插值算法逼近周期系数项，得到周期系数项B(δ)的方程表达式为：

其中B_K+1表示在当前周期中第K+1个节点的周期系数；B_K表示当前周期中第K个节点的周期系数。

步骤4：将上述公式25、27、29代入公式24中，整理得到公式30：

(H_k，1—I)X_K+1+(H₀+H_k，0)X_K+H_k，-1X_K-1+H_k，-2X_K-2

＝H_k，mX_K，-m+H_k，m-1X_K+1，-m+H_k，m-2X_K+2，-m (30)

其中，I为单位矩阵，H_k，1为当前周期中第K+1个节点响应的切削力系数；H_k，0为当前周期中第K个节点响应的切削力系数，H_k，-1为当前周期中第K-1个节点响应的切削力系数；H_k，-2为当前周期中第K-2个节点响应的切削力系数；H_k，m-2为上一个周期中第K-2个节点响应的切削力系数；H_k，m-1为上一个周期中第K-1个节点响应的切削力系数；H_k，m为上一个周期中第K个节点响应的切削力系数；

其中，

步骤B：将上述公式30进行m次的迭代计算，得到相邻周期刀齿响应之间的映射关系，用矩阵表示如公式31所示：

其中X_m为当前周期的中时刻t_m的响应；X_m，-m为上一个周期中时刻t_m，-m的响应；D₁和D₂表达式如下所示：

其中，H_1，0表示当前周期中第一个节点响应的切削力系数；H_1，1表示当前周期中第1+1个节点即第二个节点响应的切削力系数。

构建相邻周期刀齿响应的传递矩阵Ф，即Floquet转移矩阵，如公式32所示：

Φ＝(D₁)^-1D₂ (32)

步骤6：根据上述相邻周期刀齿响应传递矩阵Ф的特征值，模判定铣削系统的稳定性，当Ф最大特征值的模小于1时，该铣削动力学系统处于稳定状态；Ф最大特征值的模等于1时，该铣削动力学系统临界稳定状态，当Ф最大特征值的模大于1时，该铣削动力学系统处于不稳定状态；具体的判定如公式33所示：

其中λ()的含义为求特征值。

上述步骤3中用分段三次埃米尔特插值多项式算法逼近时滞项，得到时滞项X(δ-T)的方程表达式的具体方法为：

步骤3.1:在上一个周期的每个时间区域[t_K，-m，t_K+1，-m-m]内，构造分段三次埃尔米特插值多项式，得到公式34：

其中，为响应X_K，-m的一阶导数，为响应X_K+1，-m的一阶导数；其中系数ε₀(δ)、ε₁(δ)、β₀(δ)、β₁具体的表述如公式35所示：

当δ＝t_K，-m时，

当δ＝t_K+1，-m时，

步骤3.2：由于令τ＝t_K+1，-m—t_K，-m并将上述公式36和公式37代入公式34中；得到X(δ-T)的表达式如公式38所示：

在两自由度的系统中，该预测方法的步骤1的具体方法为：

建立两自由度再生颤振铣削动力学的时滞微分方程如公式39所示：

其中q₁(t)和q₂(t)分别表示在两自由度系统中x，y方向的坐标，h_xx(t)与h_yy(t)表示在两自由度系统中x，y方向的切削力系数；h_xy(t)表示x方向与y方向之间的耦合作用，h_yx(t)表示y方向与x方向之间的耦合作用；ξ、m_t、ω_n分别为刀具的阻尼比、质量和固有频率；

令x₁(t)＝[q₁(t) q₂(t)]，将公式18变换为空间状态形式，如公式44所示

其中，分别表示为q₁(t)，q₂(t)的一阶导数，为x₁(t)的一阶导数，A₁为常数矩阵，B₁(t)表示周期为T的系数矩阵，具体A₁和B₁(t)的表示如公式45、46所示：其中，

将两自由度的动力学方程变换为空间状态形式，再进行上述步骤2至步骤6的操作后得到两自由度的稳定域的准确预测。

本实施例利用Matlab软件进行编程，画出单自由度和两自由度的稳定性叶瓣图及收敛图，通过稳定性图来预测铣削过程中地稳定性，选取不同地径向浸入比，分别取0.1和1得到的稳定性叶瓣图如图1、2、3和4所示及选取不同轴向切深w，分别取0.2mm和0.5mm得到的收敛图如图5和6所示。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

Claims (7)

1.一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法，其特征在于，具体包括如下步骤：

步骤1：将考虑再生颤振铣削动力学的时滞微分方程变换为空间状态形式，得到该时滞微分方程的空间状态方程；

步骤2：将刀齿切削周期等分成若干个小间隔，在任意一个小间隔的时间长度内对上述时滞微分方程的空间状态方程进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程；

步骤3：利用插值多项式算法逼近上述的状态项、时滞项、周期系数项，从而得到状态项、时滞项、周期系数项的方程表达式；并将状态项、时滞项、周期系数项的方程表达式代入步骤2中的状态项、时滞项、周期系数项三者之间的关系方程中，并利用迭代法计算得到相邻周期刀齿响应之间的映射关系；根据所述的映射关系构建相邻周期刀齿响应的传递矩阵Φ；

步骤4：计算步骤3中的相邻周期刀齿响应传递矩阵Φ的特征值，判定铣削系统的稳定性，具体的判定如公式1所示：

其中λ()的含义为求特征值。

2.基于权利要求1所述的预测方法，其特征在于，在单自由度系统中所述步骤1的具体方法为：

步骤1.1：构建考虑再生颤振的单自由度铣削动力学的时滞微分方程，如公式2所示：

其中ξ、m_t、ω_n分别为刀具的阻尼比、质量和固有频率；h(t)表示周期为T的切削力系数；q(t)为刀具的坐标，为q(t)的一阶导数，为q(t)的二阶导数；w为轴向切深；

步骤1.2：令X(t)＝[q(t)p(t)]^T，其中将公式2转换成空间状态形式如公式3所示：

其中A₀为常数矩阵 为X(t)的一阶导数；B(t)表示周期为T的系数矩阵，

3.基于权利要求2所述的预测方法，其特征在于，所述步骤2的具体方法为：将刀齿切削周期T等分成m个小间隔，则具有y＝m+1个节点，步长将公式2在任意一个时间区间[t_K，t_K+1]内进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程，如公式4所示，其中t_K表示第K个节点的时刻；所述K＝1，2，...，y；

其中X(δ)为状态项、X(δ-T)为时滞项、B(δ)为周期系数项，X_k+1表示时刻t_K+1的响应；X_k表示时刻t_K的响应；δ为时间变量，0≤δ≤τ或者t_K≤δ≤t_K+1。

4.基于权利要求3所述的预测方法，其特征在于，所述步骤3中确定状态项、时滞项、周期系数项的方程表达式的具体方法为：

在当前周期的时间区域[t_K-2，t_K+1]中，利用牛顿插值多项式算法将该时间区域中的四个时刻t_K-2、t_K-1、t_K、t_K+1以及这四个时刻的响应X_K-2、X_K-1、X_K、X_K+1逼近状态项，K＞2，得到状态项X(δ)的方程表达式为：

X(δ)＝a₁X_K+1+b₁X_K+c₁X_K-1+d₁X_K-2 (5)

其中，a₁、b₁、b₁、d₁为系数，具体的表示如公式6所示；

在上一个周期的时间区域[t_K，-m，t_K+2，-m]中，利用分段三次埃米尔特插值多项式算法将三个时刻t_K，-m、t_K+1，-m、t_K+2，-m以及这三个时刻的响应X_K，-m、X_K+1，-m、X_K+2，-m逼近时滞项，得到时滞项X(δ-T)的方程表达式如公式7所示，所述t_K，-m、t_K+1，-m、t_K+2，-m分别表示上一个周期中节点K、K+1、K+2的时刻：

X(δ-T)＝a₂X_K，-m+b₂X_K+1，-m+c₂X_K+2，-m (7)

其中a₂、b₂、c₂为系数，具体表达式如公式8所示；

用线性插值算法逼近周期系数项，得到周期系数项B(δ)的方程表达式为：

其中B_K+1表示在当前周期中第K+1个节点的周期系数；B_K表示当前周期中第K个节点的周期系数。

5.基于权利要求4所述的预测方法，其特征在于，步骤3中的构建相邻周期刀齿响应的传递矩阵Φ具体方法如下：

步骤A：将公式5、7、9代入状态项、时滞项、周期系数项三者之间的关系方程中，整理得到公式10；

(H_k，1-I)X_K+1+(H₀+H_k，0)X_K+H_k，-1X_K-1+H_k，-2X_K-2

＝H_k，mX_K，-m+H_k，m-₁X_K+1，-m+H_k，m-2X_K+2，-m (10)

其中，I为单位矩阵，H_k，1为当前周期中第K+1个节点响应的切削力系数；H_k，0为当前周期中第K个节点响应的切削力系数，H_k，-1为当前周期中第K-1个节点响应的切削力系数；H_k，-2为当前周期中第K-2个节点响应的切削力系数；H_k，m-2为上一个周期中第K-2个节点响应的切削力系数；H_k，m-1为上一个周期中第K-1个节点响应的切削力系数；H_k，m为上一个周期中第K个节点响应的切削力系数；

其中，

步骤B：将上述公式10进行m次的迭代计算，得到相邻周期刀齿响应之间的映射关系，用矩阵表示如公式11所示：

其中X_m为当前周期的中时刻t_m的响应；X_m，-m为上一个周期中时刻t_m，-_m的响应；D₁和D₂表达式如下所示：

步骤C：构建相邻周期刀齿响应的传递矩阵Φ，如公式12所示：

Φ＝(D₁)^-1D₂ (12)。

6.基于权利要求4所述的预测方法，其特征在于，上述用分段三次埃米尔特插值多项式算法逼近时滞项，得到时滞项X(δ-T)的方程表达式的具体方法为：

步骤S1：在上一个周期的每个时间区域[t_K，-m，t_K+1，-m]内，构造分段三次埃尔米特插值多项式，得到公式13：

其中，为响应X_K，-m的一阶导数，为响应X_K+1，-m的一阶导数；系数ε₀(δ)、ε₁(δ)、β₀(δ)、β₁具体的表述如公式14所示：

当δ＝t_K，-m时，

当δ＝t_K+1，-m时，

步骤S2：由于令τ＝t_K+1，-m-t_K，-m；并将上述公式15和公式16代入公式13中；得到X(δ-T)的表达式如公式17所示：

7.基于权利要求1所述的预测方法，其特征在于，在两自由度系统中所述步骤1的具体方法为：

构建两自由度再生颤振铣削动力学的时滞微分方程如公式18所示：

其中q₁(t)和q₂(t)分别表示在两自由度系统中x，y方向上刀具的坐标，h_xx(t)与h_yy(t)表示在两自由度系统中x，y方向上周期为T的切削力系数；h_xy(t)表示x方向与y方向之间的耦合作用，h_yx(t)表示y方向与x方向之间的耦合作用；ξ、m_t、ω_n分别为刀具的阻尼比、质量和固有频率；

令将公式18变换为空间状态形式，如公式19所示

其中，分别表示为q₁(t)，q₂(t)的一阶导数，为x₁(t)的一阶导数，A₁为常数矩阵，B₁(t)表示周期为T的系数矩阵，具体A₁和B₁(t)的表示如公式20、21所示：

其中w为轴向切深。