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

摘要

本发明涉及先进制造领域，具体涉及一种基于三次多项式的铣削稳定性预测方法，本发明采用三次多项式来逼近动力学方程中的状态项、时滞项和周期系数项，采用多个已知时间点及其响应来拟合所需项，减小了计算方法的局部误差，从而提高了预测方法的精度；同时在获得稳定性叶瓣图的过程中，引入H矩阵，而不是直接代入F矩阵进行计算，减小了F矩阵计算过程中的迭代次数，从而节约计算方法的时间，提高计算效率。

主权项

一种基于三次多项式的铣削稳定性预测方法，其特征在于其包括以下步骤：①建立铣刀在单自由度铣削过程中的动力学方程：x·(t)=Ax(t)+B(t)x(t)-B(t)x(t-τ)---(1)]]>其中，为常系数矩阵，为随时间周期变化的系数矩阵，x(t)表示刀具在t时刻的状态响应，ωn表示刀尖点的固有频率，ζ表示相对阻尼，mt表示模态质量，w表示轴向切削深度，τ表示时滞；h(t)表示瞬时切屑厚度，其表达式为：h(t)=Σj=1Ng(φj(t))sin((φj(t))[Ktcos(φj(t))+Knsin(φj(t))]---(2)]]>式(2)中，N表示铣刀的刀齿数目，Kt和Kn分别为切向和法向的切削力系数，φj(t)为第j刀齿的角位移，表达式为φj(t)＝(2πΩ/60)t+(j‑1)·2π/N，窗函数g(φj(t))定义式为：式(3)中，φst和φex分别为第j刀齿的切入和切出角，当采用顺铣时，φst＝arccos(2ae/D‑1)，φex＝π；当采用逆铣时，φst＝0，φex＝arccos(1‑2ae/D)，ae/D为径向浸入比，即径向切深/刀具直径的比值；②将单自由度的铣削过程动力学方程(1)的时滞项τ平均分为m个小区间，则时间步长为其中任意一个时间小区间表示为[ti，ti+1]，i＝1，2，3，…m，将方程(1)在时间小区间[ti，ti+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)，具体过程如下：在构建三次勒让德多项式的过程中，需要用到前四项勒让德多项式：P0(s)＝1  (5)P1(s)＝S  (6)P2(s)=12(3s2-1)---(7)]]>P3(s)=12(5s3-3s)---(8)]]>在区间[ti‑2，ti+1]上构建三次勒让德多项式，Pn(s)是区间[‑1，1]上关于权函数ω(x)≡1的正交函数系，在应用勒让德多项式的过程中，需要将区间[ti‑2，ti+1]变换到[‑1，1]上，再利用勒让德多项式进行拟合，令当t在区间[ti‑2，ti+1]上变化时，对应的S在[‑1，1]上变化，变换后的表达式为：P~n(s)=Pn(t)=Pn(2s-(ti+1+ti-2)ti+1-ti-2),(n=0,1,...,3)---(9)]]>变换后的个体表达式为：P0(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)]]>针对时域区间[ti‑2，ti+1]，利用该区间上的四个时间点ti‑2、ti‑1、ti、ti+1及其对应的响应xi‑2、xi‑1、xi、xi+1对状态项进行拟合，可以表示为：x(t)＝a0·P0(s)+a1·P1(s)+a2·P2(s)+a3·P3（s)   (14)各项的系数表达式为其中ωj为权值函数，ωj≡1，令ti‑1＝‑Δt，ti＝0，ti+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)]]>再分别将ti‑1＝‑Δt，ti＝0，ti+1＝Δt、(16)、(17)、(18)式代入(15)式得a0=xi-2+xi-1+xi+xi-14---(19)]]>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)分别用一次勒让德多项式拟合的多项式来表示，拟合过程中将时间区间[ti，ti+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)]]>其中F0＝eAΔ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)]]>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)]]>其中Pi＝[I‑H20Bi‑H21Bi+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)]]>系统的离散映射可以表示为Ψ＝MmMm‑1…M1，Ψ即为系统的Floquet转移矩阵，其中Mi=M11iM12iM13i...0M1kiM1,k+1iI00...0000I0...000.....................00000I0---(33)]]>⑤计算Floquet转移矩阵ψ的特征值，通过特征值的模判定系统的稳定性，具体的判定准则如下：