1.一种五轴数控机床加工中基于精细积分的球头铣刀颤振稳定域叶瓣图建模方法，其特征在于以下步骤：
步骤一，建立球头铣刀刀具-工件动力学方程，表示如下：
其中，mtx，ξx，ωnx分别是为刀具系统x方向的模态质量，阻尼系数，固有频率；mty，ξy，ωny分别是为刀具系统y方向的模态质量，阻尼系数，固有频率；Ftx(t)和Fty(t)分别为在x，y方向上作用在铣刀上的动态切削力；
步骤二，求解球头铣刀刀齿上的动态切削力Ftx(t)和Fty(t)，具体为：
其中，
步骤三、五轴数控机床平面加工过程中球头铣刀与工件的接触区域半解析建模
3.1)进行刀具路径规划，设置加工参数
球头铣刀在五轴数控机床上进行平面铣削，球头铣刀球头半径为R，前倾角为α度，侧倾角为β度，相邻切削刀轨间距为L_xl，轴向切削深度为L_jg；
3.2)前倾角为α度、侧倾角为β度时球头铣刀与工件接触区域边界由a，b，c，d四条线组成，其中，a号线为铣刀球头与工件加工表面的交线，b号线为铣刀球头与工件过渡表面的交线，c号线为本次走刀与上一次走刀在已加工表面留下的加工残余而形成的线，d号线为铣刀球头与上一次走刀在工件留下的加工痕迹的交线；另外，e号线为上一次走刀在工件加工表面留下的加工痕迹，根据球头铣刀特点，将接触区域边界求解问题转化为这些边界在垂直于刀具轴线平面的投影方程的求解问题；
3.3)由刀具前倾角和侧倾角为零度时的三维直角坐标系X0-Y0-Z0推导出刀具前倾角为α、侧倾角为零度时的三维直角坐标系X-Y-Z，具体为：
建立刀具前倾角和侧倾角都为零度时的三维直角坐标系X0-Y0-Z0，其中刀具顶点为原点O0，刀轴线为Z0轴，X0轴与刀具进给方向相同；下面将由三维直角坐标X0-Y0-Z0推导出刀具前倾角为α，侧倾角为零度时的三维直角坐标系X-Y-Z，由于侧倾角都为零度，故只需要在X0-Z0平面讨论即可；
在X0-Z0平面内，将Z0轴以O0为原点，顺时针倾斜α度，从Y0负方向向原点看，α值为正数，顺时针；α值为负数，逆时针，此处以α为正进行阐述，倾斜之后的轴线就是刀具前倾角为α度，侧倾角为零度时的刀具轴线，直线方程z0＝tan(90-α)×x0与直线z0＝R的交点cir_0，即为刀具前倾角为α度，侧倾角为零度时的刀具球心，经计算cir_0的坐标为xcir_o＝R/tan(90-α)，zcir_0＝R；
以cir_0为原点，R为半径建立圆的方程(x0-xcir_0)2+(z0-ycir_0)2＝R2，此圆方程与X0轴的切点即为前倾角为α度，侧倾角为零度时的球头铣刀的刀触点dc_0，经计算dc_0的坐标xdc_0＝R/tan(90-α)，zdc_0＝0；
(x0-xcir_0)2+(z0-ycir_0)2＝R2方程与z0＝tan(90-α)×x0方程两个交点中y0值较小的点即为前倾角为α度，侧倾角为0度时的球头铣刀的刀位点dw_0；经计算xdw_0＝R/tan(90-α)-R×sinα，zdw_0＝R-R×cosα；
以dc_0为原点，建立三维直角坐标系X1-Y1-Z1，其中，X1-Z1平面与X0-Z0平面相重合，X1与X0相互重合且方向相同，Y1与Y0平行且方向相同；
以dw_0为原点，以方程z0＝tan(90-α)×x0为Z轴建立三维直角坐标系X-Y-Z，其中，X-Z平面与X0-Z0平面相重合，Y与Y0平行且方向相同，Z轴正方向远离工件，在X-Y-Z坐标系下，铣刀球头轮廓方程为x2+y2+(z-R)2＝R2；
3.4)通过X-Y-Z三维坐标系，建立刀具前倾角为α、侧倾角为β度时的Xm-Ym-Zm坐标系，获取Xm-Ym平面上点(xm,ym)与X-Y平面上点(x,y)的关系，具体为：
将X-Y-Z坐标系中的X-Y平面沿着Z轴正方向平移R距离，使坐标系的原点O与铣刀球心重合，获得新的坐标系Xp-Yp-Zp；
由Xp-Yp-Zp坐标系获得方式可知，Xp-Yp平面中的点(xp,yp)与X-Y平面中的点(x,y)的关系为:
xp＝x，yp＝y               (3)
以Xp轴为轴线，将Yp-Zp平面逆时针旋转β度，Xp轴正半轴向原点看，β值为正数，顺时针；β值为负数，逆时针，以下以β为负值进行阐述，旋转之后，Yp轴变为Yf轴，Zp轴变为Zf轴，形成一个新的坐标系Xf-Yf-Zf，其中Xf与Xp轴重合且方向相同；
由Xf-Yf-Zf坐标系获得方式可知，Xp-Yp-Zp坐标系中点(xp,yp,zp)与Xf-Yf-Zf坐标系中的点(xf,yf,zf)关系为：
令zf＝zp＝0，则可得：
xf＝xp，yf＝ypcosβ           (4)
将坐标系Xf-Yf-Zf中的Xf-Yf平面沿着Zf负方向平移R距离，得到Xm-Ym-Zm坐标系；
由Xm-Ym-Zm坐标系获得方式，可知Xm-Ym平面中的点(xm,ym)与Xf-Yf平面中的点(xf,yf)的关系为：
xm＝xf，ym＝yf              (5)
由公式(3)，(4)和(5)可得X-Y平面上点与Xm-Ym平面上点的关系：
xm＝x，ym＝ycosβ              (6)
由式(6)可知，Xm-Ym平面上点与X-Y平面上点的关系，因此，a，b，c，d四条边界线在Xm-Ym坐标系下的投影方程的求解问题转化为这四条边界线在X-Y坐标系下的投影方程的求解问题；
3.5)确定a，b，c，d四条线在X-Y坐标系下的投影方程
3.5.1)a号线在X-Y坐标系下投影方程求解过程为：
a号线为铣刀球头与工件加工表面的交线；工件加工表面在X-Y-Z坐标系下的方程需要通过X0-Y0-Z0坐标系与X-Y-Z坐标系的关系获得；
在X0-Z0坐标系下，在工件加工表面任取A(xA0,zA0)、B(xB0,zB0)两点，由于工件加工面为平面，故zA0＝zB0＝L_jg；通过X–Z坐标系与X0-Z0坐标系的关系，可获得A、B两点在X–Z坐标系的值(xA,zA)，(xB,zB)，具体获取方式如下：
获取A、B在X–Z坐标系下坐标之后，计算出经过A、B两点在X–Z坐标系下的直线方程z＝kx+b，其中 此方程在X–Z坐标系下为直线方程，在X-Y-Z下可
表示为工件加工表面方程；
在X-Y-Z下，将工件加工表面方程z＝kx+b与铣刀球头方程x2+y2+(z-R)2＝R2方程联立，即可得到a号线在X-Y二维坐标系下投影方程x2+y2+(kx+b-R)2＝R2；
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在Xm-Ym坐标系下的投影方程
利用x＝xm，y＝ym/cosβ分别替换掉a、b、c、d在X-Y坐标系下投影方程中的x、y，获得的新的方程为a、b、c、d在Xm-Ym坐标系下的投影方程，它们共同围成的区域即为前倾角为α、侧倾角为β时球头铣刀与工件接触区域在Xm-Ym坐标系下的投影；
步骤四、精细积分法对刀具-工件动力学方程时域数值求解
4.1)精细积分法对刀具-工件动力学方程进行时域数值求解
由公式(2)将公式(1)表示为如下形式：
将式(7)表示为如下所示的哈密顿系统：
其中，
将式(8)中A(t)v(t)-A(t)v(t-T)用f(t)来表示，则对于非齐次方程(8)，由常微分理论可知，一般解为：
将时滞周期T均分为m份，即T＝m·τ，在[tp,tp+1]中，将f(t)表示为如下形式：
f(t)＝r0+r1(t-tp)                      (10)
其中，r0＝f(tp)＝A(tp)v(tp)-A(tp)v(tp-m·τ)；
由式(9)和式(10)可将v(tp+1)表示为：
v(tp+1)＝T·[v(tp)+A0-1(r0+A0-1r1)]-A0-1(r0+A0-1r1+r1·τ)           (11)其中， 将τ均分为Λ＝220份，则
由于 足够小，故用Taylor展开式去近似表示为
由式(12)和式(13)可知，
Ta＝A0Δt+(A0Δt)2/2！+(A0Δt)3/3！+(A0Δt)4/4！             (15)
将r0和r1进一步分别表示为
r0＝Apvp-Apvp-m                            (16)
将(16)和(17)带入到式(11)中，可得：
其中，MM＝TA0-1-A0-1，NN＝TA0-2-A0-2-A0-1τ；
若(I-NN/τ·Ap+1)可逆，则式(18)可表示为：
其中，
4.2)确定axx,p、axy,p、ayx,p、ayy,p和axx,p+1、axy,p+1、ayx,p+1、ayy,p+1以Nf＝2阐述axx,p、axy,p、ayx,p、ayy,p的确定方法；在tp时刻，将球头铣刀全部切削刃投影到Xm-Ym坐标系下，其中， n为刀具转速，j从1到Nf，k从0到π/2；由铣刀刀齿数和构建出的刀具-工件的接触区域可知在时滞周期内的任意时刻，只有一条切削刃参与切削，设此切削刃为第一条切削刃，该切削刃在单齿切削周期内到达与b号线相切位置、b投影线与a投影线交点位置、与d号线相切位置、a投影线与d投影线交点位置、d投影线与c投影线交点位置、c投影线与b投影线交点位置、时滞周期切削结束的时刻分别记为t1、t2、t3、t4、t5、t6、t7；
当tp＝0时，第一条切削刃没有参与切削；当tp在0-t1时间内，第一条切削刃没有参与切削，此时kmax,1和kmin,1值都为零；当tp在t1-t2时间内，第一条切削刃与b号线投影线相交，两个交点所对应的轴向角中的最大值即为tp时刻该切削刃在公式(20)中所对应的kmax，1，最小值为tp时刻该切削刃在公式(20)中所对应的kmin,1；同理，能够确定tp分别在t2-t3，t3-t4，t4-t5，t5-t6时间段内，第一条切削刃所对应kmax,1和kmin,1值；当tp在t6-t7时间内，第一条切削刃没有参与切削，此时kmax,1和kmin,1值都为零；
确定出各个切削刃在tp参与切削的最大轴向角和最小轴向角之后，通过式(20)得到axx,p、axy,p、ayx,p、ayy,p值；
通过以上步骤得到tp+1时刻时axx,p+1、axy,p+1、ayx,p+1、ayy,p+1值；
步骤五、叶瓣图构建
5.1)切削稳定性判定方法
建立系数矩阵Cp，该矩阵满足离散映射：
vp+1＝Cpvp                         (21)
其中，vp是个(2m+4)维的向量，表示为：
矩阵Cp为(2m+4)维矩阵，表示为：
矩阵PK为4×4矩阵等于式(19)中的(I-NN/τ·Ap+1)-1(T+MM·Ap-NN/τ·Ap)，矩阵RK1为4×2矩阵，等于式(19)中的-(I-NN/τ·Ap+1)-1NN/τ·Ap+1的前两列，RK2为4×2矩阵，等于式-1
(19)中的(I-NN/τ·Ap+1) (NN/τ·Ap-MM·Ap)的前两列；
通过使用一系列离散Cp，p＝0,1,2…,m-1，构建时滞周期T内的过渡矩阵Φ，亦即：
vp＝Φv0                (24)
式中，Φ定义为：Φ＝Cm-1Cm-2…C1C0；
由Floquet理论可知，传递函数Φ所有特征值模的最大值小于1，等于1和大于1，分别表示在该刀具转速n和轴向切深L_jg下，切削处于稳定状态，临界状态和不稳定状态；
5.2)叶瓣图构建
主轴转速不变，不断改变刀具轴向切深，判断其所对应的切削状态，获得在[0,L_jg]范围内的临界轴向切深；改变轴向切深后的接触区域边界投影中只有a号线投影方程变化，b号线和c号线投影方程仍与L_jg时相同；改变主轴转速，获得不同主轴转速在[0,L_jg]范围内所对应的临界轴向切削深度；最终，构建出临界轴向切削深度随主轴转速变化的函数关系，即为颤振稳定域叶瓣图。
2.根据权利要求1所述的五轴数控机床加工中基于精细积分的球头铣刀颤振稳定域叶瓣图建模方法，其特征在于，所述的步骤二求解球头铣刀刀齿上的动态切削力Ftx(t)和Fty(t)的具体步骤如下：
2.1)建立球头铣刀第j刀刃上第i个切削微元坐标表达式如下：
其中，R为球头铣刀半径；μ切削刃螺旋角；t为切削过程中的时间，单位为s；k为第j刀刃上第i个切削微元的轴向接触角，在一个切削刃上所能取值范围为[0，π/2]；ψji(k)为第j刀刃上第i个切削微元径向滞后角；φ10(t)为第一个切削刃端点处转动的角度，n为刀具转速，单位为r/min；φji(t)为第j刃上的第i个微元处瞬时径向接触角；Nf为刀齿数目；xji(t)，yji(t)，zji(t)表示球头铣刀第j刀刃上第i个切削微元在所建立的坐标系下的坐标值；
2.2)球头铣刀瞬时动态切削力计算
球头铣刀第j刀刃上第i个切削微元轴向角为k所受的切向力dFt,ji(φji(t),k)、径向力dFr,ji(φji(t),k)、轴向力dFa,ji(φji(t),k)表示为：
其中，h(φji(t),kji)为第j刀刃上第i个切削微元瞬时切削厚度；Ktc为切向力系数，Krc为径向力系数，Kac为轴向力系数；db＝R·dk，R为球头铣刀半径；
2.3)瞬时动态切削厚度计算
刀具瞬时铣削厚度包括瞬时静态切削厚度和瞬时动态切削厚度，瞬时静态切削厚度与颤振无关，故将其忽略，仅考虑刀具x和y方向的振动，球头铣刀第j刀刃上第i个切削微元，轴向角为k的瞬时动态铣削厚度为：
其中，x(t)-x(t-T)，y(t)-y(t-T)表示当前时刻t和前一个刀齿切削(t-T)时刻在x和y方向的动态振动矢量；T为时滞周期，在高速切削条件下，时滞周期等于单齿切削周期，即时滞周期 Nf为铣刀刀齿数；n为刀具转速，单位为r/min；
2.4)刀具瞬时动态切削力计算
通过坐标变换，获得第j刀刃上第i个切削微元x，y方向动态切削力：
由式(26)，式(27)和式(28)得到如下结果：
其中，axx,ji(t)、axy,ji(t)、ayx,ji(t)、ayy,ji(t)如下所示；
确定每个切削刃在t时刻的参与切削片段的数目和每个参与切削的片段所对应的最大轴向角和最小轴向角，通过公式(2)就得到球头铣刀上的动态切削力；
3.根据权利要求1或2所述的五轴数控机床加工中基于精细积分的球头铣刀颤振稳定域叶瓣图建模方法，其特征在于，所述的步骤三的3.5)确定a，b，c，d，四条线在X-Y坐标系下的投影方程的具体步骤为：
3.5.1)b号线投影方程求解过程
b号线为球头铣刀垂直于进给方向的轮廓圆截面与工件过渡表面的交线，此圆投影到X-Y二维坐标系下为一椭圆，椭圆方程为
3.5.2)c号线投影方程求解过程
c号线为本次走刀与前一走刀在工件已加工表面留下的加工残余而形成的线，在Y1-Z1坐标系下，获得铣刀球头投影方程y12+(z1-R)2＝R2，将其沿着Y1负方向平移L_xl，得到方程(y1+L_xl)2+(z1-R)2＝R2，即为与本次刀触点相对应的上一次加工时的铣刀球头轮廓方程；
两圆交点M所对的Y1值为-L_xl/2，即为c号线在X1-Y1二维坐标系下的投影，由于X-Y-Z坐标系与X1-Y1-Z1坐标系的X-Z平面与X1-Z1平面相重合，Y轴与Y1轴相互投影方程平行，故c号线在X-Y二维坐标系下的投影方程为y＝-L_xl/2；
3.5.3)d号线求解过程
d号线投影方程求解过程与a号和e号线投影方程有关；
1)e号线投影方程求解过程
e号线为前一次走刀在工件加工表面留下的加工痕迹；在Y1-Z1坐标系下，将z1＝L_jg带入到与本次刀位点相对应的上一次加工时的铣刀球头轮廓方程(y1+L_xl)2+(z1-R)2＝R2中，得到 的最 大y 1即 为 e号 线在 X1 -Y 1 坐标 系下的 投 影 ，经计 算该 值 为由于X-Y-Z坐标系与X1-Y1-Z1坐标系的X-Z平面与X1-Z1平面
相重合，Y轴与Y1轴相互平行，因此，e号线在X-Y二维坐标系下的投影方程为
2)d号线投影方程求解过程
d号线上的N点在X-Y坐标系下的投影可由a号线和e号线在X-Y坐标系下投影方程联立得到，经计算N点x值满足方程 y值为
利用分层思想可知，d号线上的其它点可以认为是不同轴向切削深度所对应的a号线与e 号 线的 交 点 ，因 此 将 N 点 坐 标 中的 L _ jg 用变 量 k _ v 替 换 ，可 得 ，单纯从数
学角度观察N点坐标，可知在其它值一定的情况下，N点的坐标值只与变量k_v值有关，改变k_v值，就能得到d号线上其它点的投影坐标；
k_v能够取得的最大取值为L_jg，能取的最小取值为e号线在Y1-Z1坐标系投影所对应的Z1坐标，即M点所对用的Z1值 为了精准表示d号线投影方程，决定在X-Y
平面上沿Y轴等间距在d号线投影上取7个点的坐标值，具体如下所示：
d号线在X1-Y1平面投影所对应的最大Y1值为 最小Y1值为
从Y1的最大值和最小值中等间距取7值，通过方程
取x大于零的解
和 获取7组所对应的x，y值，即可得到d号线在X-Y平面上沿Y轴等
间距分布的7个点，对这7个坐标值进行牛顿插值即获取在X-Y坐标系下的六次多项式方程。