CN110188311A - 基于刀齿切削时程精细积分的高速加工稳定域预测方法 - Google Patents

Abstract

本发明提供了一种基于刀齿切削时程精细积分的高速加工稳定域预测方法，包括：首先构建考虑再生颤振的高速铣削动力学模型；然后离散主轴旋转周期，在每个离散的小间隔上将时滞微分方程中Duhamel积分简化为状态项和时滞项的积分；其次采用分段三次埃尔米特插值多项式拟合时滞项，同时采用多个已知时间节点响应及其导数来逼近所需项，降低了计算方法的拟合误差，从而提高了预测方法的精度；而且在求解稳定性叶瓣图的过程中，运用精细积分算法，避免了计算过程中的矩阵求逆运算，极大的减少了计算方法所消耗时间，提高了计算效率。

Description

基于刀齿切削时程精细积分的高速加工稳定域预测方法

技术领域

本发明属于先进制造技术领域，尤其涉及一种基于刀齿切削时程精细积分的高速加工稳定域预测方法。

背景技术

高速铣削具有在加工过程中保持相对低的切削力、保证高加工精度的同时又能获得高材料去除率的优势，被广泛应用于航空、航天和汽车工业制造的精加工领域。高速铣削的诸多优势的发挥是以高速条件下的无振动稳定切削过程的实现为前提。由于缺乏系统有效的高速切削振动和加工不稳定性控制的理论指导和实用工具，相关企业在制造构件时大多仍沿用经验或试切法确定铣削参数，致使重金购置的高档数控机床低速、低效使用，不能充分发挥高速铣削技术的优越性，造成巨大的浪费。通过对高速铣削过程的稳定域进行预测，科学合理优化铣削参数，以达到高速铣削条件下的无振动稳定加工的目的，对充分发挥数控机床特别是高档机床的加工性能和提高数控加工的质量和效率具有重要意义。因此，高精度、高效的铣削稳定域预测方法显得尤为重要。

目前，文献“Quo Q,SunY,JiangY.On the accurate calculation ofmillingstability limits using third-order full-discretization method[J].International Journal of Machine Tools and Manufacture,2012,62:61-66.”中公开了一种基于三阶全离散法的铣削稳定性预测方法，该方法通过采用三阶牛顿插值多项式逼近Duhamel积分，具有较好的计算效率，但计算精度很差。然后文献“Liu Y,Zhang D,WuB.An efficient full-discretization method for prediction of milling stability[J].International Journal of Machine Tools and Manufacture,2012,63:44-48.”中公开了一种基于全离散法的高效铣削稳定性预测方法，该方法具有较高的计算精度，但计算的效率较低。

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

发明内容

发明目的：为解决上述稳定性预测方法不能同时达到实际要求的预测精度和效率，本发明提供了一种基于刀齿切削时程精细积分的高速加工稳定域预测方法。

技术方案：本发明提供一种基于刀齿切削时程精细积分的高速加工稳定域预测方法，该方法包括如下步骤：

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

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

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

步骤4：利用精细积分算法对传递矩阵方程进一步求解，得到求解后的传递矩阵方程；

步骤5：根据求解后的传递矩阵方程，得到相邻周期刀齿响应之间的映射关系；并构建相邻周期刀齿响应的传递矩阵Ф；

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

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

进一步的，所述得到时滞微分方程的空间状态方程的方法如下所示：

所述考虑再生颤振铣削动力学的时滞微分方程如公式2所示：

其中，M、C、q(t)、w和K_c(t)分别表示刀具的模态质量、模态阻尼、模态刚度矩阵、模态坐标、轴向切削深度和周期系数矩阵；T为时滞量，T＝60/NΩ，N为铣刀齿数，Ω为主轴转速，为q(t)的一阶导数；为q(t)的二阶导数；

令令x(t)＝[q(t)p(t)]^T，并通过空间变换，将公式2变换成如下形式：

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

将周期T等分成m个小间隔，则具有m+1个节点，步长τ＝T/m；在任意一个时间区间[t_K,t_K+1]内对公式2进行积分计算，得到公式4，其中t_K表示第K个节点的时刻；所述K＝1,2,…,m+1；

其中，x_K+1表示时刻t_K+1的响应；x_K表示时刻t_K的响应；δ为时间变量，0≤δ≤τ或者t_K≤δ≤t_K+1；

令B(δ)x(δ)＝f(δ)为状态项，B(δ)x(δ-T)＝g(δ-T)为时滞项，得到状态项与时滞项二者之间的关系方程，如公式5：

进一步，所述得到传递矩阵方程的具体方法为：

步骤3.1：通过构建三次牛顿插值多项式逼近状态项f(δ)，得到f(δ)的表达式，具体方法如下所示：

在时间区域[t_K-2,t_K+1]中，利用牛顿插值多项式算法将该时间区域中的四个时刻t_K-2、t_K-1、t_K、t_K+1以及这四个时刻的响应f_K-2、f_K-1、f_K、f_K+1对状态项进行插值计算，得到状态项f(δ)的方程表达式为：

f(δ)＝a₁f_K+1+b₁f_K+c₁f_K-1+d₁f_K-2 (6)

其中，a₁、b₁、b₁、d₁如下所示；

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

g(δ-T)＝a₂g_K-m+b₂g_K+1-m+c₂g_K+2-m (7)

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

步骤3.2：将公式6和公式7代入公式5中，得到公式8，即传递矩阵方程：

其中

进一步的，所述利用精细积分算法对传递矩阵方程进一步求解具体为:利用精细积分对进行求解，具体的方法为：

步骤4.1：对任意一个步长为τ的间隔进行精细划分，划分为2^p段精细区段，则任意一段精细则：

对T1进行泰勒级数展开：

其中Ta为增量，I为单位矩阵；

将公式10，代入公式9,得到公式11：

经过P次迭代计算，得到：

把公式12代入到公式8，中得到：

H_K,-2x_K-2+H_K,-1x_K-1+(H_K,0+T1)x_K+(H_K,1-I)＝H_K,mx_K-m+H_K,m-1x_K+1-m+H_K,m-2x_K+2-m (13)

其中H_K,1＝(h_K,1+i_K,1)B_K+1；

H_K，0＝(h_K，0+i_K，0)B_K；

H_K，-1＝(h_K，-1+i_K，-1)B_K-1；

H_K，-2＝(h_K，-2+i_K，-2)B_K-2；

H_K，m＝(h_K，m+i_K，m)B_K；

H_K，m-1＝(h_K，m-1+i_K，m-1)B_K+1；

H_K，m-2＝(h_K，m-2+i_K，m-2)B_K+2；

其中B_K+1、B_K、B_K-1、B_K-2、B_K+2为节点K+1、K、K-1、K-2、K+2处的周期系数。

进一步的，所述构建相邻周期刀齿响应的传递矩阵Φ的具体方法为：

对求解后的传递矩阵方程进行若干次迭代计算，得到相邻相邻周期刀齿响应之间的映射关系如下所示：

其中

根据映射关系构建相邻周期刀齿响应传递矩阵Φ：

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

进一步的，所述利用精细积分算法对传递矩阵方程进一步求解具体为利用精细积分对进行求解，具体的方法为：

对任意一个步长为τ的间隔进行精细划分，得到2^p段精细区段，则任意一段精细 利用Q项帕德级数近似举证指数得到：

其中：

在使用精细积分时，只计算和存储增量部分，故将N_Q(A₀σ)和D_Q(A₀σ)表示为如下形式：

其中

则

则

将公式24代入公式23，得到T1的表达式，并将T1代入公式8中，得到公式13。

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

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

其中，为响应g_K-m的一阶导数，为响应g_K+1-m的一阶导数；系数μ₀(δ)、μ₁(δ)、如下所示：

当δ＝t_K-m时，

当δ＝t_K+1-m时，

步骤S2：由于令τ＝t_K+1-m-t_K-m；并将上述公式27、公式26和μ₀(δ)、μ₁(δ)、代入公式25中；得到g(δ-T)的表达式如公式28所示：

有益效果：本发明将空间状态动力学微分方程中Duhamel积分简化为状态项和时滞项的积分，通过三次牛顿多项式和分段三次埃尔米特多项式分别逼近状态项和时滞项，减少计算方法的拟合误差，从而提高了本方法的预测精度；同时在求解稳定性叶瓣图的过程中，运用精细积分法，避免计算过程中的矩阵求逆运算，在不牺牲计算方法的精度的前提下减少计算方法的时间，极大的提高计算效率。经测试，当离散切削周期为40段时，本发明方法和背景技术方法仿真需要的时间分别为118秒、214秒和492秒，效率分别提高了45％和76％；当离散切削周期为80段时，本发明方法和背景技术方法仿真需要的时间分别为304秒、572秒和1117秒，效率分别提高了47％和73％；当离散切削周期为120段时，本发明方法和背景技术方法仿真需要的时间分别为590秒、1177秒和2043秒，效率分别提高了50％和71％。

附图说明

图1为本发明的流程图；

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

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

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

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

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

图7为本发明轴向切深w为0.7mm，主轴转速Ω为5000rmp的收敛图。

具体实施方式

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

如图1所示，本实施例将考虑再生颤振铣削动力学的时滞微分方程变换为空间状态形式，得到该时滞微分方程的空间状态方程；将主轴转速离散成K_n个点，将轴向切深离散成K_w个点；令R＝0,1,2,…,K_n；令J＝0,1,2,…K_w；将点(R，J)的刀齿切削周期等分成若干个小间隔；任意一个小间隔的时间长度内对上述时滞微分方程的空间状态方程进行积分计算，得到状态项、时滞项二者之间的关系方程；利用插值多项式算法逼近上述的状态项、时滞项，从而得到状态项、时滞项方程表达式；并将状态项、时滞项的方程表达式代入状态项、时滞项二者之间的关系方程中得到传递矩阵方程，所述传递矩阵方程为用于构建相邻周期刀齿响应的传递矩阵的方程；利用精细积分算法对传递矩阵方程进一步求解，得到矩阵系数D₁和D₂；根据矩阵系数D₁和D₂，得到相邻周期刀齿响应之间的映射关系；并构建相邻周期刀齿响应的传递矩阵Ф。

计算相邻周期刀齿响应传递矩阵Ф的特征值，根据Floquet理论，判定铣削系统的稳定性。

本实施例具体的方法为：

构建考虑再生颤振的动力学模型，其动力学模型可由式(29)描述，此模型可适用双自由度动力学模型：

其中，M、C、q(t)、w和K_c(t)分别表示刀具的模态质量、模态阻尼、模态刚度矩阵、模态坐标、轴向切削深度和周期系数矩阵；T为时滞量，T＝60/NΩ，N为铣刀齿数，单位为rpm；为q(t)的一阶导数；为q(t)的二阶导数。

令和x(t)＝[q(t)p(t)]^T，通过变换，式(29)变换为如下的空间状态形式：

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

将刀齿切削周期T等分为m个小区间，则具有m+1个节点，步长τ＝T/m；在任意一个时间区间[t_K,t_K+1]内对公式30进行积分计算，得到公式31，其中t_K表示第K个节点的时刻；因为本实施例中将刀齿切削周期T等分为m个小区间，在本实施例中第一个节点和最后一个节点为同一个位置；

令B(δ)x(δ)＝f(δ)为状态项，B(δ)x(δ-T)＝g(δ-T)式(31)可表示为

通过构建三次牛顿插值多项式来逼近步骤式(32)的状态项f(δ)和埃尔米特插值多项式来逼近步骤中式(32)的时滞项g(δ-T)，具体过程如下：

在间隔[t_K-2,t_K+1]中，将该间隔上的四个节点t_K-2、t_K-1、t_K、t_K+1及对应的响应f_K-2、f_K-1、f_K、f_K+1对状态项f(δ)进行插值，表示如下：

f(δ)＝a₁f_K+1+b₁f_K+c₁f_K-1+d₁f_K-2 (33)

其中

在间隔[t_K-m，t_K+2-m]中，将该间隔上的三个节点的时刻t_K-m、t_K+1-m、t_K+2-m及对应的响应g_K-m、g_K+1-m、g_K+2-m对时滞项g(δ-T)进行插值，表示如下：

g(δ-T)＝a₂g_K-m+b₂g_K+1-m+c₂g_K+2-m (34)

其中

将式(34)和(33)带入式(32)，整理得：

其中

其中将基本区段τ进行精细划分。记精细区段为σ＝τ/2^P，则

一般取P＝20，精细区段σ就已经是非常小的区段。当σ非常小时，T1可以采用泰勒级数展开有限项进行近似；

将式(37)带入式(36)得：

因此，增量Ta可以通过P次迭代得到，具体算法如下：

for(i＝0；i＜P；i++)

将式(39)代入到式(35)得：

H_K，-2X_K-2+H_K，-1X_K-1+(H_K，0+T1)x_K+(H_K，1-I)＝H_K，mX_K-m+H_K，m-1X_K+1-m+H_K，m-2X_K+2-m (40)

其中

H_K，1＝(h_K，1+i_K，1)B_K+1

H_K，0＝(h_K，0+i_K，0)B_K

H_K，-1＝(h_K，-1+i_K，-1)B_K-1

H_K，-2＝(h_K，-2+i_K，-2)B_K-2

H_K，m＝(h_K，m+i_K，m)B_K

H_K，m-1＝(h_K，m-1+i_K，m-1)B_K+1

H_K，m-2＝(h_K，m-2+i_K，m-2)B_K+2

根据公式(40)，求得相邻周期刀齿响应之间得映射关系，通过矩阵表示如下：

其中

铣削动力学系统的相邻周期刀齿响应传递矩阵Ф表示为

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

计算相邻周期刀齿响应的传递矩阵Ф的特征值，根据Floquet理论，判定铣削系统的稳定性，具体判定准则如下：

本实施例所述运用基于埃尔米特插值多项式的方法来拟合式(32)的时滞项g(δ-T)。具体过程如下：

在每个小间隔[t_k-m,t_k+1-m]上构造分段三次埃尔米特插值多项式

其中

由求导公式可知，g(t)在δ＝t_k-m和δ＝t_k+1-m的导数可表示为

令τ＝t_K+1-m—t_K-m，将式(45)、(46)带入式(44)得

本实施例中的求解也可以采用基于帕德级数近似的精细积分方法，具体过程如下：

将基本区段τ进行精细划分，记精细区段为σ＝τ/2^P，可采Q项帕德级数近似举证指数，即：

其中

同样，为了避免舍入误差，在使用精细积分算法时，也只计算和存储增量部分。因此，将N_Q(A₀σ)和D_Q(A₀σ)表示为如下形式：

其中

则

其中

将公式56代入公式55，得到T1的表达式，并将T1代入公式35中，得到公式40。

本实施例所述的一种基于刀齿切削时程精细积分的高速加工稳定域预测方法，铣削系统的最典型的两种模型为单自由度动力学模型和双自由度动力学模型，具体如下：

单自由度动力学模型：具有X方向的单自由度铣削加工动力学模型可由下列方程表示

其中，m_t、ζ、ω_n、和和X(t)分别是刀具的质量、阻尼比、固有频率和坐标；T为时滞量，T＝60/NΩ且N为铣刀齿数，Ω为铣刀主轴转速，单位为rpm；h(t)表示周期为T的切削力系数。

令X′(t)＝[X(t)通过变换，式(57)可改写为

其中

双自由度动力学模型：在x和y方向上具有相等的模态参数的双自由度铣削动力学模型可由下列方程表示

其中q₁(t)和q₂(t)分别表示在两自由度系统中x，y方向上刀具的坐标

其中K_t和K_n分别是切向和法向的切削力系数；φ_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是径向浸入比。

令通过矩阵变换，则式(59)可改写为

其中

本实施例中铣刀刀齿数目为2，径向力系数和法向力系数分别为7.96×10⁸和1.68×10⁸，铣刀一阶固有频率为922×2×π，相对阻尼为0.011，质量为0.14899，逆铣。将时滞周期T等分为40个小间隔，将主轴转速与径向切削深度构成地平面划分为400×200网格。

将上述方法和参数通过Matlab软件进行编程，画出单自由度和两自由度的稳定性叶瓣图及收敛图，通过稳定性叶瓣图来预测铣削加工的稳定域，选取径向浸入比分别为0.1和1得到的稳定性叶瓣图如图1、2、3和4所示及选取轴向切深w分别为0.2mm和0.7mm得到的收敛图如图5和6所示。

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

Claims (8)

1.基于刀齿切削时程精细积分的高速加工稳定域预测方法，其特征在于，包括如下步骤：

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

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

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

步骤4：利用精细积分算法对传递矩阵方程进一步求解，得到求解后的传递矩阵方程；

步骤5：根据求解后的传递矩阵方程，得到相邻周期刀齿响应之间的映射关系；并构建相邻周期刀齿响应的传递矩阵Φ：

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

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

2.根据权利要求1所述的方法，其特征在于，所述得到时滞微分方程的空间状态方程的方法如下所示：

所述考虑再生颤振铣削动力学的时滞微分方程如公式2所示：

其中，M、C、q(t)、w和K_c(t)分别表示刀具的模态质量、模态阻尼、模态刚度矩阵、模态坐标、轴向切削深度和周期系数矩阵；T为时滞量，T＝60/NΩ，N为铣刀齿数，Ω为主轴转速，为q(t)的一阶导数；为q(t)的二阶导数；

令令x(t)＝[q(t) p(t)]^T，并通过空间变换，将公式2变换成如下形式：

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

3.根据权利要求2所述的方法，其特征在于，所述得到状态项、时滞项二者之间的关系方程的具体方法为：

将周期T等分成m个小间隔，则具有m+1个节点，步长τ＝T/m；在任意一个时间区间[t_K，t_K+1]内对公式2进行积分计算，得到公式4，其中t_K表示第K个节点的时刻；K＝1，2，3，...，m；

其中，x_K+1表示时刻t_K+1的响应；x_K表示时刻t_K的响应；δ为时间变量，0≤δ≤τ或者t_K≤δ≤t_K+1；

令B(δ)x(δ)＝f(δ)为状态项，B(δ)x(δ-T)＝g(δ-T)为时滞项，得到状态项与时滞项二者之间的关系方程，如公式5：

4.根据权利要求3所述的方法，其特征在于，所述得到传递矩阵方程的具体方法为：

步骤3.1：通过构建三次牛顿插值多项式逼近状态项f(δ)，得到f(δ)的表达式，具体方法如下所示：

在时间区域[t_K-2，t_K+1]中，利用牛顿插值多项式算法将该时间区域中的四个时刻t_K-2、t_K-1、t_K、t_K+1以及这四个时刻的响应f_K-2、f_K-1、f_K、f_K+1对状态项进行插值计算，得到状态项f(δ)的方程表达式为：

f(δ)＝a₁f_K+1+b₁f_K+c₁f_K-1+d₁f_K-2 (6)

其中，a₁、b₁、b₁、d₁如下所示；

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

g(δ-T)＝a₂g_K-m+b₂g_K+1-m+c₂g_K+2-m (7)

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

步骤3.2：将公式6和公式7代入公式5中，得到公式8，即传递矩阵方程：

其中

5.根据权利要求4所述的方法，其特征在于，所述利用精细积分算法对传递矩阵方程进一步求解具体为：利用精细积分对进行求解，具体的方法为：

步骤4.1：对任意一个步长为τ的间隔进行精细划分，划分为2^p段精细区段，则任意一段精细则：

对T1进行泰勒级数展开：

其中Ta为增量，I为单位矩阵；

将公式10，代入公式9，得到公式11：

经过P次迭代计算，得到：

把公式12代入到公式8，中得到：

H_K，-2X_K-2+H_K，-1X_K-1+(H_K，0+T1)x_K+(H_K，1-I)＝H_K，mX_K-m+H_K，m-1X_K+1-m+H_K，m-2X_K+2-m (13)

其中H_K，1＝(h_K，1+i_K，1)B_K+1；

H_K，0＝(h_K，0+i_K，0)B_K；

H_K，-1＝(h_K，-1+i_K，-1)B_K-1；

H_K，-2＝(h_K，-2+i_K，-2)B_K-2；

H_K，m＝(h_K，m+i_K，m)B_K；

H_K，m-1＝(h_K，m-1+i_K，m-1)B_K+1；

H_K，m-2＝(h_K，m-2+i_K，m-2)B_K+2；

其中B_K+1、B_K、B_K-1、B_K-2、B_K+2为节点K+1、K、K-1、K-2、K+2处的周期系数。

6.根据权利要求5所述的方法，其特征在于，所述构建相邻周期刀齿响应的传递矩阵Φ的具体方法为：

对求解后的传递矩阵方程进行若干次迭代计算，得到相邻相邻周期刀齿响应之间的映射关系如下所示：

其中

根据映射关系构建相邻周期刀齿响应传递矩阵Φ：

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

7.根据权利要求5所述的方法，其特征在于，所述利用精细积分算法对传递矩阵方程进一步求解具体为利用精细积分对进行求解，具体的方法为：

对任意一个步长为τ的间隔进行精细划分，得到2^p段精细区段，则任意一段精细 利用Q项帕德级数近似举证指数得到：

其中：

在使用精细积分时，只计算和存储增量部分，故将N_Q(A₀σ)和D_Q(A₀σ)表示为如下形式：

其中

则

则

将公式24代入公式23，得到T1的表达式，并将T1代入公式8中，得到公式13。

8.根据权利要求4所述的方法，其特征在于，上述用埃米尔特插值多项式算法逼近时滞项，得到时滞项g(δ-T)的方程表达式的具体方法为：

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

其中，为响应g_K-m的一阶导数，为响应g_K+1-m的一阶导数；系数μ₀(δ)、μ₁(δ)、如下所示：

当δ＝t_K-m时，

当δ＝t_K+1-m时，

步骤S2：由于令τ＝t_K+1-m-t_K-m；并将上述公式27、公式26和μ₀(δ)、μ₁(δ)、代入公式25中；得到g(δ-T)的表达式如公式28所示：