[发明专利]一种基于三次多项式的铣削稳定性预测方法

权利要求书

1.一种基于三次多项式的铣削稳定性预测方法，其特征在于其包括以下步骤：

①建立铣刀在单自由度铣削过程中的动力学方程：

x·(t)=Ax(t)+B(t)x(t)-B(t)x(t-τ)---(1)]]>

其中，为常系数矩阵，为随时间周期变化的系数矩阵，x(t)表示刀具在t时刻的状态响应，ω_n表示刀尖点的固有频率，ζ表示相对阻尼，m_t表示模态质量，w表示轴向切削深度，τ表示时滞；

h(t)表示瞬时切屑厚度，其表达式为：

h(t)=Σj=1Ng(φj(t))sin((φj(t))[Ktcos(φj(t))+Knsin(φj(t))]---(2)]]>

式(2)中，N表示铣刀的刀齿数目，K_t和K_n分别为切向和法向的切削力系数，φ_j(t)为第j刀齿的角位移，表达式为φ_j(t)＝(2πΩ/60)t+(j-1)·2π/N，窗函数g(φ_j(t))定义式为：

式(3)中，φ_st和φ_ex分别为第j刀齿的切入和切出角，当采用顺铣时，φ_st＝arccos(2a_e/D-1)，φ_ex＝π；当采用逆铣时，φ_st＝0，φ_ex＝arccos(1-2a_e/D)，a_e/D为径向浸入比，即径向切深/刀具直径的比值；

②将单自由度的铣削过程动力学方程(1)的时滞项τ平均分为m个小区间，则时间步长为其中任意一个时间小区间表示为[t_i，t_i+1]，i＝1，2，3，…m，将方程(1)在时间小区间[t_i，t_i+1]上进行积分，得到

xi+1=eAΔtxi+∫titi+1eA(ti+1-s)[B(s)x(s)-B(s)x(s-τ)]ds---(4)]]>

③通过构建三次勒让德多项式来拟合步骤②中式(4)的状态项X(s)、时滞状态项x(s-τ)和随时间变化的周期系数项B(s)，具体过程如下：

在构建三次勒让德多项式的过程中，需要用到前四项勒让德多项式：

P₀(s)＝1(5)

P₁(s)＝S(6)

P2(s)=12(3s2-1)---(7)]]>

P3(s)=12(5s3-3s)---(8)]]>

在区间[t_i-2，t_i+1]上构建三次勒让德多项式，P_n(s)是区间[-1，1]上关于权函数ω(x)≡1的正交函数系，在应用勒让德多项式的过程中，需要将区间[t_i-2，t_i+1]变换到[-1，1]上，再利用勒让德多项式进行拟合，

令当t在区间[t_i-2，t_i+1]上变化时，对应的S在[-1，1]上变化，变换后的表达式为：

P~n(s)=Pn(t)=Pn(2s-(ti+1+ti-2)ti+1-ti-2),(n=0,1,...,3)---(9)]]>

变换后的个体表达式为：

P₀(s)＝1(10)

P1(s)=2s-(ti+1+ti-2)ti+1-ti-2---(11)]]>

P2(s)=12[3(2s-(ti+1+ti-2)ti+1-ti-2)2-1]---(12)]]>

P3(s)=12[5·(2s-(ti+1+ti-2)ti+1-ti-2)3-3·2s-(ti+1+ti-2)ti+1-ti-2]---(13)]]>

针对时域区间[t_i-2，t_i+1]，利用该区间上的四个时间点t_i-2、t_i-1、t_i、t_i+1及其对应的响应x_i-2、x_i-1、x_i、x_i+1对状态项进行拟合，可以表示为：

x(t)＝a₀·P₀(s)+a₁·P₁(s)+a₂·P₂(s)+a₃·P₃（s) (14)

各项的系数表达式为

其中ω_j为权值函数，ω_j≡1，

令t_i-1＝-Δt，t_i＝0，t_i+1＝Δt，分别代入(11)、(12)、(13)式得

P1(s)=23Δt·s+13---(16)]]>

P2(s)=23(Δt)2·s2+23Δt·s-13---(17)]]>

P3(s)=2027(Δt)3·s3+109(Δt)2·s2-49Δt·s-1127---(18)]]>

再分别将t_i-1＝-Δt，t_i＝0，t_i+1＝Δt、(16)、(17)、(18)式代入(15)式得

a0=xi-2+xi-1+xi+xi-14---(19)]]> 2

a1=920(-xi-2-13xi-1+13xi+xi+1)---(20)]]>

a2=920(xi-2-13xi-1-13xi+xi-1)---(21)]]>

a3=7291700(-xi-2+1127xi-1-1127xi+xi+1)---(22)]]>

将(16)-(22)式代入(14)式，整理得：

x(s)=(53425+81425Δt·s-317(Δt)2·s2-2785(Δt)3·s3)·xi-2+(76425-118425Δt·s+885(Δt)2·s2+1185(Δt)3·s3)·xi-1+(179425+33425Δt·s-517(Δt)2·s2-1185(Δt)3·s3)·xi+(32425+174425Δt·s+6685(Δt)2·s2+2785(Δt)3·s3)·xi+1---(23)]]>

将时滞项x(s-τ)和周期系数项B(s)分别用一次勒让德多项式拟合的多项式来表示，拟合过程中将时间区间[t_i，t_i+1]变换到区间[-1，1]上，具体形式如下：

x(s-τ)=(1-sΔt)·xi-m+sΔt·xi-m+1---(24)]]>

B(s)=(1-sΔt)·Bi+sΔt·Bi+1---(25)]]>

④构建Floquet转移矩阵，将(23)、(24)、(25)式代入(4)式得：

xi+1=F0xi+(H14Bi+H15Bi+1)xi-2+(H16Bi+H17Bi+1)xi-1+(H18Bi+H19Bi+1)xi+(H20Bi+H21Bi+1)xi+1-(H11Bi+H12Bi+1)xi-m-(H12Bi+H13Bi+1)xi+1-m---(26)]]>

其中

F₀＝e^AΔt(27a)

F1=∫titi+1eA(ti+1-s)ds=∫0ΔteA(Δt-s)ds=(F0-I)A-1---(27b)]]>

F2=∫titi+1eA(ti+1-s)(s-ti)ds=∫0ΔteA(Δt-s)sds=(F1-(Δt)I)A-1---(27c)]]> 3

F3=∫titi+1eA(ti+1-s)(s-ti)2ds=∫0ΔteA(Δt-s)s2ds=(2F2-(Δt)2I)A-1---(27b)]]>

F4=∫titi+1eA(ti+1-s)(s-ti)3ds=∫0ΔteA(Δt-s)s3ds=(3F3-(Δt)3I)A-1---(27e)]]>

F5=∫titi+1eA(ti+1-s)(s-ti)4ds=∫0ΔteA(Δt-s)s4ds=(4F4-(Δt)4I)A-1---(27f)]]>

H11=1(Δt)2[(Δt)2·F1-2Δt·F2+F3]---(28a)]]>

H12=1(Δt)2[Δt·F2-F3]---(28b)]]>

H13=1(Δt)2·F3---(28c)]]>

H14=53425F1+28425ΔtF2-156425(Δt)2F3-1285(Δt)3F4+2785(Δt)4F5---(28d)]]>

H15=53425ΔtF2+81425(Δt)2F3-317(Δt)3F4-2785(Δt)4F5---(28e)]]>

H16=76425F1-194425ΔtF2+158425(Δt)2F3+385(Δt)3F4-1185(Δt)4F5---(28f)]]>

H17=76425ΔtF2-118425(Δt)2F3+885(Δt)3F4+1185(Δt)4F5---(28g)]]>

H18=179425F1-146425ΔtF2-158425(Δt)2F3+1485(Δt)3F4+1185(Δt)4F5---(28h)]]>

H19=179425ΔtF2+33425(Δt)2F3-517(Δt)3F4-1185(Δt)4F5---(28i)]]>

H20=32425F1+142425ΔtF2+156425(Δt)2F3-3985(Δt)3F4-2785(Δt)4F5---(28j)]]>

H21=32425ΔtF2+174425(Δt)2F3+6685(Δt)3F4+2785(Δt)4F5---(28k)]]>

方程(26)可写为

xi+1=Pi(H14Bi+H15Bi+1)xi-2+(H16Bi+H17Bi+1)xi-1+(F0+H18Bi+H19Bi+1)xi+(H20Bi+H21Bi+1)xi+1-(H11Bi+H12Bi+1)xi-m-(H12Bi+H13Bi+1)xi+1-m---(29)]]>

其中

P_i＝[I-H₂₀B_i-H₂₁B_i+1]^-1(30)

通过方程，可以得到各时间点振动位移映射关系，通过矩阵表示如下：

xi+1xixi-1...xi+1-m=M11iM12iM13i...0M1kiM1,k+1iI00...0000I0...000.....................0000II0xixi-1xi-2...xi-m---(31)]]>

其中

M11i=Pi(F0+H18Bi+H19Bi+1)---(32a)]]>

M12i=Pi(H16Bi+H17Bi+1)---(32b)]]>

M13i=Pi(H14Bi+H15Bi+1)---(32c)]]>

M1mi=-Pi(H12Bi+H13Bi+1)---(32d)]]>

Mi,m+1i=-Pi(H11Bi+H12Bi+1)---(32e)]]>

系统的离散映射可以表示为

Ψ＝M_mM_m-1…M₁，Ψ即为系统的Floquet转移矩阵，

其中

Mi=M11iM12iM13i...0M1kiM1,k+1iI00...0000I0...000.....................00000I0---(33)]]>

⑤计算Floquet转移矩阵ψ的特征值，通过特征值的模判定系统的稳定性，具体的判定准则如下：