铣削精加工的稳定性与表面位置误差综合预测方法及设备

1.本发明涉及机械加工中的铣削精加工技术领域，更具体的说是涉及铣削精加工的稳定性与表面位置误差综合预测方法及设备。


背景技术：

2.铣削加工技术被广泛应用于航空、航天、船舶和汽车等领域的复杂零件制造过程中，是获取高精度加工质量的重要基础技术之一。实际加工过程中，由于加工参数选取不合理而导致的加工过程的颤振现象是降低零件表面质量及限制生产效率的主要原因之一。因此，铣削过程的稳定性分析是实现无颤振铣削的重要前提。但是，由于存在强迫振动，铣削过程的无颤振加工参数仍然不足以保证达到所期望的加工精度，尤其在精加工过程中应优先考虑被加工零件的表面质量，在满足表面质量要求的前提下，再尽可能保证最大加工效率。
3.目前，本领域相关技术人员已经做了一些研究，如文献“li zy,jiang sl,sun yw.chatter stability and surface location error predictions in milling with mode coupling and process damping[j].proceedings of the institution of mechanical engineers,part b:journal ofengineering manufacture,2019,233(3):686-698.”公开了一种考虑模态耦合与过程阻尼的铣削动力学建模方法，并利用二阶半离散法同时预测铣削稳定性和表面位置误差，但是其在考虑刀具跳动等情况时，无法进行准确的稳定性边界和表面位置误差预测，预测精度有限；又如文献“wan m,zhang wh,dang jw,et al.a unified stability prediction method for milling process with multiple delays[j].international journal of machine tools and manufacture,2010,50(1):29-41.”公开了一种考虑刀具跳动的多时滞铣削稳定性预测模型，仅用于预测铣削稳定性，而且该方法未考虑铣削过程模态耦合效应与过程阻尼，预测结果误差较大。
[0004]
目前公开号为：“cn106808320a”的专利，公开了一种考虑刀具齿间角、螺旋角及偏心的铣削力预测方法，仅用于预测切削力；公开号为：“cn108647413a”的专利，公开了一种微细表面位置误差与稳定性综合预测方法，未考虑铣削过程中的模态耦合效应、过程阻尼和刀具跳动的影响，较实际结果预测误差较大。
[0005]
综上所述，现有的铣削稳定性和表面位置误差预测模型大都仅考虑铣削过程中的再生效应，而实际上，铣削过程往往伴随多种复杂因素的影响，大量研究表明基于传统动力学模型的稳定性和表面位置误差预报和实验结果之间存在着较大的偏差，而偏差的产生很大程度上是因为忽略了模态耦合效应、过程阻尼和刀具跳动。
[0006]
因此，针对现有技术的以上缺陷或改进需求，如何提供一种铣削精加工的稳定性与表面位置误差综合预测方法及设备是本领域技术人员亟需解决的问题。


技术实现要素：

[0007]
有鉴于此，本发明提供了一种铣削精加工的稳定性与表面位置误差综合预测方法
及设备，用以解决上述现有技术中存在的技术问题。
[0008]
为了实现上述目的，本发明提供如下技术方案：
[0009]
一方面，本发明实施例1公开了一种铣削精加工的稳定性与表面位置误差综合预测方法，包括以下步骤：
[0010]
步骤1：获取铣刀信息，所述铣刀信息包括：刀具几何参数、模态参数、切削力系数及切削条件；
[0011]
步骤2：基于铣刀信息计算任意时刻切削刃角度位置、切削半径、进给量及在x和y方向总切削力；
[0012]
步骤3：构建铣削动力学时滯微分方程组，并转化为空间状态方程；
[0013]
步骤4：将刀齿切削周期t等分成
[0014]
将刀齿切削周期等分成若干个小间隔，在任意一个小间隔的时间长度内对所述空间状态方程进行积分计算，确定主轴转速范围及其对应的离散数、轴向切深范围及其对应离散数、单个刀齿切削周期离散数；
[0015]
步骤5：利用三阶牛顿、埃尔米特和线性插值多项式近似计算所述空间状态方程中的状态项、时滞项和静态力项，构建铣削精加工系统在单个时间周期上的状态传递关系，判定铣削精加工系统的稳定性；
[0016]
步骤6：提取出相应的稳态系数向量y方向系数，预报动态加工误差；
[0017]
步骤7：基于步骤5和步骤6，获得稳定性叶瓣图和表面位置误差预测图。
[0018]
优选的，所述步骤2基于铣刀信息计算任意时刻切削刃角度位置、切削半径、进给量及在x和y方向总切削力的具体步骤为：
[0019]
步骤2.1：计算任意时刻切削刃的刀具角度位置；
[0020]
步骤2.2：计算刀具偏心情况下切削刃实际每齿的切削半径和每齿进给量；
[0021]
步骤2.3：计算在刀具切向和径向的切削力，并通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上；
[0022]
步骤2.4：通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上，根据累加作用在刀齿上的切削力计算作用在铣刀上x和y方向的总切削力。
[0023]
优选的，计算任意时刻切削刃的刀具角度位置，具体公式如下：
[0024][0025]
式中，ω为主轴转速；n为刀齿数；j为刀齿编号。
[0026]
优选的，计算刀具偏心情况下切削刃实际每齿的切削半径和每齿进给量，具体步骤如下：
[0027]
每个齿的实际切削半径为：
[0028][0029][0030][0031]
实际的每齿进给为：
[0032][0033][0034][0035]
式中，f
t
为名义每齿进给量，ar为径向切削深度；
[0036]
根据几何关系，相邻两个刀齿之间的夹角由余弦定理可获得：
[0037][0038][0039][0040]
第j齿的切入切出角表示为：
[0041][0042][0043]
通过上述得到在具偏心的情况下，切削刃实际每齿的切削半径和每齿进给量。
[0044]
计算在刀具切向和径向的切削力，并通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上，具体步骤如下：
[0045]
刀具切向和径向的切削力计算：
[0046]
作用在刀齿j上的铣削力主要包括剪切力和过程阻尼力，分别表示为
[0047]ft
＝f
ts
+f
tp
,fr＝f
rs
+f
rp
[0048]
式中，f
ts
和f
rs
为剪切效应引起的切向力和法向力，f
tp
和f
rp
为过程阻尼效应引起的切向力和法向力。基于剪切力模型，f
ts
和f
rs
分别表示为
[0049]fts
＝g(φj(t))k
tap
h(t),f
rs
＝g(φj(t))k
nap
h(t)
[0050]
式中，k
t
和kn分别为切向与径向切削力系数，ap为轴向切削深度，窗函数g(φj(t))定义为
[0051][0052]
式中，φ
st
和φ
ex
分别是刀齿的切入与切出角。基于再生颤振理论，瞬时切削厚度可以表示为：
[0053]
h(t)＝[f
t
+x(t-t)-x(t)]sin(φj(t))+[y(t-t)-y(t)]cos(φj(t))
[0054]
式中，f
t
为每齿进给量，t为时滞量且等于刀齿切削周期，即t＝60/(nω)。
[0055]
基于过程阻尼机理，径向过程阻尼力可以表示为：
[0056]frp
＝g(φj(t))k
spaps[0057]
式中，k
sp
和s分别表示压入系数和面积。过程阻尼效应产生的切向力可以从下式获得：
[0058]ftp
＝μf
rp
[0059]
式中，μ为摩擦系数。等效粘性阻尼可以表示为：
[0060][0061]
式中，v是切削速度，w是刀具磨损带宽度。
[0062]
通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上，根据累加作用在刀齿上的切削力计算作用在铣刀上x和y方向的总切削力。
[0063]
优选的，所述步骤3构建铣削动力学时滯微分方程组，并转化为空间状态方程的具体方法为：
[0064]
步骤3.1：构建铣削动力学时滯微分方程组：
[0065][0066]
式中，
[0067][0068]
步骤3.2：令借助柯西变换，转为空间状态形式，得到空间状态方程为：
[0069][0070]
式中，
[0071][0072]
式中，m为模态质量矩阵，m
xx
、m
yy
为x和y向的质量，m
xy
、m
yx
为x、y向在此坐标下的耦合质量。c为模态阻尼矩阵，c
xx
、c
yy
为x和y向的阻尼，c
xy
、c
yx
为x、y向在此坐标下的耦合阻尼。k为模态刚度矩阵，k
xx
、k
yy
为x和y向的刚度，k
xy
、k
yx
为x、y向在此坐标下的耦合刚度。c
eq
为等效粘性阻尼，k
sp
为压入系数，a
p
为轴向切削深度，w为刀具磨损带宽，v为切削速度。c
p
、h为切削力的转化系数矩阵，有：
[0073][0074][0075][0076][0077]
f0(t)为静态力向量，f
t
为每齿进给量。
[0078]
优选的，所述步骤4将刀齿切削周期等分成若干个小间隔，在任意一个小间隔的时
间长度内对所述空间状态方程进行积分计算的具体方法为：
[0079]
定义κ(t)＝a(t)x(t)和α(t-t)＝b(t)x(t-t)；
[0080]
以x(t0)为初始条件，空间状态方程的通解表示为：
[0081][0082]
将刀齿切削周期t等分为m个时间间隔，则时间步长τ＝t/m；上式在时间kτ≤t≤(k+1)τ，(k＝1,2,3,
…
m)区段上可表示为：
[0083][0084]
式中，δ＝ξ-kτ,δ∈[0,τ]；式中，x((k+1)τ)记作x
k+1
。
[0085]
将刀齿切削周期t等分成为m个时间间隔，则时间步长τ＝t/m；上式在时间kτ≤t≤(k+1)τ，(k＝1,2,3,
…
m)区段上可表示为：
[0086][0087]
式中，δ＝ξ-kτ,δ∈[0,τ]；式中，x((k+1)τ)记作x
k+1
。
[0088]
根据机床约束(如允许最大切削功率、允许最大切削力)、刀具约束(如刀具的材料和刀具耐用度)和零件约束(如被加工零件的精度、材料)等几大类约束来确定允许最高切削速度，进而按照允许最高切削速度来确定铣削加工工艺系统可以承受的允许最低主轴转速n
min
、最高主轴转速n
max
，从而可得主轴转速范围[n
min
,n
max
]。主轴转速范围的离散数kn通常设为200。
[0089]
根据机床刚度、刀具直径和加工精度要求来确定允许最大轴向切深，进而按照在粗加工时通常要求被加工的零件的材料去除率高而采用较大轴向切深和精加工时通常要求被加工的零件的加工精度要求高而采用较小轴向切深来确定可以接受的允许最低主轴转速n
min
、最高主轴转速n
max
，从而可得轴向切深范围[w
min
,w
max
]。轴向切深范围的离散数kw通常设为100。
[0090]
刀齿切削周期t的离散数m通常设为20、30、40、50、60、70、80、90、100、200等10个档次。
[0091]
颤振稳定性求解过程中可兼顾对计算精度和计算效率的实际需求情况减设或增设主轴转速范围的离散数、轴向切深范围的离散数和单个刀齿切削周期的离散数。
[0092]
优选的，所述步骤5计算所述空间状态方程中的状态项、时滞项和静态力项，构建铣削精加工系统在单个时间周期上的状态传递关系，判定铣削精加工系统的稳定性的具体步骤为：
[0093]
步骤5.1：利用三阶牛顿、埃尔米特和线性插值多项式近似分别计算所述空间状态方程中的状态项、时滞项和静态力项；
[0094]
步骤5.2：将状态项、时滞项和静态力项带入所述空间状态方程，得到如下离散动态映射：
[0095]
(f
k+1-i)x
k+1
+(t1+fk)xk+f
k-1
x
k-1
+f
k-2
x
k-2
＝f
k-m
x
k-m
+f
k-m+1
x
k-m+1
+f
k-m+2
x
k-m+2-g
kfk-g
k+1fk+1
[0096]
步骤5.3：使用扩展的状态向量zk＝[x1,x2,
…
,x
m+1
]
t
，yk＝[f1,f2,
…
,f
m+1
]
t
，则离散动态映射表示为：
[0097]
d1zk＝d2z
k-m
+d3yk[0098]
式中，
[0099][0100][0101][0102]
步骤5:.4：建立铣削精加工系统在单个时间周期上状态传递关系；
[0103]
zm＝φz0+d3yk[0104]
式中，
[0105]
φ＝(d1)-1
d2[0106]
基于floquet理论，铣削精加工系统稳定性由传递矩阵φ的特征值的模决定，具体公式为：
[0107][0108]
优选的，所述步骤5.1利用三阶牛顿插值多项式近似计算所述空间状态方程中的状态项的具体公式为：
[0109]
κ(δ)由状态向量在t
k-2
、t
k-1
、tk和t
k+1
时刻的响应κ
k-2
、κ
k-1
、κk和κ
k+1
近似为：
[0110]
κ(δ)＝a1κ
k+1
+b1κk+c1κ
k-1
+d1κ
k-2
[0111]
式中，κ
k+1
＝a
k+1
x
k+1
,κk＝akxk,κ
k-1
＝a
k-1
x
k-1
,κ
k-2
＝a
k-2
x
k-2
，式中t
k-2
＝(k-2)τ，κ((k+1)τ)记作κ
k+1
＝a((k+1)τ)x((k+1)τ)，a((k+1)τ)记作a
k+1
，插值系数a1、b1、c1和d1分别为：
[0112][0113]
[0114][0115][0116]
优选的，所述步骤5.1利用埃尔米特插值多项式近似计算所述空间状态方程中的时滞项的具体公式为：
[0117]
α(δ-t)＝a2α
k-m
+b2α
k-m+1
+c2α
k-m+2
[0118]
式中，α
k-m
＝bkx
k-m
α
k-m+1
＝b
k+1
x
k-m+1
α
k-m+2
＝b
k+2
x
k-m+2
，式中α(k-m-t)记作α
k-m
＝b(k-m)x(k-m-t)，b(kτ)记作bk，插值系数a2、b2和c2分别为：
[0119][0120][0121][0122]
优选的，所述步骤5.1利用线性插值多项式近似计算所述空间状态方程中的静态力项的具体方法为：
[0123][0124]
式中，f(kτ)记作fk。
[0125]
优选的，所述步骤6：提取出相应的稳态系数向量y方向系数，预报动态加工误差的具体方法为：
[0126]
步骤6.1：对于多周期情况，第i周期和第(i-1)周期的状态传递关系表达式为：
[0127]
zi＝φz
i-1
+d3yk[0128]
步骤6.2：在稳态切削时，稳态系数为：
[0129]
zk＝z
k-m
＝z*
[0130]
步骤6.3：将稳态系数)代入到状态传递关系表达式得到：
[0131]z*
＝(d
1-d2)-1
d3yk[0132]
步骤6.4：从步骤6.3的公式中提取出相应的稳态系数向量y方向系数，预报动态加工误差。
[0133]
另一方面，本发明还提供了一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现铣削精加工的稳定性与表面位置误差综合预测方法。
[0134]
经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种铣削精加工的稳定性与表面位置误差综合预测方法及设备，用于铣削精加工的切削稳定性与表面位置误差综合预测，考虑了铣削过程中的模态耦合、过程阻尼和刀具跳动效应的影响提高了解析模型的精确性，同时，通过建立的考虑模态耦合、过程阻尼和刀具跳动的稳定性叶瓣图和表面位置误差图，能够方便快捷的根据图形预测加工状态和加工误差，为以考虑无颤振和动态加工误差小于期望值为约束的最优选取加工参数提供理论依据。
附图说明
[0135]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
[0136]
图1为本发明的方法流程示意图；
[0137]
图2为本发明考虑模态耦合、过程阻尼和刀具跳动的稳定性叶瓣图；
[0138]
图3为本发明铣削加工稳定性和表面位置误差综合预测图；
[0139]
图4为本发明不同偏心距的稳定性叶瓣图；
[0140]
图5为本发明不同偏心角的稳定性叶瓣图。
具体实施方式
[0141]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0142]
参见附图1所示，本发明实施例公开了一种铣削精加工的稳定性与表面位置误差综合预测方法，基于铣削精系统，包括以下步骤：
[0143]
步骤1：获取铣刀信息，铣刀信息包括：刀具几何参数、模态参数、切削力系数及切削条件；
[0144]
步骤2：基于铣刀信息计算任意时刻切削刃角度位置、切削半径、进给量及在x和y方向总切削力；
[0145]
步骤3：构建铣削动力学时滯微分方程组，并转化为空间状态方程；
[0146]
步骤4：将刀齿切削周期等分成若干个小间隔，在任意一个小间隔的时间长度内对空间状态方程进行积分计算，确定主轴转速范围及其对应的离散数、轴向切深范围及其对应离散数、单个刀齿切削周期离散数；
[0147]
步骤5：利用三阶牛顿、埃尔米特和线性插值多项式近似计算空间状态方程中的状态项、时滞项和静态力项，构建铣削精加工系统在单个时间周期上的状态传递关系，判定铣削精加工系统的稳定性；
[0148]
步骤6：提取出相应的稳态系数向量y方向系数，预报动态加工误差；
[0149]
步骤7：基于步骤5和步骤6，获得稳定性叶瓣图和表面位置误差预测图。
[0150]
在一个具体实施例中，步骤2基于铣刀信息计算任意时刻切削刃角度位置、切削半径、进给量及在x和y方向总切削力的具体步骤为：
[0151]
步骤2.1：计算任意时刻切削刃的刀具角度位置；
[0152]
步骤2.2：计算刀具偏心情况下切削刃实际每齿的切削半径和每齿进给量；
[0153]
步骤2.3：计算在刀具切向和径向的切削力，并通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上；
[0154]
步骤2.4：通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上，根据累加作用在刀齿上的切削力计算作用在铣刀上x和y方向的总
切削力。
[0155]
具体的，针对具有n条切削刃的铣刀，第j个刀齿上的瞬时径向接触角φj(j＝1...n)可以表示为时间t和主轴转速ω的函数关系式：
[0156][0157]
刀具制造误差和刀具在机床上的安装误差将会导致刀具轴线与机床主轴线不重合和引起刀具偏心，多齿铣刀旋转加工存在刀具偏心。刀具偏心的存在将会改变实际每齿的切削半径和每齿进给量。
[0158]
具体的，计算刀具偏心情况下切削刃实际每齿的切削半径和每齿进给量，具体步骤如下：
[0159]
每个齿的实际切削半径为：
[0160][0161]
实际的每齿进给为：
[0162][0163]
式中，f
t
为名义每齿进给量，ar为径向切削深度。
[0164]
根据几何关系，相邻两个刀齿之间的夹角由余弦定理可获得：
[0165][0166]
第j齿的切入切出角表示为：
[0167][0168]
具体的，计算在刀具切向和径向的切削力，并通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上。作用在刀齿j上的铣削力主要包括剪切力和过程阻尼力，分别表示为：
[0169]ft
＝f
ts
+f
tp
,fr＝f
rs
+f
rp
ꢀꢀ
(6)
[0170]
式中，f
ts
和f
rs
为剪切效应引起的切向力和法向力，f
tp
和f
rp
为过程阻尼效应引起的切向力和法向力。基于剪切力模型，f
ts
和f
rs
分别表示为：
[0171]fts
＝g(φj(t))k
tap
h(t),f
rs
＝g(φj(t))k
nap
h(t)
ꢀꢀ
(7)
[0172]
式中，k
t
和kn分别为切向与径向切削力系数，a
p
为轴向切削深度，窗函数g(φj(t))定义为：
[0173][0174]
式中，φ
st
和φ
ex
分别是刀齿的切入与切出角。基于再生颤振理论，瞬时切削厚度可以表示为
[0175]
h(t)＝[f
t
+x(t-t)-x(t)]sin(φj(t))+[y(t-t)-y(t)]cos(φj(t))
ꢀꢀ
(9)
[0176]
式中f
t
为每齿进给量，t为时滞量且等于刀齿切削周期，即t＝60
/
(nω)。
[0177]
基于过程阻尼机理，径向过程阻尼力可以表示为
[0178]frp
＝g(φj(t))k
spapsꢀꢀ
(10)
[0179]
式中，k
sp
和s分别表示压入系数和面积。过程阻尼效应产生的切向力可以从下式获得
[0180]ftp
＝μf
rp
ꢀꢀ
(11)
[0181]
式中μ为摩擦系数。等效粘性阻尼可以表示为
[0182][0183]
式中，v是切削速度，w是刀具磨损带宽度。
[0184]
具体的，通过累加作用在刀齿上的切削力，并计算作用在铣刀上x和y方向的总切削力：
[0185][0186]
将公式(6)、(7)、(11)和(12)代入到公式(13)，切削合力可以改写为：
[0187][0188]
式中，
[0189][0190][0191][0192][0193]
在一个具体实施例中，综合考虑再生效应、过程阻尼、模态耦合效应和刀具跳动的刀具动力学方程表示为如下的时滞微分方程组：
[0194][0195]
式中，
[0196][0197]
令借助柯西变换，公式(19)可转为空间状态形式
[0198][0199]
式中，
[0200][0201]
在一个具体实施例中，为了简化算法推演过程，定义κ(t)＝a(t)x(t)和α(t-t)＝b(t)x(t-t)。以x(t0)为初始条件，公式(21)的通解表示为：
[0202][0203]
将刀齿切削周期t等分为m个时间间隔，则时间步长τ＝t/m。公式(23)在时间kτ≤t≤(k+1)τ，(k＝1,2,3,
…
m)区段上可以表示为
[0204][0205]
将刀齿切削周期t等分成为m个时间间隔，则时间步长τ＝t/m；上式在时间kτ≤t≤(k+1)τ，(k＝1,2,3,
…
m)区段上可表示为：
[0206][0207]
式中，δ＝ξ-kτ,δ∈[0,τ]；式中，x((k+1)τ)记作x
k+1
。
[0208]
根据机床约束(如允许最大切削功率、允许最大切削力)、刀具约束(如刀具的材料和刀具耐用度)和零件约束(如被加工零件的精度、材料)等几大类约束来确定允许最高切削速度，进而按照允许最高切削速度来确定铣削加工工艺系统可以承受的允许最低主轴转速n
min
、最高主轴转速n
max
，从而可得主轴转速范围[n
min
,n
max
]。主轴转速范围的离散数kn通常设为200。
[0209]
根据机床刚度、刀具直径和加工精度要求来确定允许最大轴向切深，进而按照在粗加工时通常要求被加工的零件的材料去除率高而采用较大轴向切深和精加工时通常要求被加工的零件的加工精度要求高而采用较小轴向切深来确定可以接受的允许最低主轴转速n
min
、最高主轴转速n
max
，从而可得轴向切深范围[w
min
,w
max
]。轴向切深范围的离散数kw通常设为100。
[0210]
刀齿切削周期t的离散数m通常设为20、30、40、50、60、70、80、90、100、200等10个档次。
[0211]
颤振稳定性求解过程中可兼顾对计算精度和计算效率的实际需求情况减设或增设主轴转速范围的离散数、轴向切深范围的离散数和单个刀齿切削周期的离散数。
[0212]
在一个具体实施例中，对于公式(24)中的积分项，分别使用三阶牛顿和埃尔米特插值逼近状态项和时滞项。对于静态力项f(δ)，则采用线性插值逼近。
[0213]
对于公式(24)中的状态项κ(δ)，使用三阶牛顿插值多项式，即κ(δ)由状态向量在
t
k-2
、t
k-1
、tk和t
k+1
时刻的响应κ
k-2
、κ
k-1
、κk和κ
k+1
近似为：
[0214]
κ(δ)＝a1κ
k+1
+b1κk+c1κ
k-1
+d1κ
k-2
ꢀꢀ
(25)
[0215]
式中，插值系数a1、b1、c1和d1分别为
[0216][0217][0218][0219][0220]
κ
k+1
＝a
k+1
x
k+1
,κk＝akxk,κ
k-1
＝a
k-1
x
k-1
,κ
k-2
＝a
k-2
x
k-2
ꢀꢀ
(30)
[0221]
在一个具体实施例中，对于公式(24)中的时滞项α(δ-t)，使用埃尔米特插值多项式进行逼近，具体如下：
[0222]
α(δ-t)＝a2α
k-m
+b2α
k-m+1
+c2α
k-m+2
ꢀꢀ
(31)
[0223]
式中，插值系数a2、b2和c2分别为：
[0224][0225][0226][0227]
α
k-m
＝bkx
k-m
α
k-m+1
＝b
k+1
x
k-m+1
α
k-m+2
＝b
k+2
x
k-m+2
ꢀꢀ
(35)
[0228]
在一个具体实施例中，对于公式(24)中的静态力项f(δ)以相应的边界值做线性逼近，即：
[0229][0230]
将公式(25)-(36)代入公式(24)可得到如下离散动态映射：
[0231][0232]
式中，
[0233][0234][0235][0236]
[0237][0238][0239][0240][0241][0242]
公式(37)整理可得：
[0243][0244]
在一个具体实施例中，使用扩展的状态向量zk＝[x1,x2,
…
,x
m+1
]
t
，yk＝[f1,f2,
…
,f
m+1
]
t
，则公式(47)可以表示为
[0245]
d1zk＝d2z
k-m
+d3ykꢀꢀ
(48)
[0246]
式中，
[0247][0248][0249][0250]
系统在单个时间周期上状态传递关系可以表示为：
[0251]
zm＝φz0+d3ykꢀꢀ
(52)
[0252]
式中，
[0253]
φ＝(d1)-1
d2
ꢀꢀ
(53)
[0254]
基于floquet理论，系统稳定性由传递矩阵φ的特征值的模决定，即：
[0255][0256]
在一个具体实施例中，提取出相应的稳态系数向量y方向系数，预报动态加工误差的具体方法为：
[0257]
对于多周期情况，第i周期和第(i-1)周期的状态传递关系可表示为
[0258]
zi＝φz
i-1
+d3ykꢀꢀ
(55)
[0259]
在稳态切削时，稳态系数可从公式(52)的不动点获得：
[0260]
zk＝z
k-m
＝z
*
ꢀꢀ
(56)
[0261]
将公式(56)代入到公式(55)得：
[0262]z*
＝(d
1-d2)-1
d3ykꢀꢀ
(57)
[0263]
从公式(57)中提取出相应的稳态系数向量y方向系数，预报动态加工误差。
[0264]
另一方面，本发明实施例1还提供了一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，处理器执行程序时实现铣削精加工的稳定性与表面位置误差综合预测方法。
[0265]
综上，在一个具体实施例中，本发明切削半径与进给量，再将作用在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上，并通过累加作用计算作用在铣刀上x和y方向的总切削力，然后在结构动力学框架下构建综合考虑模态耦合效应、过程阻尼和刀具跳动的铣削动力学时滞微分方程组，在任意一个小间隔内对其进行积分计算，并利用三阶牛顿、埃尔米特和线性插值多项式分别近似计算积分式中的状态项、时滞项和静力项，进而构建系统在单个时间周期上的状态传递关系，最后基于floquet理论和不动点理论，同时输出铣削稳定性叶瓣图和表面位置误差预测图。在精加工时，应优先考虑表面质量，在满足表面质量的前提下，再尽可能保证最大加工效率。本发明方法能实现铣削精加工的切削稳定性与表面位置误差综合预测，从而精确预测精加工过程中的表面质量，进而达到提高加工精度和质量的效果，具有很好的准确性及时效性。
[0266]
实施例2
[0267]
当给定相应参数：铣刀刀齿数目n＝3，式中刀具直径d＝25.4mm；刀具磨损带宽度w＝0.08mm，压入系数k
sp
＝1.5
×
10
14
n/m3，库仑力系数μ＝0.3；切向和法向力系数分别为k
t
＝7
×
108n/m2和kn＝0.49
×
108n/m2；刀具偏心参数为偏心距ρ＝7.2μm及偏心角度θ0＝65.09
°
，详见表1。将主轴周期划分为100个小间隔，将由主轴转速与径向切深构成的平面划分为200
×
200网格。
[0268]
表1铣削相关参数
[0269][0270]
将上述步骤和参数通过matlab软件进行编程，画出稳定性叶瓣图和表面位置误差图，通过稳定性图来预测铣削过程中的稳定性和加工过程表面位置误差。图2是按照本发明优选实施例所构建的稳定性叶瓣图，图3是按照本发明的优选实施例同时构建的稳定性叶瓣图和加工表面位置误差图，分别取不同的偏心距和偏心角得到的稳定性叶瓣图如图4和图5所示。
[0271]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0272]
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

技术特征：
1.一种铣削精加工的稳定性与表面位置误差综合预测方法，基于铣削精系统，其特征在于，包括以下步骤：步骤1：获取铣刀信息；步骤2：基于铣刀信息计算任意时刻切削刃角度位置、切削半径、进给量及在x和y方向总切削力；步骤3：构建铣削动力学时滯微分方程组，并转化为空间状态方程；步骤4：将刀齿切削周期t等分成若干个时间间隔，在任意一个时间间隔的时间长度内对所述空间状态方程进行积分计算，确定主轴转速范围及其对应的离散数、轴向切深范围及其对应离散数、单个刀齿切削周期离散数；步骤5：计算所述空间状态方程中的状态项、时滞项和静态力项，构建铣削精加工系统在单个时间周期上的状态传递关系，判定铣削精加工系统的稳定性；步骤6：提取出相应的稳态系数向量y方向系数，预报动态加工误差；步骤7：基于步骤5和步骤6，获得稳定性叶瓣图和表面位置误差预测图。2.根据权利要求1所述的一种铣削精加工的稳定性与表面位置误差综合预测方法，其特征在于，所述步骤2基于铣刀信息计算任意时刻切削刃角度位置、切削半径、进给量及在x和y方向总切削力的具体步骤为：步骤2.1：计算任意时刻切削刃的刀具角度位置；步骤2.2：计算刀具偏心情况下切削刃实际每齿的切削半径和每齿进给量；步骤2.3：计算在刀具切向和径向的切削力，并通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上；步骤2.4：通过笛卡尔坐标系转换将切削力在刀具切向和径向的切削力转换到机床坐标系的x轴和y轴上，根据累加作用在刀齿上的切削力计算作用在铣刀上x和y方向的总切削力。3.根据权利要求1所述的一种铣削精加工的稳定性与表面位置误差综合预测方法，其特征在于，所述步骤3构建铣削动力学时滯微分方程组，并转化为空间状态方程的具体方法为：步骤3.1：构建铣削动力学时滯微分方程组：式中，m为模态质量矩阵，m
xx
、m
yy
为x和y向的质量，m
xy
、m
yx
为x、y向在此坐标下的耦合质量；c为模态阻尼矩阵，c
xx
、c
yy
为x和y向的阻尼，c
xy
、c
yx
为x、y向在此坐标下的耦合阻尼；k为模态刚度矩阵，k
xx
、k
yy
为x和y向的刚度，k
xy
、k
yx
为x、y向在此坐标下的耦合刚度；c
eq
为等效粘性阻尼，c
eq
＝(k
sp
a
p
w2)/(4v)，k
sp
为压入系数，a
p
为轴向切削深度，w为刀具磨损带宽，v为切削速度；c
p
、h为切削力的转化系数矩阵，有：
f0(t)为静态力向量，f
t
为每齿进给量；步骤3.2：令借助柯西变换，转为空间状态形式，得到空间状态方程为：式中，式中，a0、a(t)、b(t)、f(t)由合并同类项得出，a0是x(t)分化出的常系数矩阵，不随时间变化，而a(t)、b(t)、f(t)皆随时间变化。4.根据权利要求1所述的一种铣削精加工的稳定性与表面位置误差综合预测方法，其特征在于，所述步骤4，包括：将刀齿切削周期t等分成若干个时间间隔后，在任意一个时间间隔的时间长度内对所述空间状态方程进行积分计算的具体方法为：定义κ(t)＝a(t)x(t)和α(t-t)＝b(t)x(t-t)，以x(t0)为初始条件，空间状态方程的通解表示为：将刀齿切削周期t等分成为m个时间间隔，则时间步长τ＝t/m；上式在时间kτ≤t≤(k+1)τ，(k＝1,2,3,
…
m)区段上表示为：式中，δ＝ξ-kτ，δ∈[0,τ]，x((k+1)τ)记作x
k+1
。5.根据权利要求1所述的一种铣削精加工的稳定性与表面位置误差综合预测方法，其特征在于，所述步骤5计算所述空间状态方程中的状态项、时滞项和静态力项，构建铣削精加工系统在单个时间周期上的状态传递关系，判定铣削精加工系统的稳定性的具体步骤为：步骤5.1：利用三阶牛顿、埃尔米特和线性插值多项式分别计算所述空间状态方程中的状态项、时滞项和静态力项；步骤5.2：将状态项、时滞项和静态力项带入所述空间状态方程，得到如下离散动态映射：(f
k+1-i)x
k+1
+(t1+f
k
)x
k
+f
k-1
x
k-1
+f
k-2
x
k-2
＝f
k-m
x
k-m
+f
k-m+1
x
k-m+1
+f
k-m+2
x
k-m+2-g
k
f
k-g
k+1
f
k+1
步骤5.3：使用扩展的状态向量zk＝[x1,x2,
…
,x
m+1
]
t
，yk＝[f1,f2,
…
,f
m+1
]
t
，则离散动
态映射表示为：d1z
k
＝d2z
k-m
+d3y
k
式中，式中，式中，步骤5.4：建立铣削精加工系统在单个时间周期上状态传递关系；z
m
＝φz0+d3y
k
式中，φ＝(d1)-1
d2基于floquet理论，铣削精加工系统稳定性由传递矩阵φ的特征值的模决定，具体公式为：6.根据权利要求5所述的一种铣削精加工的稳定性与表面位置误差综合预测方法，其特征在于，所述步骤5.1利用三阶牛顿插值多项式计算所述空间状态方程中的状态项的具体公式为：κ(δ)由状态向量在t
k-2
、t
k-1
、t
k
和t
k+1
时刻的响应κ
k-2
、κ
k-1
、κ
k
和κ
k+1
近似为：κ(δ)＝a1κ
k+1
+b1κ
k
+c1κ
k-1
+d1κ
k-2
式中，κ
k+1
＝a
k+1
x
k+1
,κ
k
＝a
k
x
k
,κ
k-1
＝a
k-1
x
k-1
,κ
k-2
＝a
k-2
x
k-2
，式中t
k-2
＝(k-2)τ，κ((k+1)τ)记作κ
k+1
＝a((k+1)τ)x((k+1)τ)，a((k+1)τ)记作a
k+1
，插值系数a1、b1、c1和d1分别为：
7.根据权利要求5所述的一种铣削精加工的稳定性与表面位置误差综合预测方法，其特征在于，所述步骤5.1利用埃尔米特插值多项式计算所述空间状态方程中的时滞项的具体公式为：α(δ-t)＝a2α
k-m
+b2α
k-m+1
+c2α
k-m+2
式中，α
k-m
＝b
k
x
k-m
α
k-m+1
＝b
k+1
x
k-m+1
α
k-m+2
＝b
k+2
x
k-m+2
，式中α(k-m-t)记作α
k-m
＝b(k-m)x(k-m-t)，b(kτ)记作b
k
，插值系数a2、b2和c2分别为：分别为：分别为：8.根据权利要求5所述的一种铣削精加工的稳定性与表面位置误差综合预测方法，其特征在于，所述步骤5.1利用线性插值多项式计算所述空间状态方程中的静态力项的具体方法为：式中，f(kτ)记作f
k
。9.根据权利要求1所述的一种铣削精加工的稳定性与表面位置误差综合预测方法，其特征在于，所述步骤6：提取出相应的稳态系数向量y方向系数，预报动态加工误差的具体方法为：步骤6.1：对于多周期情况，第i周期和第(i-1)周期的状态传递关系表达式为：z
i
＝φz
i-1
+d3y
k
步骤6.2：在稳态切削时，稳态系数为：z
k
＝z
k-m
＝z
″
步骤6.3：将稳态系数代入到状态传递关系表达式得到：z
*
＝(d
1-d2)-1
d3y
k
；步骤6.4：从步骤6.3的公式中提取出相应的稳态系数向量y方向系数，预报动态加工误差。10.一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1至9任一所述的铣削精加工的稳定性与表面位置误差综合预测方法。

技术总结
本发明公开了铣削精加工的稳定性与表面位置误差综合预测方法及设备，属于机械加工中的铣削精加工技术领域。包括：获取铣刀信息；计算任意时刻切削刃角度位置、切削半径、进给量及在X和Y方向总切削力；构建铣削动力学时滯微分方程组并转化为空间状态方程；对空间状态方程进行积分计算，确定主轴转速范围及其对应的离散数、轴向切深范围及其对应离散数、单个刀齿切削周期离散数；计算空间状态方程中的状态项、时滞项和静态力项，构建铣削精加工系统在单个时间周期上的状态传递关系，判定铣削精加工系统的稳定性；提取出相应的稳态系数向量Y方向系数，预报动态加工误差；获得稳定性叶瓣图和表面位置误差预测图。本发明具有很好的准确性及时效性。确性及时效性。确性及时效性。


技术研发人员：杨雪峰 杨文安 黄超 周芸翔
受保护的技术使用者：南京航空航天大学
技术研发日：2023.06.19
技术公布日：2023/8/31