CN105843177A - 铣削加工主轴转速正弦调制参数优化方法 - Google Patents

Abstract

本发明提供了一种铣削加工主轴转速正弦调制参数优化方法，包括：对变转速铣削加工系统进行动力学建模，建立变时滞二阶微分动力学方程；建立主轴转速正弦调制参数与时滞变量之间的映射关系并求取时滞；对动力学方程进行状态空间变换，得到变换后的状态空间方程；在相邻两个Floquet周期内对状态空间方程进行变步长离散；利用变步长数值积分法判定加工系统的稳定性；以获得最大加工效率为目标，以无颤振加工和主轴变速极限为约束，建立约束优化模型；利用智能优化算法获得优化后的变转速铣削正弦调制参数。本发明采用优化后的正弦调制参数进行变转速加工可以极大地提高加工效率。

Description

铣削加工主轴转速正弦调制参数优化方法

技术领域

本发明涉及机械加工技术领域，具体地，涉及一种铣削加工主轴转速正弦调制参数优化方法。

背景技术

铣削是最为常见的机械加工方式之一，常用于叶轮叶片等复杂曲面类零件的切削加工。铣削过程中如果加工参数选择不当则极易发生再生型颤振，再生型颤振属于自激振动，它严重影响加工质量并对刀具和机床主轴产生不同程度的危害。

在铣削加工过程中采用主轴转速连续变化(如正弦变化)的策略可以有效地打破再生颤振的发生机制，提高无颤振稳定铣削的加工效率。然而主轴转速调制是一把双刃剑，不恰当的调制参数不但不会提高稳定铣削参数极限，反而可能使加工振动更为剧烈。为了更好的发挥变转速铣削的颤振抑制作用，首先要对变转速铣削加工系统进行动力学建模，然后采用合理的优化方法对转速调制参数进行优化。目前，针对变转速铣削转速调制参数优化的方法较少，且均是通过扫描调制参数空间绘制高维参数化图谱然后选择最优调制参数的方法进行变速参数优化，存在计算精度低、优化时间长等缺点，因此提出高精高效的铣削转速调制参数智能优化方法，对于避免加工颤振、提高加工质量具有十分重要的意义和背景。

发明内容

针对现有技术中的缺陷，本发明的目的是提供一种铣削加工主轴转速正弦调制参数优化方法。

根据本发明提供的铣削加工主轴转速正弦调制参数优化方法，包括如下步骤：

步骤1：对变转速铣削加工系统进行动力学建模，建立变时滞二阶微分动力学方程；

步骤2：建立主轴转速正弦调制参数与时滞变量之间的映射关系并求取时滞；

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

步骤4：在动力学方程的相邻两个Floquet周期内对状态空间方程进行变步长离散；

步骤5：利用变步长数值积分法判定所述变转速铣削加工系统的稳定性；

步骤6：以获得最大加工效率为目标，以无颤振加工和机床主轴变速极限为约束，建立转速调制参数约束优化模型；

步骤7：利用智能优化算法对正弦调制参数进行优化，获得优化后的调速参数。

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

对变转速铣削加工系统进行动力学建模后得到的动力学方程为：

其中，

其中，M为模态质量矩阵，C为模态阻尼矩阵，K为模态刚度矩阵，为的加速度状态向量，为第t个时刻的速度状态向量，q(t)为第t个时刻的位移状态向量，a_p为轴向切深，K_c(t)为第t个时刻的切削系数矩阵，t为时刻，τ(t)为第t个时刻的时滞变量，f₀(t)为t时刻与动态切厚无关的切削力分量f₀的值，h_xx(t)、h_xy(t)、h_yx(t)和h_yy(t)均为切削力系数函数,下标xx、xy、yx、yy分别表示切削力系数函数在系数矩阵中的(1,1)、(1,2)、(2,1)、(2,2)位置)，j为刀齿变量，N为刀齿数目，g(φ_j(t))表示开关函数，φ_j(t)表示圆弧角，K_tc、K_te、K_nc、K_ne分别表示切向切削力系数、切向刃口力系数、径向切削力系数、径向刃口力系数，f₁₁(t)和f₂₁(t)分别为t时刻与动态切厚无关的切削力分量函数，下标11和21分别表示分量函数位于切削力分量矩阵的(1,1)、(2,1)位置，f_z(t)为t时刻每齿进给量，下标z表示进给量按每齿计算。

优选地，在所述步骤2，时滞变量τ(t)与主轴转速正弦调制参数RVA和RVF的映射关系如下：

其中，τ′(t)为第t个时刻的时滞变量的导数，Ω₀为名义主轴转速，RVA为正弦调幅参数，T为正弦调制周期，t为时间，τ(t)为t时刻的时滞变量；正弦调制周期T由正弦调频参数RVF和名义主轴转速Ω₀根据公式求得。

利用经典四阶龙格-库塔法求取时滞变量如下：

其中，z(t,τ(t))＝τ′(t)，τ_n+1、τ_n分别为在t_n+1时刻、t_n时刻的时滞变量值，h为离散步长。

优选地，所述步骤3，具体步骤为：对所述动力学方程进行状态空间变换，获得其状态空间方程：

其中，为状态向量，x(t)为t时刻x方向的振动位移，y(t)为t时刻y方向的振动位移，为为t时刻x方向的振动速度，为t时刻y方向的振动速度， 为状态向量的导数，由于a_pF₀(t)与动态切厚x(t)-x(t-τ(t))无关，因此不影响变转速铣削的稳定性。

优选地，所述步骤4，包括如下步骤：

步骤4.1，假定存在互质的整数p和q，使得如下公式成立：

pT＝qτ₀

其中，τ₀为名义时滞且则动力学方程的Floquet周期为pT；N为刀齿数目，Ω₀为名义主轴转速；

步骤4.2，根据如下两个公式分别求取一个Floquet周期内每个刀齿的切入时刻和切出时刻

根据以上两个公式求取的所有时刻点按升序排列成切入切出时间点集合φ_st为切入角，φ_ex为切出角，为刀齿j的切入时刻对应的刀具旋转角，为刀齿j的切出时刻对应的刀具旋转角；

步骤4.3，对于集合内的任意两个相邻元素构成的时间段如果刀齿处于切削状态即则根据离散误差要求对此时间段进行进一步离散，进入步骤4.4，否则不离散；其中，h_xx为切削力系数函数矩阵的(1,1)位置处元素，下标xx表示h_xx位于切削力系数函数矩阵的(1,1)位置，为集合中的第k个元素，为集合中的第(k+1)个元素；

步骤4.4，增加离散点和其中，δ是一个设定的相对于远小于离散步长的参数；

至此，周期pT已经被离散为许多不等距的小片段。

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

步骤5.1：对于任意离散区间[t_k,t_k+1]，k＝0,1,…,m-1，t_k为第k个离散时刻，t_k+1为第(k+1)个离散时刻，m为离散时刻点的个数，所述状态空间方程的解析解为：

由梯形公式得：

其中，B_k、B_k+1分别为B(t_k)、B(t_k+1)的简写；x_k为第k个离散时刻对应的状态向量，x_k+1为第(k+1)个离散时刻对应的状态向量，ξ为时刻变量，B(ξ)为ξ时刻B矩阵的值，x(ξ)为ξ时刻x向量的值，x(ξ-τ(ξ))为(ξ-τ(ξ))时刻x向量的值，x(t_k-τ_k)为(t_k-τ_k)时刻x向量的值，

步骤5.2：为了利用相邻两个Floquet周期[-pT,pT]上的离散点对上述公式中的时滞项x(t_k-τ_k)和x(t_k+1-τ_k+1)进行插值表示，假设

其中，p_k,q_k∈{1-m,2-m,…,0,…,m-1}；τ_k为第k个离散时刻对应的时滞，为在离散时刻序列中下标为p_k的时刻，为在离散时刻序列中下标为q_k的时刻；

步骤5.3：应用三点拉格朗日插值公式对时滞项进行插值表示，如下：

其中，

整理得，

其中， 其中，I为单位矩阵

步骤5.4：定义s＝max{k-q_k,k-p_k},(k＝0,1,...,m-1)，则获得如下离散映射关系：

y_k+1＝G_ky_k

其中，y_k是一个维数为(2s+4)×1的向量：

其中，x_k为k时刻x方向振动位移，y_k为k时刻y方向振动位移，为k时刻x方向振动速度，为k时刻y方向振动速度，col表示列向量；

G_k是一个(2s+4)×(2s+4)的系数矩阵：

其中，代表矩阵E_i第j行第k列的元素，E₀＝F_k+1 ^-1F_k， 如果不同的元素出现在矩阵G_k的相同位置，则将这些元素的代数和作为G_k此位置处的元素；

相邻两个Floquet周期[-pT,0]和[0,pT]之间的状态转移矩阵Φ为：

Φ＝G_m-1G_m-2…G₀

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

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

以获得最大加工效率为目标，以无颤振加工和机床主轴变速极限为约束，建立正弦调制参数约束优化模型。

max a_p

s.t.max{|ei(Φ)|}≤1

a_max≤a_lim

RVA^L≤RVA≤RVA^U

RVF^L≤RVF≤RVF^U

其中，max{|ei(Φ)|}表示状态转移矩阵Φ的模长最大的特征值的模，a_lim表示机床的主轴的变速性能(单位转/秒²)，RVA^L和RVA^U分别表示幅值调制参数RVA的下界和上界，RVF^L和RVF^U分别表示频率调制参数RVF的下界和上界，a_max表示转速调制的最大加速度，对于正弦调速可得：

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

利用罚函数理论处理约束优化模型中的不等式约束，建立新的约束优化模型如下：

min-a_p+σ₁·(max{0,max{|ei(Φ)|}-1})²+σ₂·(max{0,a_max-a_lim})²

s.t.RVA^L≤RVA≤RVA^U

RVF^L≤RVF≤RVF^U

其中，σ₁和σ₂为两个惩罚因子；再利用智能优化算法对正弦调制参数进行优化，获得优化后的变速参数。

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

1、本发明提出的铣削加工主轴转速正弦调制参数优化方法，首次采用了智能优化算法对调制参数进行优化，并在变转速铣削加工系统的稳定性约束环节采用了比现有算法具有更高计算效率精度的变步长数值积分法，因此极大地缩短了转速正弦调制参数优化时间、最大程度保证了优化结果的精度和可靠性；

2、本发明与恒转速铣削加工相比，采用优化后的正弦调制参数进行变转速铣削加工可以显著提高无颤振加工效率。

附图说明

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

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

图2为变转速铣削加工系统在“正弦调频-正弦调幅-轴向切深”参数空间中的三维稳定性图谱。

具体实施方式

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

请同时参阅图1至图2。

对变转速铣削加工系统进行动力学建模，建立变时滞二阶微分动力学方程；

建立主轴转速正弦调制参数与时滞变量之间的映射关系并求取时滞；

对所述动力学方程进行状态空间变换，得到状态空间方程；

在相邻两个Floquet周期内对状态空间方程进行变步长离散；

利用变步长数值积分法判定加工系统的稳定性；

以获得最大加工效率为目标，以无颤振加工和机床主轴变速极限为约束，建立转速调制参数约束优化模型；

利用智能优化算法对正弦调制参数进行优化，获得优化后的调速参数。

具体地，本实施例提供了一种铣削加工主轴转速正弦调制参数优化方法，包括如下步骤：

步骤1：对变转速铣削加工系统进行动力学建模，建立变时滞二阶微分动力学方程；

步骤2：建立主轴转速正弦调制参数与时滞变量之间的映射关系并求取时滞；

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

步骤4：在相邻两个Floquet周期内对状态空间方程进行变步长离散；

步骤5：利用变步长数值积分法判定加工系统的稳定性；

步骤6：以获得最大加工效率为目标，以无颤振加工和机床主轴变速极限为约束，建立转速调制参数约束优化模型；

步骤7：利用智能优化算法对正弦调制参数进行优化，获得优化后的调速参数。

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

两自由度铣削加工系统的动力学方程可以表述为：

其中，M为模态质量矩阵，C为模态阻尼矩阵，K为模态刚度矩阵，为加速度状态向量，为速度状态向量，q(t)为位移状态向量，F(t)为切削力矩阵；q(t)＝[x(t),y(t)]^T，x(t)和y(t)分别为x和y方向的振动位移；F(t)＝[F_x,F_y]，F_x和F_y分别为x和y方向的切削分力；

作用在刀齿j上的切向力F_t,j和径向力F_r,j分别表述为：

其中，a_p为轴向切深，K_tc、K_te、K_nc和K_ne分别表示切向切削力系数、切向刃口力系数、径向切削力系数和径向刃口力系数，h_j(t)为未变形切削厚度，可由以下公式求得：

其中，f_z(t)为每齿进给量，τ(t)为时滞变量，φ_j(t)表示刀齿j在t时刻的角位移：

其中，Ω(s)表示角位移变量，N为刀齿齿数；

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

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

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

将公式(5)代入公式(1)，得到变转速加工系统的动力学方程如下：

其中，

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

正弦调制后的主轴转速Ω(t)随时间t变化，可由以下公式表述：

Ω(t)＝Ω₀+Ω₀RVAS(t) (8)

其中，Ω₀为名义主轴转速，RVA为正弦调幅参数，S(t)为调制规律，可由以下公式求得：

其中，T为正弦调制周期；

正弦调幅参数RVA和正弦调频参数RVF可分别由公式(10)和公式(11)求得：

其中，Ω_A为调制幅值，f为调制频率，f＝1/T；

变转速铣削动力学方程(7)中的时滞τ(t)是时变的，其满足如下公式：

根据公式(12)很难获得τ(t)的解析表达式，因此使用高效高精数值方法对其进行求解；

公式(12)两边对时间t求导，可得：

定义τ′(t)＝z(t,τ(t))，然后利用经典四阶龙格-库塔方法求解时滞τ(t)：

其中，τ_n+1和τ_n分别为在t_n+1时刻和t_n时刻的时滞变量值，h为离散步长；公式(14)的初始值(t₀,τ(t₀))可根据公式(12)由数值方法求得。

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

对动力学方程(7)进行状态空间变换，获得其状态空间表达式：

其中，为状态向量，为状态向量的导数，由于a_pF₀(t)与动态切厚x(t)-x(t-τ(t))无关，因此不影响变转速铣削的稳定性。

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

步骤4.1，假定存在互质的整数p和q，使得如下公式成立：

pT＝qτ₀ (16)

其中，τ₀为名义时滞且则动力学方程的Floquet周期为pT；

步骤4.2，根据如下两个公式分别求取一个周期内每个刀齿的切入时刻和切出时刻

根据公式(17)和公式(18)求取的所有时刻点按升序排列成切入切出时间点集合

步骤4.3，对于集合内的任意两个相邻元素构成的时间段如果刀齿处于切削状态即则根据离散误差要求对此时间段进行进一步离散，否则不离散；

步骤4.4，增加离散点和其中，δ是一个相对于离散步长很小的数；

至此，周期pT已经被离散为许多不等距的小片段。

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

对于任意离散区间[t_k,t_k+1](k＝0,1,…,m-1)，加工系统动力学方程状态空间表达式的解析解为：

由梯形公式可得：

其中，B_k和B_k+1分别为B(t_k)和B(t_k+1)的简写；

为了利用相邻两个Floquet周期[-pT,pT]上的离散点对上述公式中的时滞项x(t_k-τ_k)和x(t_k+1-τ_k+1)进行插值表示，我们假设

其中，p_k,q_k∈{1-m,2-m,…,0,…,m-1}；

应用三点拉格朗日插值公式对时滞项进行插值表示，如下：

其中，

将公式(23)和公式(24)代入公式(20)可得：

其中，

定义s＝max{k-q_k,k-p_k},(k＝0,1,...,m-1)，则可获得如下离散映射关系：

y_k+1＝G_ky_k (26)

其中，y_k是一个维数为(2s+4)×1的向量：

G_k是一个(2s+4)×(2s+4)的系数矩阵：

其中，代表矩阵E_i第j行第k列的矩阵，E₀＝F_k+1 ^-1F_k， 如果不同的元素出现在矩阵G_k的相同位置，则将这些元素的代数和作为G_k此位置处的元素；

相邻两个Floquet周期[-pT,0]和[0,pT]之间的状态转移矩阵为：

Φ＝G_m-1G_m-2…G₀ (29)

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

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

以获得最大加工效率为目标，以无颤振加工和机床主轴变速极限为约束，建立正弦调制参数约束优化模型。

max a_p

s.t.max{|ei(Φ)|}≤1

a_max≤a_lim

RVA^L≤RVA≤RVA^U

RVF^L≤RVF≤RVF^U

其中，max{|ei(Φ)|}表示状态转移矩阵Φ的模长最大的特征值的模，a_lim表示机床的主轴的变速性能(单位转/秒²)，RVA^L和RVA^U分别表示幅值调制参数RVA的下界和上界，RVF^L和RVF^U分别表示频率调制参数RVF的下界和上界，a_max表示转速调制的最大加速度，对于正弦调速可得：

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

利用罚函数理论处理约束优化模型中的不等式约束，建立新的约束优化模型如下：

min-a_p+σ₁·(max{0,max{|ei(Φ)|}-1})²+σ₂·(max{0,a_max-a_lim})²

s.t.RVA^L≤RVA≤RVA^U

RVF^L≤RVF≤RVF^U

其中，σ₁和σ₂为两个惩罚因子。再利用智能优化算法对正弦调制参数进行优化，获得优化后的变速参数。

下面结合具体加工实例说明本发明的具体实施方案，实例参数引自文献1Bayly,P.V.,Mann,B.P.,Schmitz,T.L.,Peters,D.A.,Stepan,G.,and Insperger,T.,"Effectsof radial immersion and cutting direction on chatter instability in end-milling,"Proc.ASME International Mechanical Engineering Congress andExposition,Proceedings,pp.351-363。铣刀直径D＝1.27×10^-2m，齿数N＝2，径向切深比a_r/D＝0.1，固有频率f_n＝922Hz，阻尼比刚度k＝1.34×10⁶N/m，切削力系数K_tc＝6×10⁸N/m²、K_rc＝2×10⁸N/m²，机床主轴的加速度极限为a_lim＝700rev/s²，选定的名义转速为Ω₀＝6000rpm。

将已知参数代入发明内容中的步骤1-步骤7。

利用遗传算法获得的优化结果为：幅值调制参数RVA＝0.168，频率调制参数RVF＝0.043，无颤振稳定加工轴向切深极限为3.8mm。并以变步长数值积分法获得的三维“正弦调频-正弦调幅-轴向切深”参数空间中的稳定性图谱作为验证，如图2所示。结果显示智能优化算法优化的结果与高维稳定性图谱获得的结果一致。与相同条件下恒转速铣削无颤振稳定加工轴向切深极限0.8mm相比，采用本发明优化后的正弦调制参数可以将加工效率提高375％。

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

Claims (8)

1.一种铣削加工主轴转速正弦调制参数优化方法，其特征在于，包括如下步骤：

步骤1：对变转速铣削加工系统进行动力学建模，建立变时滞二阶微分动力学方程；

步骤2：建立主轴转速正弦调制参数与时滞变量之间的映射关系并求取时滞；

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

步骤4：在动力学方程的相邻两个Floquet周期内对状态空间方程进行变步长离散；

步骤5：利用变步长数值积分法判定所述变转速铣削加工系统的稳定性；

步骤6：以获得最大加工效率为目标，以无颤振加工和机床主轴变速极限为约束，建立转速调制参数约束优化模型；

步骤7：利用智能优化算法对正弦调制参数进行优化，获得优化后的调速参数。

2.根据权利要求1所述的铣削加工主轴转速正弦调制参数优化方法，其特征在于，所述步骤1，具体为：

对变转速铣削加工系统进行动力学建模后得到的动力学方程为：

其中，

其中，M为模态质量矩阵，C为模态阻尼矩阵，K为模态刚度矩阵，为的加速度状态向量，为第t个时刻的速度状态向量，q(t)为第t个时刻的位移状态向量，a_p为轴向切深，K_c(t)为第t个时刻的切削系数矩阵，t为时刻，τ(t)为第t个时刻的时滞变量，f₀(t)为t时刻与动态切厚无关的切削力分量f₀的值，h_xx(t)、h_xy(t)、h_yx(t)和h_yy(t)均 为切削力系数函数,下标xx、xy、yx、yy分别表示切削力系数函数在系数矩阵中的(1,1)、(1,2)、(2,1)、(2,2)位置)，j为刀齿变量，N为刀齿数目，g(φ_j(t))表示开关函数，φ_j(t)表示圆弧角，K_tc、K_te、K_nc、K_ne分别表示切向切削力系数、切向刃口力系数、径向切削力系数、径向刃口力系数，f₁₁(t)和f₂₁(t)分别为t时刻与动态切厚无关的切削力分量函数，下标11和21分别表示分量函数位于切削力分量矩阵的(1,1)、(2,1)位置，f_z(t)为t时刻每齿进给量，下标z表示进给量按每齿计算。

3.根据权利要求1所述的铣削加工主轴转速正弦调制参数优化方法，其特征在于，在所述步骤2，时滞变量τ(t)与主轴转速正弦调制参数RVA和RVF的映射关系如下：

其中，τ′(t)为第t个时刻的时滞变量的导数，Ω₀为名义主轴转速，RVA为正弦调幅参数，T为正弦调制周期，t为时间，τ(t)为t时刻的时滞变量；正弦调制周期T由正弦调频参数RVF和名义主轴转速Ω₀根据公式求得；

利用经典四阶龙格-库塔法求取时滞变量如下：

其中，z(t,τ(t))＝τ′(t)，τ_n+1、τ_n分别为在t_n+1时刻、t_n时刻的时滞变量值，h为离散步长。

4.根据权利要求2所述的铣削加工主轴转速正弦调制参数优化方法，其特征在于，所述步骤3，具体步骤为：对所述动力学方程进行状态空间变换，获得其状态空间方程：

其中，为状态向量，x(t)为t时刻x方向的振动位移，y(t)为t时刻y方向的振动位移，为为t时刻x方向的振动速度，为t时刻y方向的振动速度， 为状态向量的导数，由于a_pF₀(t)与动态切厚x(t)-x(t-τ(t))无关，因此不影响变转速铣削的稳定性。

5.根据权利要求1所述的铣削加工主轴转速正弦调制参数优化方法，其特征在于，所述步骤4，包括如下步骤：

步骤4.1，假定存在互质的整数p和q，使得如下公式成立：

pT＝qτ₀

其中，τ₀为名义时滞且则动力学方程的Floquet周期为pT；N为刀齿数目，Ω₀为名义主轴转速；

步骤4.2，根据如下两个公式分别求取一个Floquet周期内每个刀齿的切入时刻和切出时刻

根据以上两个公式求取的所有时刻点按升序排列成切入切出时间点集合φ_st为切入角，φ_ex为切出角，为刀齿j的切入时刻对应的刀具旋转角，为刀齿j的切出时刻对应的刀具旋转角；

步骤4.3，对于集合内的任意两个相邻元素构成的时间段如果刀齿处于切削状态即则根据离散误差要求对此时间段进行进一步离散，进入步骤4.4，否则不离散；其中，h_xx为切削力系数函数矩阵的(1,1)位置处元素，下标xx表示h_xx位于切削力系数函数矩阵的(1,1)位置，为集合中的第k个元素，为集合中的第(k+1)个元素；

步骤4.4，增加离散点和其中，δ是一个设定的相对于远小于离散步长的参数；

至此，周期pT已经被离散为许多不等距的小片段。

6.根据权利要求5所述的铣削加工主轴转速正弦调制参数优化方法，其特征在于，所述步骤5，具体为：

步骤5.1：对于任意离散区间[t_k,t_k+1]，k＝0,1,…,m-1，t_k为第k个离散时刻，t_k+1为第(k+1)个离散时刻，m为离散时刻点的个数，所述状态空间方程的解析解为：

由梯形公式得：

其中，B_k、B_k+1分别为B(t_k)、B(t_k+1)的简写；x_k为第k个离散时刻对应的状态向量，x_k+1为第(k+1)个离散时刻对应的状态向量，ξ为时刻变量，B(ξ)为ξ时刻B矩阵的值，x(ξ)为ξ时刻x向量的值，x(ξ-τ(ξ))为(ξ-τ(ξ))时刻x向量的值，x(t_k-τ_k)为(t_k-τ_k)时刻x向量的值，

步骤5.2：为了利用相邻两个Floquet周期[-pT,pT]上的离散点对上述公式中的时滞项x(t_k-τ_k)和x(t_k+1-τ_k+1)进行插值表示，假设

其中，p_k,q_k∈{1-m,2-m,…,0…,m-1}；τ_k为第k个离散时刻对应的时滞，为在离散时刻序列中下标为p_k的时刻，为在离散时刻序列中下标为q_k的时刻；

步骤5.3：应用三点拉格朗日插值公式对时滞项进行插值表示，如下：

其中，

整理得，

其中， 其中，I为单位矩阵

步骤5.4：定义s＝max{k-q_k,k-p_k},(k＝0,1,…,m-1)，则获得如下离散映射关系：

y_k+1＝G_ky_k

其中，y_k是一个维数为(2s+4)×1的向量：

其中，x_k为k时刻x方向振动位移，y_k为k时刻y方向振动位移，为k时刻x方向振动速度，为k时刻y方向振动速度，col表示列向量；

G_k是一个(2s+4)×(2s+4)的系数矩阵：

其中，代表矩阵E_i第j行第k列的元素，E₀＝F_k+1 ^-1F_k， 如果不同的元素出现在矩阵G_k的相同位置，则将这些元素的代数和作为G_k此位置处的元素；

相邻两个Floquet周期[-pT,0]和[0,pT]之间的状态转移矩阵Φ为：

Φ＝G_m-1G_m-2…G₀

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

7.根据权利要求2所述的铣削加工主轴转速正弦调制参数优化方法，其特征在于，所述步骤6，具体为：

以获得最大加工效率为目标，以无颤振加工和机床主轴变速极限为约束，建立正弦调制参数约束优化模型。

其中，max{|ei(Φ)|}表示状态转移矩阵Φ的模长最大的特征值的模，a_lim表示机床的主轴的变速性能，单位转/秒²，RVA^L和RVA^U分别表示幅值调制参数RVA的下界和上界，RVF^L和RVF^U分别表示频率调制参数RVF的下界和上界，a_max表示转速调制的最大加速度，对于正弦调速可得：

8.根据权利要求7所述的铣削加工主轴转速正弦调制参数优化方法，其特征在于，所述步骤7，具体为：

利用罚函数理论处理约束优化模型中的不等式约束，建立新的约束优化模型如下：

其中，σ₁和σ₂为两个惩罚因子；再利用智能优化算法对正弦调制参数进行优化，获得优化后的变速参数。