一种采用传递矩阵缩减技术的铣削过程的颤振预测方法

1.本技术涉及铣削过程的颤振预测技术领域，特别地涉及一种采用传递矩阵缩减技术的铣削过程的颤振预测方法。


背景技术：

2.铣削颤振通常由再生效应引起，是一种不稳定的现象，不利于实现高性能加工。常用的避免颤振的方法之一是使用稳定叶瓣图(sld)选择无颤振参数，它可以描述不同切削深度和主轴组合下铣削过程的稳定域和不稳定域之间的边界速度。铣削过程的控制方程通常被表述为周期性延迟微分方程(dde)。但是用于评估稳定性边界的封闭形式的解决技术是不可用的，除非预设了一些强简化，这对预测精度是不利的。因此，在过去的几十年中，已经提出了各种近似数值方法。
3.tlusty通过一系列时域模拟提出了一种稳定叶瓣的评估方法。campomanes 和altintas提出了一种改进的铣削时域模型来模拟非常小的径向切削宽度下的振动切削条件。尽管这种方法可以准确地解释难以解析建模的非线性效应，但计算量显著增加，所以它们效率低下。
4.为此，频域方法和时域半解析方法被提出。许多研究人员没有解决频域中的 dde，而是专注于解决时域中的dde。除了不同的离散化方法外，同一种离散化方法的收敛速度也会影响计算的准确性和效率。因此，许多工作来通过使用高阶数值方法来获得高收敛速度。高阶方法通常不仅会导致更高的收敛速度，而且会导致更大的计算复杂度。
5.在时域半解析方法中，总是建立以转移矩阵形式表示的离散动态映射，以将一个周期的运动与早一个周期的相应运动联系起来。然后基于floquet理论，通过判断状态转移矩阵的特征值的最大模数是否大于1，可以检测出铣削颤振。这意味着如果可以降低转移矩阵的维数，则可以提高计算效率。本技术提出一种基于5点高斯求积的切削稳定性预测方法，以提高收敛速度。更重要的是，首先提出了一种有效的降低转移矩阵维数的技术，并给出了该技术的基本原理。该技术可以进一步提高所提出方法或同类数值积分方法的计算效率。


技术实现要素：

6.针对上述现有技术中的问题，本技术提出了一种采用传递矩阵缩减技术的铣削过程的颤振预测方法，包括以下步骤：
7.步骤s1、建立再生效应的二自由度铣削系统的运动控制式；
8.步骤s2、将所述运动控制式转化为状态空间形式方程；
9.步骤s3、采用五点gauss-legendre求积法则在对应时间区间内逼近所述状态空间形式方程的解的积分项；
10.步骤s4、通过传递矩阵缩减技术降低传递矩阵的维数；
11.步骤s5、通过floquet理论确定颤振稳定性。
12.在一个实施方式中，所述步骤s1包括建立再生效应的二自由度铣削系统的运动控制式：
[0013][0014]
其中q分别是加速度、速度、位移向量；m、c、k分别是质量、阻尼和刚度矩阵；a
p
和t代表轴向切削深度和时滞量；用于等间距铣刀时，t是单齿切削周期，t＝60/nfω；其中，nf和ω分别代表刀具齿数和主轴转速；kc代表随时间周期变化的铣削力系数矩阵；
[0015]
kc(t)＝kc(t-t)，定义为：
[0016][0017]
其中
[0018][0019][0020][0021]
；
[0022]
式中k
t
和kn代表线性化切向与法向切削力系数；φj(t)表示第j刀齿的角位移：
[0023][0024]
g(φj(t))是一个窗口函数，如果第j个切削刃在切削中则等于1，如果它不在切削中则等于0；
[0025][0026]
其中，φ
st
(t)和φ
ex
(t)分别为切入和切出角。
[0027]
在一个实施方式中，所述步骤s2包括将所述运动控制式转化为状态空间形式方程：
[0028][0029]
其中
[0030][0031]
[0032][0033]
将b(t)(ψ(t)-ψ(t-t))项视为齐次方程的非齐次项，所述状态空间形式方程的解能够表示为
[0034][0035]
其中ψ(t0)表示初始时刻t＝t0的状态值，a是一个常数矩阵。
[0036]
在一个实施方式中，所述步骤s3包括将时间间隔[t0+tf,t0+t]等分为n 个小区间，时间步长记为δt＝(t-tf)/n；相应的采样时间点为ti＝t0+tf+ (i-1)δt，其中i＝1,..n+1和ψ(t1)＝ψ(t0+tf)；由于积分符号下的函数ψ(t)是未知的，所以所述状态空间形式方程的解是一种积分方程；采用五点gauss-legendre求积法则在时间区间[t
i-4
,ti](i＝5,6,...,n+1)内逼近所述状态空间形式方程的解的积分项：
[0037]
；
[0038]
其中五点gauss-legendre求积法则的加权系数ωj＝[0.2369269,0.4786287,0.5688889,0.4786287,0.2369269]， j＝0、1、2、3、4；表示积分时刻，它是通过将高斯正交点从区间[-1,1] 缩放到时间区间[t
i-4
,ti]获得的，如下所示：
[0039][0040]
在区间[-1,1]中的高斯节点ζj＝[-0.9061798,-0.5384693,0,+0.5384693, +0.9061798]；
[0041]
状态项通过边界值ψ(t
i-4+j
)使用二次插值法来近似；由于和位于时间间隔[t
i-4
,t
i-3
]，这两个值用ψ(t
i-4+j
)插值；
[0042][0043]
相似的，和利用ψ(t
i-2
)插值，
[0044][0045][0046]
对延时状态项进行插值：
[0047]
[0048][0049]
得到
[0050]
ψi＝e
4δta
ψ
i-4 +2δya
p
w0e0b
i-4
[a0(ψ
i-4-ψ
i-4
，t)+b0(ψ
i-3-ψ
i-3，t
)+c0(ψ
i-2-ψ
i-2，t
)] +2δya
p
w1e1b
i-3
[a1(ψ
i-4-ψ
i-4，t
)+b1(ψ
i-3-ψ
i-3，t
)+c1(ψ
i-2-ψ
i-2，t
)] +2δya
p
w2e2b
i-2
(ψ
i-2-ψ
i-2，t
) +2δya
p
w3e3b
i-1
[a3(ψ
i-2-ψ
i-2，t
)+b3(ψ
i-1-ψ
i-1，t
)+c3(ψ
i-ψ
i，t
)] +2δya
p
w4e4bi[a4(ψ
i-2-ψ
i-2，t
)+b4(ψ
i-1-ψ
i-1，t
)+c4(ψ
i-ψ
i，t
)]；
[0051]
将符号简化以使以下公式推导更简洁，即：
[0052][0053][0054]
ψ
i-j
＝ψ(t
i-j
)(j＝0,1,2,3,4)
[0055]
ψ
i-j,t
＝ψ(t
i-j-t)(j＝0,1,2,3,4)
[0056][0057]
并且
[0058][0059]
当前齿周期的状态能够链接到前一个齿周期的状态，用于恒螺距和螺旋刀具；中间状态(i＝5,6,...,n+1)表示为以下一般形式：
[0060][0061]
其中d
1i
＝[0 0 0 0 i]，d
2i
＝[e
4δta
0 0 0],d
3i
表示为
[0062][0063]
在时间区间[t
i-1
,ti](i＝2，3，4)内，采用两点gauss-legendre求积法则逼近所述状态空间形式方程的解的积分项，并采用线性插值法对状态项进行插值；中间状态(i＝2、3、4)用一般形式表示如下：
[0064][0065]
其中d
1i
＝[0i]，d
2i
＝[e
4δta
0]，d
3i
表示为：
[0066][0067]
在两点gauss-legendre求积法则中和都等于1；
[0068][0069][0070]
ψ
i-j
＝ψ(t
i-j
)(j＝0,1)
[0071]
ψ
i-j,t
＝ψ(t
i-j-t)(j＝0,1)
[0072][0073]
其中
[0074]
考虑一个时期内的所有状态，构造以下离散动态映射：
[0075][0076]
其中i是单位矩阵；
[0077][0078]
[0079][0080]
那么大的传递矩阵可以写成
[0081]
φ＝l-1
n；
[0082]
其中l＝i-c-δta
p
p；
[0083]
并且n＝-δta
p
p+g。
[0084]
在一个实施方式中，所述步骤s4包括：
[0085]
步骤s41：通过方程m＝i-c-δta
p
p和n＝-δta
p
p+g计算(4n+4)
×ꢀ
(4n+4)矩阵m和n；
[0086]
步骤s42：去除矩阵n的第(4i-1)和(4i)列，其中i＝1,2,...,n，得到 (4n+4)
×
(2n+4)矩阵
[0087]
步骤s43：进行左矩阵除法运算得到矩阵
[0088]
步骤s44：删除h矩阵的第(4i-1)和(4i)行，其中i＝1,2,...,n，得到 (2n+4)
×
(2n+4)矩阵
[0089]
在一个实施方式中，所述步骤s5包括根据floquet理论求解矩阵的特征值以评估铣削稳定性。
[0090]
上述技术特征可以各种适合的方式组合或由等效的技术特征来替代，只要能够达到本技术的目的。
[0091]
本技术提供的一种采用传递矩阵缩减技术的铣削过程的颤振预测方法，与现有技术相比，至少具备有以下有益效果：
[0092]
1、本技术具有较高的计算精度和效率。利用五点gauss-legendre求积法则来近似以积分方程形式表示的铣削系统运动方程。然后建立由转移矩阵表示的当前和前一个齿周期状态之间的离散动态映射，本技术所提出的方法比1st sdm、 2nd fdm、la和2pgim具有更高的收敛速度。由于这个优势，允许减小离散化参数。由于仅逼近强迫振动，因此所提出方法的离散化参数可以进一步降低。
[0093]
2、本技术的二自由度铣削系统的传递矩阵的维数通过矩阵缩减技术从(4n+4)
ꢀ×
(4n+4)缩减到(2n+4)
×
(2n+4)。当离散化参数或动力学方程的维数很大时，该技术对计算效率非常重要。此外，该技术适用于其他基于数值积分的颤振预测方法，以进一步减少计算时间。将转移矩阵的维数减少了近一半。
附图说明
[0094]
在下文中将基于实施例并参考附图来对本技术进行更详细的描述。其中：
[0095]
图1表示铣削动力学模型；
[0096]
图2为降低转移矩阵维数的技术说明；
[0097]
图3表示一阶半离散化方法(1st sdm)、二阶全离散化方法(2nd fdm)、加速度线性逼近法(la)、两点高斯积分法(2pgim)和所提出的方法的不同离散化参数n的特征值收敛情况；
[0098]
图4表示本技术所提出的方法与la花费的计算时间对比；
[0099]
图5表示所提出的方法预测的sld；
[0100]
图6表示将所提出方法与1st sdm、2nd fdm、la和2pgim在离散化参数 n设置为20，400
×
200大小的参数网格上预测的sld与参考sld高度一致；
[0101]
图7表示实验装置；
[0102]
图8表示将所提方法得到的稳定性边界与实验结果进行比较；
[0103]
图9表示矩阵缩减技术计算效率。
具体实施方式
[0104]
下面将结合附图对本技术作进一步说明。如图1所示，前一个齿周期留下的波浪形表面光洁度在随后的齿周期中被去除。切屑厚度的再生效应是铣削过程中产生颤振的主要原因。
[0105]
本技术提供了一种采用传递矩阵缩减技术的铣削过程的颤振预测方法，包括下述步骤：
[0106]
(步骤1)：考虑再生效应的二自由度(dof)铣削系统的运动由下式控制：
[0107][0108]
其中q分别是加速度，速度，位移向量。m,c,k分别是质量，阻尼和刚度矩阵。a
p
和t代表轴向切削深度和时滞量。用于等间距铣刀时，t是单齿切削周期，t＝60/nfω其中，nf和ω分别代表刀具齿数和主轴转速。kc代表随时间周期变化的铣削力系数矩阵。
[0109]
(步骤2)：由步骤1kc(t)＝kc(t-t)定义为：
[0110][0111]
其中
[0112][0113][0114][0115][0116]
式中k
t
和kn代表线性化切向与法向切削力系数。φj(t)表示图1所示第j刀齿的角位移：
[0117][0118]
(步骤3)：用窗口函数g(φj(t))判断切削刃是否在切削中，在则等于1，如果它不在切削中则等于0
[0119][0120]
其中，φ
st
(t)和φ
ex
(t)分别为切入和切出角。
[0121]
(步骤4)：步骤(1)中的运动方程可以转化为状态空间形式
[0122][0123]
其中
[0124][0125][0126][0127]
(步骤5)：将b(t)(ψ(t)-ψ(t-t))项视为齐次方程的非齐次项，步骤(4)的解可以表示为
[0128][0129]
其中ψ(t0)表示初始时刻t＝t0的状态值。
[0130]
(步骤6)：当铣刀离开工件时，铣削系统经历自由振动过程。因此，b(t)变为零矩阵。在自由振动持续时间结束时，即t＝t0+tf时，状态变为：
[0131]
(步骤7)：当铣刀去除工件材料时，铣削系统经历强迫振动过程。步骤(4) 是一种积分方程。我们采用五点gauss-legendre求积法则在时间区间[t
i-4
, ti](i＝5,6,...,n+1)内逼近步骤(4)方程的积分项：
[0132][0133]
其中五点gauss-legendre求积法则的加权系数ωj＝[0.2369269,0.4786287,0.5688889,0.4786287, 0.2369269](j＝0,1,2,3,4)。
[0134]
(步骤8)：表示积分时刻它是通过将高斯正交点从区间[-1,1]缩放到时间区间[t
i-4
,ti]获得的，如下所示：
[0135][0136]
在区间[-1,1]中的高斯节点ζj＝[[-0.9061798,-0.5384693,0,+0.5384693, +0.9061798]。
[0137]
(步骤9)：状态项可以通过边界值ψ(t
i-4+j
)使用二次插值法来近似。由于和位于时间间隔[t
i-4
,t
i-3
]，这两个值用ψ(t
i-4+j
)插值。和利用ψ(t
i-2
)插值，延时状态项可以用类似的方法进行插值。
[0138][0139]
[0140][0141][0142]
(步骤10)：将步骤(9)的四个公式代入步骤(7)的公式中得到
[0143]
ψi＝e
4δta
ψ
i-4 +2δya
p
w0e0b
i-4
[a0(ψ
i-4-ψ
i-4，t
)+b0(ψ
i-3-ψ
i-3，t
)+c0(ψ
i-2-ψ
i-2，t
)] +2δya
p
w1e1b
i-3
[a1(ψ
i-4-ψ
i-4，t
)+b1(ψ
i-3-ψ
i-3，t
)+c1(ψ
i-2-ψ
i-2，t
)] +2δya
p
w2e2b
i-2
(ψ
i-2-ψ
i-2，t
) +2δya
p
w3e3b
i-1
[a3(ψ
i-2-ψ
i-2，t
)+b3(ψ
i-1-ψ
i-1，t
)+c3(ψ
i-ψ
i，t
)] +2δya
p
w4e4bi[a4(ψ
i-2-ψ
i-2，t
)+b4(ψ
i-1-ψ
i-1，t
)+c4(ψ
i-ψ
i，t
)]
[0144]
可以用符号化简使得公式推导更简洁即
[0145][0146][0147]
ψ
i-j
＝ψ(t
i-j
)(j＝0,1,2,3,4)
[0148]
ψ
i-j,t
＝ψ(t
i-j-t)(j＝0,1,2,3,4)
[0149][0150]
并且
[0151][0152]
(步骤11)基于步骤10中的公式，在使用恒螺距的螺旋刀具当前齿周期的状态可以链接到前一个齿周期的状态。中间状态(i＝5,6,...,n+1)可以表示为以下一般形式
[0153][0154]
其中d
1i
＝[0 0 0 0i]，d
2i
＝[e
4δta
0 0 0],d
3i
表示为
[0155][0156]
(步骤12)在时间区间[t
i-1
,ti](i＝2，3，4)内，采用两点gauss
‑ꢀ
legendre求积法则逼近(步骤5)的积分项，并采用线性插值法对状态项进行插值。经过上述类似的推导，中间状态(i＝2、3、4)可以用一般形式表示如下：
[0157][0158]
其中d
1i
＝[0i]，d
2i
＝[e
4δta
0]，d
3i
表示为
[0159][0160]
(步骤13)在两点gauss-legendre求积法则中和都等于1。
[0161][0162][0163]
ψ
i-j
＝ψ(t
i-j
)(j＝0,1)
[0164]
ψ
i-j,t
＝ψ(t
i-j-t)(j＝0,1)
[0165][0166]
其中
[0167]
(步骤14)考虑一个周期的所有状态，即步骤6步骤11步骤12构造离散动态映射
[0168][0169]
其中i是单位矩阵
[0170]
[0171][0172]
[0173]
(步骤15)那么大的转移矩阵可以写成
[0174]
φ＝l-1
n；
[0175]
其中l＝i-c-δta
p
p；
[0176]
并且n＝-δta
p
p+g。
[0177]
(步骤16)如图2所示，通过步骤15计算(4n+4)
×
(4n+4)矩阵l和n。
[0178]
(步骤17)去除矩阵n的第(4i-1)和(4i)列，其中i＝1,2,...,n，得到 (4n+4)
×
(2n+4)矩阵
[0179]
(步骤18)进行左矩阵除法运算得到矩阵
[0180]
(步骤19)删除h矩阵的第(4i-1)和(4i)行，其中i＝1,2,...,n，得到 (2n+4)
×
(2n+4)矩阵
[0181]
(步骤20)根据floquet理论求解矩阵的特征值以评估铣削稳定性。
[0182]
在一个实施例中，采用双刃铣刀。固有频率f
t
＝ωt/(2π)相对阻尼ξ
t
和模态质量m
t
分别为922hz、0.011和0.03993kg。切削力系数k
t
和kn分别为 6
×
108n/m2和2
×
108n/m2。顺铣切削，三种不同的轴向切削深度a
p
＝0.2,0.5, 1mm,径向浸入比ae/d＝100％，主轴转速ω＝5000rpm。
[0183]
首先，将所提方法的收敛速度与现有方法进行比较。与1st sdm、2nd fdm、加速度线性逼近法(la)和2pgim的单自由度系统相比，所提出方法的临界特征值的误差如图3(a)—3(c)所示。所提出的方法的收敛速度远高于其他方法。该方法的临界特征值误差非常小。该方法有一个更明显的优势，即当离散化次数选择在 30—45之间时，临界特征值的误差可以达到一个满意的值。
[0184]
其次，将所提方法的计算时间与现有方法进行了比较。使用1st sdm、2nd fdm、la、2pgim和所提出的方法在400
×
200大小的主轴速度ω和切削深度a
p
参数网格上计算sld。为了更清楚地展示所提出的方法与其他方法之间的差异，所提出的单自由度铣削稳定性方法的matlab代码采用了分离自由振动和强迫振动以及减小传递矩阵维数的方法。
[0185]
从图3(a)—3(c)、图4(a)—图4(f)可以看出，所提出的方法的收敛速度远高于 1st sdm、2nd fdm、la、2pgim提出的其他方法。与2pgim相比，该方法的计算时间可以减少近20％。所提出的方法花费的计算时间不到2nd fdm的35％。该方法预测的sld与参考sld非常吻合，计算时间分别为5.75秒和11.03秒。图4(a)—图4(f)表明所提出的方法与la花费的计算时间大致相同，这两种方法比其他方法更有效。与2pgim相比，该方法的计算时间可以减少近20％。所提出的方法花费的计算时间不到2nd fdm的35％。更重要的是，该方法比其他方法具有更好的准确性，因为该方法的收敛速度更高。允许显著减小离散化参数n。因此提高了计算效率。图5(a)—5(b)表示所提出的方法预测的sld，具有较小的离散化参数n，与参考sld非常吻合，表明所提出的方法是高效和准确的。
[0186]
在一个实施例中，系统参数的选择与单自由度铣削模型案例中使用的相同，并且假设它们在进给方向和法向方向上相等。离散化参数n设置为20。图6(a)— 6(d)将所提出方法与1st sdm、2nd fdm、la和2pgim在离散化参数n设置为 20，400
×
200大小的参数网格上预测的sld与参考sld高度一致。同时，所提出的方法消耗的时间最少，效率最高。表明所提出的方法非常适用于双自由度铣削过程的稳定性预测。在6(a)—6(d)中比较了通过这些
方法在400
×
200大小的参数网格上预测的sld，用于径向浸入比ae/d＝0.05、0.1、0.5和1的顺铣。可以看出，该方法的结果与参考sld高度一致。同时，所提出的方法消耗的时间最少，效率最高。这些例子表明，所提出的方法非常适用于双自由度铣削过程的稳定性预测。
[0187]
该方法在计算时间上的优势可以总结如下：
[0188]
通过矩阵缩减技术，2自由度铣削系统的传递矩阵的维数从(4n+4)
×
(4n +4)减少到(2n+4)
×
(2n+4)。对于单自由度铣削系统，传递矩阵的维数为(n+2)
×
(n+2)，小于la(n+3)
×
(n+3)和2pgim(2n+2)
×
(2n+ 2).因此，在求解矩阵特征值的过程中将减少计算时间。
[0189]
c、p、g与主轴转速有关，与轴向切深无关，计算方便。
[0190]
此外，避免了floquet转移矩阵φ的递归构造，从而节省了大量时间
[0191]
通过研究单自由度和双自由度铣削模型的示例，讨论了所提出方法的效率和准确性。提出的和现有方法的所有代码均在matlab 2019b中编写，并在个人计算机[intel(r) xeon(r) w-2123 cpu,3.60ghz,32gb]上运行。
[0192]
为了验证所提出的方法，在主轴转速范围为3500到4600rpm的条件下进行了一组逆铣试验。实验装置如图7所示，工件材料为铝7050，刀具为四齿立铣刀，使用测功机kistler 9255b和麦克风收集切割力和声音。使用测力仪kistler9255b和麦克风收集切割力和声音。切割过程是否稳定可以从力和声音测量的频谱中确定。工件材料为铝7050，刀具为四齿立铣刀，直径d＝12mm。刀具模态参数通过冲击试验获得，列于表1。径向切深选择为1.2mm，即径向浸入比 ae/d＝0.1。因此，径向和切向铣削力系数分别为kn＝355.0n/m2和k
t
＝ 754.8n/m2。sld是通过在上述铣削参数下使用所提出的方法确定的。实验结果和sld如图8所示，表明铣削系统的预测稳定性极限与实验结果吻合良好，这意味着所提出的方法可以准确地预测稳定性。图9表示矩阵缩减技术可以减少30％以上的传递矩阵求解特征值的计算时间。并且，传递矩阵缩减技术还可以提高其他基于数值积分的颤振预测方法的计算效率。
[0193]
虽然在本文中参照了特定的实施方式来描述本技术，但是应该理解的是，这些实施例仅仅是本技术的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本技术的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其他所述实施例中。