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

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

Info

Publication number
CN109740264B
CN109740264B CN201910012032.6A CN201910012032A CN109740264B CN 109740264 B CN109740264 B CN 109740264B CN 201910012032 A CN201910012032 A CN 201910012032A CN 109740264 B CN109740264 B CN 109740264B
Authority
CN
China
Prior art keywords
term
time
equation
period
coefficient
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910012032.6A
Other languages
English (en)
Other versions
CN109740264A (zh
Inventor
杨文安
黄超
蔡旭林
刘伟超
郭宇
廖文和
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN201910012032.6A priority Critical patent/CN109740264B/zh
Publication of CN109740264A publication Critical patent/CN109740264A/zh
Application granted granted Critical
Publication of CN109740264B publication Critical patent/CN109740264B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Complex Calculations (AREA)
  • Numerical Control (AREA)

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所示：
Figure BDA0001937687870000021
其中λ()的含义为求特征值。
进一步的，在单自由度系统中所述步骤1的具体方法为：
步骤1.1：构建考虑再生颤振的单自由度铣削动力学的时滞微分方程，如公式2所示：
Figure BDA0001937687870000022
其中ξ、mt、ωn分别为刀具的阻尼比、质量和固有频率；h(t)表示周期为T的切削力系数；q(t)为刀具的坐标，
Figure BDA0001937687870000031
为q(t)的一阶导数，
Figure BDA0001937687870000032
为q(t)的二阶导数；w为轴向切深；
步骤1.2：令X(t)＝[q(t) p(t)]T，其中
Figure BDA0001937687870000033
将公式2转换成空间状态形式如公式3所示：
Figure BDA0001937687870000034
其中A0为常数矩阵
Figure BDA0001937687870000035
Figure BDA0001937687870000036
为X(t)的一阶导数；B(t)表示周期为T的系数矩阵，
Figure BDA0001937687870000037
进一步的，所述步骤2的具体方法为：将刀齿切削周期T等分成m个小间隔，则具有y＝m+1个节点，步长
Figure BDA0001937687870000038
将公式2在任意一个时间区间[tK，tK+1]内进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程，如公式4所示，其中tK表示第K个节点的时刻；所述K＝1,2,…,y；
Figure BDA0001937687870000039
其中X(δ)为状态项、X(δ-T)为时滞项、B(δ)为周期系数项，Xk+1表示时刻tK+1的响应；Xk表示时刻tK的响应；δ为时间变量，0≤δ≤τ或者tK≤δ≤tK+1。
进一步的，所述步骤3中确定状态项、时滞项、周期系数项的方程表达式的具体方法为：
在当前周期的时间区域[tK-2，tK+1]中，利用牛顿插值多项式算法将该时间区域中的四个时刻tK-2、tK-1、tK、tK+1以及这四个时刻的响应XK-2、XK-1、XK、XK+1逼近状态项，K>2，得到状态项X(δ)的方程表达式为：
X(δ)＝a1XK+1+b1XK+c1XK-1+d1XK-2 (5)
其中，a1、b1、b1、d1为系数，具体的表示如公式6所示；
Figure BDA0001937687870000041
在上一个周期的时间区域[tK，-m，tK+2，-m]中，利用分段三次埃米尔特插值多项式算法将三个时刻tK，-m、tK+1，-m、tK+2，-m以及这三个时刻的响应XK，-m、XK+1，-m、XK+2，-m逼近时滞项，得到时滞项X(δ-T)的方程表达式如公式7所示，所述tK，-m、tK+1，-m、tK+2，-m分别表示上一个周期中节点K、K+1、K+2的时刻：
X(δ-T)＝a2XK，-m+b2XK+1，-m+c2XK+2，-m (7)
其中a2、b2、c2为系数，具体表达式如公式8所示；
Figure BDA0001937687870000042
用线性插值算法逼近周期系数项，得到周期系数项B(δ)的方程表达式为：
Figure BDA0001937687870000043
其中BK+1表示在当前周期中第K+1个节点的周期系数；BK表示当前周期中第K个节点的周期系数。
进一步的，步骤3中的构建相邻周期刀齿响应的传递矩阵Ф具体方法如下：
步骤A：将公式5、7、9代入状态项、时滞项、周期系数项三者之间的关系方程中，整理得到公式10；
(HK，1—I)XK+1+(H0+Hk，0)XK+Hk，-1XK-1+Hk，-2XK-2
＝Hk，mXK，-m+Hk，m-1XK+1，-m+Hk，m-2XK+2，-m (10)
其中，
Figure BDA0001937687870000044
I为单位矩阵，Hk，1为当前周期中第K+1个节点响应的切削力系数；Hk，0为当前周期中第K个节点响应的切削力系数，Hk，-1为当前周期中第K-1个节点响应的切削力系数；Hk，-2为当前周期中第K-2个节点响应的切削力系数；Hk，m-2为上一个周期中第K-2个节点响应的切削力系数；Hk，m-1为上一个周期中第K-1个节点响应的切削力系数；Hk，m为上一个周期中第K个节点响应的切削力系数；
Figure BDA0001937687870000051
Figure BDA0001937687870000052
Figure BDA0001937687870000053
Figure BDA0001937687870000054
Figure BDA0001937687870000055
Figure BDA0001937687870000056
Figure BDA0001937687870000057
其中，
Figure BDA0001937687870000058
Figure BDA0001937687870000059
Figure BDA00019376878700000510
Figure BDA00019376878700000511
Figure BDA00019376878700000512
步骤B：将上述公式10进行m次的迭代计算，得到相邻周期刀齿响应之间的映射关系，用矩阵表示如公式11所示：
Figure BDA00019376878700000513
其中Xm为当前周期的中时刻tm的响应；Xm，-m为上一个周期中时刻tm，-m的响应；D1和D2表达式如下所示：
Figure BDA0001937687870000061
Figure BDA0001937687870000062
步骤C：构建相邻周期刀齿响应的传递矩阵Ф，如公式12所示：
Φ＝(D1)-1D2 (12)。
进一步的，上述用分段三次埃米尔特插值多项式算法逼近时滞项，得到时滞项X(δ-T)的方程表达式的具体方法为：
步骤S1:在上一个周期的每个时间区域[tK，-m，tK+1，-m]内，构造分段三次埃尔米特插值多项式，得到公式13：
Figure BDA0001937687870000063
其中，
Figure BDA0001937687870000064
为响应XK，-m的一阶导数，
Figure BDA0001937687870000065
为响应XK+1，-m的一阶导数；系数ε0(δ)、ε1(δ)、β0(δ)、β1具体的表述如公式14所示：
Figure BDA0001937687870000066
当δ＝tK，-m时，
Figure BDA0001937687870000071
当δ＝tK+1，-m时，
Figure BDA0001937687870000072
步骤S2：由于
Figure BDA0001937687870000073
令τ＝tK+1，-m—tK，-m；并将上述公式15和公式16代入公式13中；得到X(δ-T)的表达式如公式17所示：
Figure BDA0001937687870000074
进一步的，在两自由度系统中所述步骤1的具体方法为：
构建两自由度再生颤振铣削动力学的时滞微分方程如公式18所示：
Figure BDA0001937687870000075
其中q1(t)和q2(t)分别表示在两自由度系统中x，y方向上刀具的坐标，hxx(t)与hyy(t)表示在两自由度系统中x，y方向上周期为T的切削力系数；hxy(t)表示x方向与y方向之间的耦合作用，hyx(t)表示y方向与x方向之间的耦合作用；ξ、mt、ωn分别为刀具的阻尼比、质量和固有频率；
令x1(t)＝[q1(t) q2(t)
Figure BDA0001937687870000076
]，将公式18变换为空间状态形式，如公式19所示
Figure BDA0001937687870000077
其中，
Figure BDA0001937687870000078
分别表示为q1(t)，q2(t)的一阶导数，
Figure BDA0001937687870000079
为x1(t)的一阶导数，A1为常数矩阵，B1(t)表示周期为T的系数矩阵，具体A1和B1(t)的表示如公式20、21所示：A1为常数矩阵，B1(t)表示周期为T的系数矩阵，具体A1和B1(t)的表示如公式20、21所示：
Figure BDA00019376878700000710
Figure BDA0001937687870000081
其中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所示：
Figure BDA0001937687870000082
其中，ζ、ωn、mt、w和q(t)分别是刀具的阻尼比、固有频率、质量、轴向切深和坐标，
Figure BDA0001937687870000083
为q(t)的一阶导数，
Figure BDA0001937687870000084
为q(t)的二阶导数，h(t)表示周期为T的切削力系数，h(t)的表达式如公式23所示。本实施例中周期T与时滞量Ts相等，Ts＝60/NΩ，N为铣刀齿数，本实施例中N的值为2，Ω为铣刀主轴转速，铣刀的一阶固有频率为922*2π，相对阻尼为0.11，质量为0.03993。
Figure BDA0001937687870000091
其中，st和sn分别是切向和法向的切削力系数，本实施例中分别取6*108和2*108；φj(t)是第j刀齿的角位移，表达式为φj(t)＝(2πΩ/60)t+2π(j-1)/N；函数g(φj(t))表达式为：
Figure BDA0001937687870000092
其中φ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，
Figure BDA0001937687870000093
Figure BDA0001937687870000094
得到公式23：
Figure BDA0001937687870000095
其中A0为常数矩阵
Figure BDA0001937687870000096
B(t)表示周期为T的系数矩阵，
Figure BDA0001937687870000097
步骤2：将刀齿切削周期T等分成m个小间隔，本实施例中m取值为40；则具有y＝m+1个节点，步长
Figure BDA0001937687870000098
将公式23在任意一个时间区间[tK，tK+1]内进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程，如公式24所示，其中tK表示第K个节点的时刻；所述K＝1,2,…,y；
Figure BDA0001937687870000099
其中X(δ)为状态项、X(δ-T)为时滞项、B(δ)为周期系数项，Xk+1表示时刻tK+1的响应；Xk表示时刻tK的响应；δ为时间变量，0≤δ≤τ或者tK≤δ≤tK+1。
步骤3：在当前周期的时间区域[tK-2，tK+1]中，利用牛顿插值多项式算法将该时间区域中的四个时刻tK-2、tK-1、tK、tK-2以及这四个时刻的响应XK-2、XK-1、XK、XK+1拟合状态项，K>2，得到状态项的X(δ)的方程表达式为：
X(δ)＝a1XK+1+b1XK+c1XK-1+d1XK-2 (25)
其中，a1、b1、b1、d1为系数，具体的表示如公式26所示；
Figure BDA0001937687870000101
在上一个周期的时间区域[tK，-m，tK+2，-m]中，利用分段三次埃米尔特插值多项式算法将三个时刻tK，-m、tK+1，-m、tK+2，-m以及这三个时刻的响应XK，-m、XK+1，-m、XK+2，-m逼近时滞项，得到时滞项X(δ-T)的方程表达式如公式27所示，所述tK，-m、tK+1，-m、tK+2，-m分别表示上一个周期中节点K、K+1、K+2的时间：
X(δ-T)＝a2XK，-m+b2XK+1，-m+c2XK+2，-m (27)
其中a2、b2、c2为系数，具体表达式如公式28所示；
Figure BDA0001937687870000102
用线性插值算法逼近周期系数项，得到周期系数项B(δ)的方程表达式为：
Figure BDA0001937687870000103
其中BK+1表示在当前周期中第K+1个节点的周期系数；BK表示当前周期中第K个节点的周期系数。
步骤4：将上述公式25、27、29代入公式24中，整理得到公式30：
(Hk，1—I)XK+1+(H0+Hk，0)XK+Hk，-1XK-1+Hk，-2XK-2
＝Hk，mXK，-m+Hk，m-1XK+1，-m+Hk，m-2XK+2，-m (30)
其中，
Figure BDA0001937687870000104
I为单位矩阵，Hk，1为当前周期中第K+1个节点响应的切削力系数；Hk，0为当前周期中第K个节点响应的切削力系数，Hk，-1为当前周期中第K-1个节点响应的切削力系数；Hk，-2为当前周期中第K-2个节点响应的切削力系数；Hk，m-2为上一个周期中第K-2个节点响应的切削力系数；Hk，m-1为上一个周期中第K-1个节点响应的切削力系数；Hk，m为上一个周期中第K个节点响应的切削力系数；
Figure BDA0001937687870000111
Figure BDA0001937687870000112
Figure BDA0001937687870000113
Figure BDA0001937687870000114
Figure BDA0001937687870000115
Figure BDA0001937687870000116
Figure BDA0001937687870000117
其中，
Figure BDA0001937687870000118
Figure BDA0001937687870000119
Figure BDA00019376878700001110
Figure BDA00019376878700001111
Figure BDA00019376878700001112
步骤B：将上述公式30进行m次的迭代计算，得到相邻周期刀齿响应之间的映射关系，用矩阵表示如公式31所示：
Figure BDA00019376878700001113
其中Xm为当前周期的中时刻tm的响应；Xm，-m为上一个周期中时刻tm，-m的响应；D1和D2表达式如下所示：
Figure BDA0001937687870000121
Figure BDA0001937687870000122
其中，H1，0表示当前周期中第一个节点响应的切削力系数；H1，1表示当前周期中第1+1个节点即第二个节点响应的切削力系数。
构建相邻周期刀齿响应的传递矩阵Ф，即Floquet转移矩阵，如公式32所示：
Φ＝(D1)-1D2 (32)
步骤6：根据上述相邻周期刀齿响应传递矩阵Ф的特征值，模判定铣削系统的稳定性，当Ф最大特征值的模小于1时，该铣削动力学系统处于稳定状态；Ф最大特征值的模等于1时，该铣削动力学系统临界稳定状态，当Ф最大特征值的模大于1时，该铣削动力学系统处于不稳定状态；具体的判定如公式33所示：
Figure BDA0001937687870000123
其中λ()的含义为求特征值。
上述步骤3中用分段三次埃米尔特插值多项式算法逼近时滞项，得到时滞项X(δ-T)的方程表达式的具体方法为：
步骤3.1:在上一个周期的每个时间区域[tK，-m，tK+1，-m-m]内，构造分段三次埃尔米特插值多项式，得到公式34：
Figure BDA0001937687870000124
其中，
Figure BDA0001937687870000125
为响应XK，-m的一阶导数，
Figure BDA0001937687870000126
为响应XK+1，-m的一阶导数；其中系数ε0(δ)、ε1(δ)、β0(δ)、β1具体的表述如公式35所示：
Figure BDA0001937687870000131
当δ＝tK，-m时，
Figure BDA0001937687870000132
当δ＝tK+1，-m时，
Figure BDA0001937687870000133
步骤3.2：由于
Figure BDA0001937687870000134
令τ＝tK+1，-m—tK，-m并将上述公式36和公式37代入公式34中；得到X(δ-T)的表达式如公式38所示：
Figure BDA0001937687870000135
在两自由度的系统中，该预测方法的步骤1的具体方法为：
建立两自由度再生颤振铣削动力学的时滞微分方程如公式39所示：
Figure BDA0001937687870000136
其中q1(t)和q2(t)分别表示在两自由度系统中x，y方向的坐标，hxx(t)与hyy(t)表示在两自由度系统中x，y方向的切削力系数；hxy(t)表示x方向与y方向之间的耦合作用，hyx(t)表示y方向与x方向之间的耦合作用；ξ、mt、ωn分别为刀具的阻尼比、质量和固有频率；
Figure BDA0001937687870000137
Figure BDA0001937687870000141
Figure BDA0001937687870000142
Figure BDA0001937687870000143
令x1(t)＝[q1(t) q2(t)
Figure BDA0001937687870000144
]，将公式18变换为空间状态形式，如公式44所示
Figure BDA0001937687870000145
其中，
Figure BDA0001937687870000146
分别表示为q1(t)，q2(t)的一阶导数，
Figure BDA0001937687870000147
为x1(t)的一阶导数，A1为常数矩阵，B1(t)表示周期为T的系数矩阵，具体A1和B1(t)的表示如公式45、46所示：其中，
Figure BDA0001937687870000148
Figure BDA0001937687870000149
将两自由度的动力学方程变换为空间状态形式，再进行上述步骤2至步骤6的操作后得到两自由度的稳定域的准确预测。
本实施例利用Matlab软件进行编程，画出单自由度和两自由度的稳定性叶瓣图及收敛图，通过稳定性图来预测铣削过程中地稳定性，选取不同地径向浸入比，分别取0.1和1得到的稳定性叶瓣图如图1、2、3和4所示及选取不同轴向切深w，分别取0.2mm和0.5mm得到的收敛图如图5和6所示。
另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

Claims (3)

1.一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法，其特征在于，具体包括如下步骤:
步骤1：将考虑再生颤振铣削动力学的时滞微分方程变换为空间状态形式，得到该时滞微分方程的空间状态方程；
步骤2：将刀齿切削周期等分成若干个小间隔，在任意一个小间隔的时间长度内对上述时滞微分方程的空间状态方程进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程；
步骤3：利用插值多项式算法逼近上述的状态项、时滞项、周期系数项，从而得到状态项、时滞项、周期系数项的方程表达式；并将状态项、时滞项、周期系数项的方程表达式代入步骤2中的状态项、时滞项、周期系数项三者之间的关系方程中，并利用迭代法计算得到相邻周期刀齿响应之间的映射关系；根据所述的映射关系构建相邻周期刀齿响应的传递矩阵Ф；
步骤4：计算步骤3中的相邻周期刀齿响应传递矩阵Ф的特征值，判定铣削系统的稳定性，具体的判定如公式1所示：
Figure FDA0003649399400000011
其中λ()的含义为求特征值；
在单自由度系统中所述步骤1的具体方法为：
步骤1.1：构建考虑再生颤振的单自由度铣削动力学的时滞微分方程，如公式2所示：
Figure FDA0003649399400000012
其中ξ、mt、ωn分别为刀具的阻尼比、质量和固有频率；h(t)表示周期为T的切削力系数；q(t)为刀具的坐标，
Figure FDA0003649399400000013
为q(t)的一阶导数，
Figure FDA0003649399400000014
为q(t)的二阶导数；w为轴向切深；
步骤1.2：令X(t)＝[q(t) p(t)]T，其中
Figure FDA0003649399400000015
将公式2转换成空间状态形式如公式3所示：
Figure FDA0003649399400000016
其中A0为常数矩阵
Figure FDA0003649399400000021
Figure FDA0003649399400000022
为X(t)的一阶导数；B(t)表示周期为T的系数矩阵，
Figure FDA0003649399400000023
所述步骤2的具体方法为：将刀齿切削周期T等分成m个小间隔，则具有y＝m+1个节点，步长
Figure FDA0003649399400000024
将公式2在任意一个时间区间[tK,tK+1]内进行积分计算，得到状态项、时滞项、周期系数项三者之间的关系方程，如公式4所示，其中tK表示第K个节点的时刻；所述K＝1,2,…,y；
Figure FDA0003649399400000025
其中X(δ)为状态项、X(δ-T)为时滞项、B(δ)为周期系数项，Xk+1表示时刻tK+1的响应；Xk表示时刻tK的响应；δ为时间变量，0≤δ≤τ或者tK≤δ≤tK+1；
所述步骤3中确定状态项、时滞项、周期系数项的方程表达式的具体方法为：
在当前周期的时间区域[tK-2,tK+1]中，利用牛顿插值多项式算法将该时间区域中的四个时刻tK-2、tK-1、tK、tK+1以及这四个时刻的响应XK-2、XK-1、XK、XK+1逼近状态项，K>2，得到状态项X(δ)的方程表达式为：
X(δ)＝a1XK+1+b1XK+c1XK-1+d1XK-2 (5)
其中，a1、b1、b1、d1为系数，具体的表示如公式6所示；
Figure FDA0003649399400000026
在上一个周期的时间区域[tK，-m,tK+2，-m]中，利用分段三次埃米尔特插值多项式算法将三个时刻tK，-m、tK+1，-m、tK+2，-m以及这三个时刻的响应XK,-m、XK+1,-m、XK+2,-m逼近时滞项，得到时滞项X(δ-T)的方程表达式如公式7所示，所述tK，-m、tK+1，-m、tK+2，-m分别表示上一个周期中节点K、K+1、K+2的时刻：
X(δ-T)＝a2XK,-m+b2XK+1,-m+c2XK+2,-m (7)
其中a2、b2、c2为系数，具体表达式如公式8所示；
Figure FDA0003649399400000031
用线性插值算法逼近周期系数项，得到周期系数项B(δ)的方程表达式为：
Figure FDA0003649399400000032
其中BK+1表示在当前周期中第K+1个节点的周期系数；BK表示当前周期中第K个节点的周期系数；
步骤3中的构建相邻周期刀齿响应的传递矩阵Ф具体方法如下：
步骤A：将公式5、7、9代入状态项、时滞项、周期系数项三者之间的关系方程中，整理得到公式10；
(Hk,1-I)XK+1+(H0+Hk,0)XK+Hk,-1XK-1+Hk,-2XK-2
＝Hk,mXK,-m+Hk,m-1XK+1,-m+Hk,m-2XK+2,-m (10)
其中，
Figure FDA0003649399400000033
I为单位矩阵，Hk,1为当前周期中第K+1个节点响应的切削力系数；Hk,0为当前周期中第K个节点响应的切削力系数，Hk,-1为当前周期中第K-1个节点响应的切削力系数；Hk,-2为当前周期中第K-2个节点响应的切削力系数；Hk,m-2为上一个周期中第K-2个节点响应的切削力系数；Hk,m-1为上一个周期中第K-1个节点响应的切削力系数；Hk,m为上一个周期中第K个节点响应的切削力系数；
Figure FDA0003649399400000034
Figure FDA0003649399400000035
Figure FDA0003649399400000036
Figure FDA0003649399400000037
Figure FDA0003649399400000041
Figure FDA0003649399400000042
Figure FDA0003649399400000043
其中，
Figure FDA0003649399400000044
Figure FDA0003649399400000045
Figure FDA0003649399400000046
Figure FDA0003649399400000047
Figure FDA0003649399400000048
步骤B：将上述公式10进行m次的迭代计算，得到相邻周期刀齿响应之间的映射关系，用矩阵表示如公式11所示：
Figure FDA0003649399400000049
其中Xm为当前周期的中时刻tm的响应；Xm,-m为上一个周期中时刻tm,-m的响应；D1和D2表达式如下所示：
Figure FDA00036493994000000410
Figure FDA00036493994000000411
步骤C：构建相邻周期刀齿响应的传递矩阵Ф，如公式12所示：
Ф＝(D1)-1D2 (12)。
2.基于权利要求1所述的预测方法，其特征在于，上述用分段三次埃米尔特插值多项式算法逼近时滞项，得到时滞项X(δ-T)的方程表达式的具体方法为：
步骤S1:在上一个周期的每个时间区域[tK,-m,tK+1,-m]内，构造分段三次埃尔米特插值多项式，得到公式13：
Figure FDA0003649399400000051
其中，
Figure FDA0003649399400000052
为响应XK,-m的一阶导数，
Figure FDA0003649399400000053
为响应XK+1,-m的一阶导数；系数ε0(δ)、ε1(δ)、β0(δ)、β1具体的表述如公式14所示：
Figure FDA0003649399400000054
当δ＝tK,-m时，
Figure FDA0003649399400000055
当δ＝tK+1,-m时，
Figure FDA0003649399400000056
步骤S2：由于
Figure FDA0003649399400000057
令τ＝tK+1,-m—tK,-m；并将上述公式15和公式16代入公式13中；得到X(δ-T)的表达式如公式17所示：
Figure FDA0003649399400000058
3.基于权利要求1所述的预测方法，其特征在于，在两自由度系统中所述步骤1的具体方法为：
构建两自由度再生颤振铣削动力学的时滞微分方程如公式18所示：
Figure FDA0003649399400000061
其中q1(t)和q2(t)分别表示在两自由度系统中x，y方向上刀具的坐标，hxx(t)与hyy(t)表示在两自由度系统中x，y方向上周期为T的切削力系数；hxy(t)表示x方向与y方向之间的耦合作用，hyx(t)表示y方向与x方向之间的耦合作用；ξ、mt、ωn分别为刀具的阻尼比、质量和固有频率；
令
Figure FDA0003649399400000062
将公式18变换为空间状态形式，如公式19所示
Figure FDA0003649399400000063
其中，
Figure FDA0003649399400000064
分别表示为q1(t)，q2(t)的一阶导数，
Figure FDA0003649399400000065
为x1(t)的一阶导数，A1为常数矩阵，B1(t)表示周期为T的系数矩阵，具体A1和B1(t)的表示如公式20、21所示：
Figure FDA0003649399400000066
Figure FDA0003649399400000067
其中w为轴向切深。
CN201910012032.6A 2019-01-07 2019-01-07 一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法 Active CN109740264B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910012032.6A CN109740264B (zh) 2019-01-07 2019-01-07 一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910012032.6A CN109740264B (zh) 2019-01-07 2019-01-07 一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法

Publications (2)

Publication Number Publication Date
CN109740264A CN109740264A (zh) 2019-05-10
CN109740264B true CN109740264B (zh) 2022-08-05

Family

ID=66363603

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910012032.6A Active CN109740264B (zh) 2019-01-07 2019-01-07 一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法

Country Status (1)

Country Link
CN (1) CN109740264B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111158315B (zh) * 2019-11-22 2023-03-10 河南理工大学 一种基于样条-牛顿公式的铣削稳定性预测方法
CN111611725B (zh) * 2020-06-18 2022-05-13 南昌航空大学 一种基于Cotes数值积分的铣削稳定域预测方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104657606A (zh) * 2015-02-10 2015-05-27 北京理工大学 一种基于三次多项式的铣削稳定性预测方法
CN104680248A (zh) * 2015-02-10 2015-06-03 北京理工大学 一种基于勒让德多项式的铣削稳定性预测方法
CN106156477A (zh) * 2015-04-28 2016-11-23 河南理工大学 薄壁件动态铣削稳定性叶瓣图高精度预测方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104657606A (zh) * 2015-02-10 2015-05-27 北京理工大学 一种基于三次多项式的铣削稳定性预测方法
CN104680248A (zh) * 2015-02-10 2015-06-03 北京理工大学 一种基于勒让德多项式的铣削稳定性预测方法
CN106156477A (zh) * 2015-04-28 2016-11-23 河南理工大学 薄壁件动态铣削稳定性叶瓣图高精度预测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
An updated full-discretization milling stability prediction method based on the higher-order Hermite-Newton interpolation polynomial;Yongjian Ji等;《The International Journal of Advanced Manufacturing Technology》;20171121;第2227–2242页 *
Hyper-third order full-discretization methods in milling stability prediction;C.G. Ozoegwu等;《International Journal of Machine Tools & Manufacture》;20150531;第1-9页 *
Third-order updated full-discretization method for milling stability prediction;Zhenghu Yan等;《The International Journal of Advanced Manufacturing Technology》;20170331;第2299-2309页 *

Also Published As

Publication number Publication date
CN109740264A (zh) 2019-05-10

Similar Documents

Publication Publication Date Title
CN110188311B (zh) 基于刀齿切削时程精细积分的高速加工稳定域预测方法
CN109740264B (zh) 一种利用牛顿和埃尔米特插值法的铣削稳定域预测方法
Li et al. Prediction of chatter stability for milling process using Runge-Kutta-based complete discretization method
Tsai et al. A real-time predictor-corrector interpolator for CNC machining
CN106156477B (zh) 薄壁件动态铣削稳定性叶瓣图高精度预测方法
CN104647132B (zh) 一种基于磁悬浮轴承电主轴的铣削颤振主动控制方法
CN108983615B (zh) 基于反双曲正弦吸引律的离散双周期重复控制器
CN108958041B (zh) 一种基于双曲正割吸引律的离散双周期重复控制方法
CN106843147B (zh) 一种基于Hamming公式预测铣削稳定性的方法
Huang et al. An efficient linear approximation of acceleration method for milling stability prediction
Qin et al. A predictor-corrector-based holistic-discretization method for accurate and efficient milling stability analysis
Pourmahmood Aghababa Control of fractional-order systems using chatter-free sliding mode approach
CN106647625A (zh) 一种基于Gear公式预测铣削稳定性的方法
Ding et al. Milling stability analysis using the spectral method
CN110064965B (zh) 一种铣削系统稳定性状态获取方法
Tao et al. Milling Stability Prediction with Multiple Delays via the Extended Adams‐Moulton‐Based Method
CN110162733B (zh) 基于整体离散策略的铣削稳定性分析方法
CN111176209B (zh) 型腔螺旋铣削加工进给率与转速离线规划方法
Huang et al. An Efficient Third‐Order Full‐Discretization Method for Prediction of Regenerative Chatter Stability in Milling
CN112016203B (zh) 基于分段Hermite插值多项式和整体离散策略预测铣削稳定性的方法
CN104852491A (zh) 一种永磁同步电机磁钢分段斜极角的设计方法
CN114818182A (zh) 一种响应驱动的薄壁结构宽频减振动力吸振器参数设计方法
CN111158315B (zh) 一种基于样条-牛顿公式的铣削稳定性预测方法
CN114895566B (zh) 一种采用传递矩阵缩减技术的铣削过程的颤振预测方法
CN111914368A (zh) 考虑螺旋角效应的变齿距变转速铣削颤振主被动抑制方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant