[发明专利]五轴数控机床加工中球头铣刀颤振稳定域叶瓣图建模方法

权利要求书

1.一种五轴数控机床加工中基于精细积分的球头铣刀颤振稳定域叶瓣图建模方法，其特征在于以下步骤：

步骤一，建立球头铣刀刀具-工件动力学方程，表示如下：

其中，m_tx，ξ_x，ω_nx分别是为刀具系统x方向的模态质量，阻尼系数，固有频率；m_ty，ξ_y，ω_ny分别是为刀具系统y方向的模态质量，阻尼系数，固有频率；F_tx(t)和F_ty(t)分别为在x，y方向上作用在铣刀上的动态切削力；

步骤二，求解球头铣刀刀齿上的动态切削力F_tx(t)和F_ty(t)，具体为：

其中，

步骤三、五轴数控机床平面加工过程中球头铣刀与工件的接触区域半解析建模

3.1)进行刀具路径规划，设置加工参数

球头铣刀在五轴数控机床上进行平面铣削，球头铣刀球头半径为R，前倾角为α度，侧倾角为β度，相邻切削刀轨间距为L_xl，轴向切削深度为L_jg；

3.2)前倾角为α度、侧倾角为β度时球头铣刀与工件接触区域边界由a，b，c，d四条线组成，其中，a号线为铣刀球头与工件加工表面的交线，b号线为铣刀球头与工件过渡表面的交线，c号线为本次走刀与上一次走刀在已加工表面留下的加工残余而形成的线，d号线为铣刀球头与上一次走刀在工件留下的加工痕迹的交线；另外，e号线为上一次走刀在工件加工表面留下的加工痕迹，根据球头铣刀特点，将接触区域边界求解问题转化为这些边界在垂直于刀具轴线平面的投影方程的求解问题；

3.3)由刀具前倾角和侧倾角为零度时的三维直角坐标系X₀-Y₀-Z₀推导出刀具前倾角为α、侧倾角为零度时的三维直角坐标系X-Y-Z，具体为：

建立刀具前倾角和侧倾角都为零度时的三维直角坐标系X₀-Y₀-Z₀，其中刀具顶点为原点O₀，刀轴线为Z₀轴，X₀轴与刀具进给方向相同；下面将由三维直角坐标X₀-Y₀-Z₀推导出刀具前倾角为α，侧倾角为零度时的三维直角坐标系X-Y-Z，由于侧倾角都为零度，故只需要在X₀-Z₀平面讨论即可；

在X₀-Z₀平面内，将Z₀轴以O₀为原点，顺时针倾斜α度，从Y₀负方向向原点看，α值为正数，顺时针；α值为负数，逆时针，此处以α为正进行阐述，倾斜之后的轴线就是刀具前倾角为α度，侧倾角为零度时的刀具轴线，直线方程z₀＝tan(90-α)×x₀与直线z₀＝R的交点cir_0，即为刀具前倾角为α度，侧倾角为零度时的刀具球心，经计算cir_0的坐标为x_{cir_o}＝R/tan(90-α)，z_{cir_0}＝R；

以cir_0为原点，R为半径建立圆的方程(x₀-x_{cir_0})²+(z₀-y_{cir_0})²＝R²，此圆方程与X₀轴的切点即为前倾角为α度，侧倾角为零度时的球头铣刀的刀触点dc_0，经计算dc_0的坐标x_{dc_0}＝R/tan(90-α)，z_{dc_0}＝0；

(x₀-x_{cir_0})²+(z₀-y_{cir_0})²＝R²方程与z₀＝tan(90-α)×x₀方程两个交点中y₀值较小的点即为前倾角为α度，侧倾角为0度时的球头铣刀的刀位点dw_0；经计算x_{dw_0}＝R/tan(90-α)-R×sinα，z_{dw_0}＝R-R×cosα；

以dc_0为原点，建立三维直角坐标系X₁-Y₁-Z₁，其中，X₁-Z₁平面与X₀-Z₀平面相重合，X₁与X₀相互重合且方向相同，Y₁与Y₀平行且方向相同；

以dw_0为原点，以方程z₀＝tan(90-α)×x₀为Z轴建立三维直角坐标系X-Y-Z，其中，X-Z平面与X₀-Z₀平面相重合，Y与Y₀平行且方向相同，Z轴正方向远离工件，在X-Y-Z坐标系下，铣刀球头轮廓方程为x²+y²+(z-R)²＝R²；

3.4)通过X-Y-Z三维坐标系，建立刀具前倾角为α、侧倾角为β度时的X_m-Y_m-Z_m坐标系，获取X_m-Y_m平面上点(x_m,y_m)与X-Y平面上点(x,y)的关系，具体为：

将X-Y-Z坐标系中的X-Y平面沿着Z轴正方向平移R距离，使坐标系的原点O与铣刀球心重合，获得新的坐标系X_p-Y_p-Z_p；

由X_p-Y_p-Z_p坐标系获得方式可知，X_p-Y_p平面中的点(x_p,y_p)与X-Y平面中的点(x,y)的关系为:

x_p＝x，y_p＝y (3)

以X_p轴为轴线，将Y_p-Z_p平面逆时针旋转β度，X_p轴正半轴向原点看，β值为正数，顺时针；β值为负数，逆时针，以下以β为负值进行阐述，旋转之后，Y_p轴变为Y_f轴，Z_p轴变为Z_f轴，形成一个新的坐标系X_f-Y_f-Z_f，其中X_f与X_p轴重合且方向相同；

由X_f-Y_f-Z_f坐标系获得方式可知，X_p-Y_p-Z_p坐标系中点(x_p,y_p,z_p)与X_f-Y_f-Z_f坐标系中的点(x_f,y_f,z_f)关系为：

令z_f＝z_p＝0，则可得：

x_f＝x_p，y_f＝y_pcosβ (4)

将坐标系X_f-Y_f-Z_f中的X_f-Y_f平面沿着Z_f负方向平移R距离，得到X_m-Y_m-Z_m坐标系；

由X_m-Y_m-Z_m坐标系获得方式，可知X_m-Y_m平面中的点(x_m,y_m)与X_f-Y_f平面中的点(x_f,y_f)的关系为：

x_m＝x_f，y_m＝y_f (5)

由公式(3)，(4)和(5)可得X-Y平面上点与X_m-Y_m平面上点的关系：

x_m＝x，y_m＝ycosβ (6)

由式(6)可知，X_m-Y_m平面上点与X-Y平面上点的关系，因此，a，b，c，d四条边界线在X_m-Y_m坐标系下的投影方程的求解问题转化为这四条边界线在X-Y坐标系下的投影方程的求解问题；

3.5)确定a，b，c，d四条线在X-Y坐标系下的投影方程

3.5.1)a号线在X-Y坐标系下投影方程求解过程为：

a号线为铣刀球头与工件加工表面的交线；工件加工表面在X-Y-Z坐标系下的方程需要通过X₀-Y₀-Z₀坐标系与X-Y-Z坐标系的关系获得；

在X₀-Z₀坐标系下，在工件加工表面任取A(x_A0,z_A0)、B(x_B0,z_B0)两点，由于工件加工面为平面，故z_A0＝z_B0＝L_jg；通过X–Z坐标系与X₀-Z₀坐标系的关系，可获得A、B两点在X–Z坐标系的值(x_A,z_A)，(x_B,z_B)，具体获取方式如下：

获取A、B在X–Z坐标系下坐标之后，计算出经过A、B两点在X–Z坐标系下的直线方程z＝kx+b，其中此方程在X–Z坐标系下为直线方程，在X-Y-Z下可表示为工件加工表面方程；

在X-Y-Z下，将工件加工表面方程z＝kx+b与铣刀球头方程x²+y²+(z-R)²＝R²方程联立，即可得到a号线在X-Y二维坐标系下投影方程x²+y²+(kx+b-R)²＝R²；

3.5.2)b号线在X-Y坐标系下投影方程为

3.5.3)c号线在X-Y坐标系下的投影方程为y＝-L_xl/2；

3.5.4)d号线求解过程在X-Y坐标系下的投影方程为六次多项式方程；

3.6)确定a，b，c，d在X_m-Y_m坐标系下的投影方程

利用x＝x_m，y＝y_m/cosβ分别替换掉a、b、c、d在X-Y坐标系下投影方程中的x、y，获得的新的方程为a、b、c、d在X_m-Y_m坐标系下的投影方程，它们共同围成的区域即为前倾角为α、侧倾角为β时球头铣刀与工件接触区域在X_m-Y_m坐标系下的投影；

步骤四、精细积分法对刀具-工件动力学方程时域数值求解

4.1)精细积分法对刀具-工件动力学方程进行时域数值求解

由公式(2)将公式(1)表示为如下形式：

将式(7)表示为如下所示的哈密顿系统：

其中，

将式(8)中A(t)v(t)-A(t)v(t-T)用f(t)来表示，则对于非齐次方程(8)，由常微分理论可知，一般解为：

将时滞周期T均分为m份，即T＝m·τ，在[t_p,t_p+1]中，将f(t)表示为如下形式：

f(t)＝r₀+r₁(t-t_p) (10)

其中，r₀＝f(t_p)＝A(t_p)v(t_p)-A(t_p)v(t_p-m·τ)；

由式(9)和式(10)可将v(t_p+1)表示为：

v(t_p+1)＝T·[v(t_p)+A₀^-1(r₀+A₀^-1r₁)]-A₀^-1(r₀+A₀^-1r₁+r₁·τ) (11)

其中，将τ均分为Λ＝2²⁰份，则

由于足够小，故用Taylor展开式去近似表示为

由式(12)和式(13)可知，

T_a＝A₀Δt+(A₀Δt)²/2！+(A₀Δt)³/3！+(A₀Δt)⁴/4！ (15)

将r₀和r₁进一步分别表示为

r₀＝A_pv_p-A_pv_p-m (16)

将(16)和(17)带入到式(11)中，可得：

其中，MM＝TA₀^-1-A₀^-1，NN＝TA₀^-2-A₀^-2-A₀^-1τ；

若(I-NN/τ·A_p+1)可逆，则式(18)可表示为：

其中，

4.2)确定a_xx,p、a_xy,p、a_yx,p、a_yy,p和a_xx,p+1、a_xy,p+1、a_yx,p+1、a_yy,p+1

以N_f＝2阐述a_xx,p、a_xy,p、a_yx,p、a_yy,p的确定方法；在t_p时刻，将球头铣刀全部切削刃投影到X_m-Y_m坐标系下，其中，n为刀具转速，j从1到N_f，k从0到π/2；由铣刀刀齿数和构建出的刀具-工件的接触区域可知在时滞周期内的任意时刻，只有一条切削刃参与切削，设此切削刃为第一条切削刃，该切削刃在单齿切削周期内到达与b号线相切位置、b投影线与a投影线交点位置、与d号线相切位置、a投影线与d投影线交点位置、d投影线与c投影线交点位置、c投影线与b投影线交点位置、时滞周期切削结束的时刻分别记为t₁、t₂、t₃、t₄、t₅、t₆、t₇；

当t_p＝0时，第一条切削刃没有参与切削；当t_p在0-t₁时间内，第一条切削刃没有参与切削，此时k_max,1和k_min,1值都为零；当t_p在t₁-t₂时间内，第一条切削刃与b号线投影线相交，两个交点所对应的轴向角中的最大值即为t_p时刻该切削刃在公式(20)中所对应的k_max，1，最小值为t_p时刻该切削刃在公式(20)中所对应的k_min,1；同理，能够确定t_p分别在t₂-t₃，t₃-t₄，t₄-t₅，t₅-t₆时间段内，第一条切削刃所对应k_max,1和k_min,1值；当t_p在t₆-t₇时间内，第一条切削刃没有参与切削，此时k_max,1和k_min,1值都为零；

确定出各个切削刃在t_p参与切削的最大轴向角和最小轴向角之后，通过式(20)得到a_xx,p、a_xy,p、a_yx,p、a_yy,p值；

通过以上步骤得到t_p+1时刻时a_xx,p+1、a_xy,p+1、a_yx,p+1、a_yy,p+1值；

步骤五、叶瓣图构建

5.1)切削稳定性判定方法

建立系数矩阵C_p，该矩阵满足离散映射：

v_p+1＝C_pv_p (21)

其中，v_p是个(2m+4)维的向量，表示为：

矩阵C_p为(2m+4)维矩阵，表示为：

矩阵PK为4×4矩阵等于式(19)中的(I-NN/τ·A_p+1)^-1(T+MM·A_p-NN/τ·A_p)，矩阵RK1为4×2矩阵，等于式(19)中的-(I-NN/τ·A_p+1)^-1NN/τ·A_p+1的前两列，RK2为4×2矩阵，等于式(19)中的(I-NN/τ·A_p+1)^-1(NN/τ·A_p-MM·A_p)的前两列；

通过使用一系列离散C_p，p＝0,1,2…,m-1，构建时滞周期T内的过渡矩阵Φ，亦即：

v_p＝Φv₀ (24)

式中，Φ定义为：Φ＝C_m-1C_m-2…C₁C₀；

由Floquet理论可知，传递函数Φ所有特征值模的最大值小于1，等于1和大于1，分别表示在该刀具转速n和轴向切深L_jg下，切削处于稳定状态，临界状态和不稳定状态；

5.2)叶瓣图构建

主轴转速不变，不断改变刀具轴向切深，判断其所对应的切削状态，获得在[0,L_jg]范围内的临界轴向切深；改变轴向切深后的接触区域边界投影中只有a号线投影方程变化，b号线和c号线投影方程仍与L_jg时相同；改变主轴转速，获得不同主轴转速在[0,L_jg]范围内所对应的临界轴向切削深度；最终，构建出临界轴向切削深度随主轴转速变化的函数关系，即为颤振稳定域叶瓣图。