本發(fā)明屬于機械制造技術(shù)領(lǐng)域，涉及到五軸數(shù)控機床銑削切削過程中顫振穩(wěn)定域建模方法，尤其是球頭銑刀在五軸數(shù)控機床銑削切削過程中顫振穩(wěn)定域葉瓣圖建模方法。

背景技術(shù)：

機械制造業(yè)作為國家經(jīng)濟增長支柱，已經(jīng)發(fā)展了上百年，建立了比較系統(tǒng)的理論體系，積累了豐富的實踐經(jīng)驗，但隨著科學(xué)技術(shù)水平的提高，迫切的需要其向高速、高效和高精度方向發(fā)展。

高速銑削加工是先進制造技術(shù)中最重要的基礎(chǔ)技術(shù)之一，是目前最重要、應(yīng)用最普遍的加工方式，然而在高速銑削過程中不合適的切削參數(shù)導(dǎo)致的切削顫振嚴重地影響到了加工效率、精度、質(zhì)量以及穩(wěn)定性，是制約高速銑削技術(shù)快速發(fā)展的關(guān)鍵因素。

國內(nèi)外學(xué)者對切削顫振進行了大量探索，研究了多種顫振形成的機理，其中再生型顫振是人們認為產(chǎn)生切削過程中產(chǎn)生顫振最直接、最根本的原因。如圖1所示，再生型顫振理論指出，由于機床結(jié)構(gòu)振動，當(dāng)?shù)毒哌M行切削時，工件的已加工表面會留下表面振紋，當(dāng)?shù)毒咴僖淮吻邢鞯竭@些遺留有振紋的工件表面時，瞬時切削厚度由名義切削厚度和動態(tài)切削厚度疊加組成，這種切削厚度的變化引起切削力的波動，反過來又引起切削刀具和工件之間的相對振動，使刀具和工件在切削過程中產(chǎn)生振動位移，從而再次在工件已加工表面留下振紋，根據(jù)前后兩次振紋之間的相位差，在靠近但不等于加工系統(tǒng)主結(jié)構(gòu)模態(tài)的顫振頻率處，隨著加工系統(tǒng)的切削厚度不斷增長，造成切削力和振動位移的不斷上升現(xiàn)象，形成強烈的自激振動，這種自激振動就是再生顫振，圖2為再生顫振機理模型。

避免切削顫振最有效辦法是在加工前構(gòu)建銑削顫振穩(wěn)定域葉瓣圖，即在給定切削條件下，繪制出軸向臨界切削深度隨主軸轉(zhuǎn)速變化的函數(shù)關(guān)系。葉瓣圖的構(gòu)建能夠為加工前切削參數(shù)的選擇提供指導(dǎo)，可以有效防止加工中顫振的發(fā)生。而球頭銑刀是典型的點接觸加工，具有很好的法矢自適應(yīng)性，是高速銑削加工中應(yīng)用最廣的刀具之一，因此構(gòu)造出針對球頭銑刀加工過程中的穩(wěn)定性葉瓣圖意義重大。

目前，國內(nèi)外還沒有針對球頭銑刀在五軸數(shù)控機床加工過程中的顫振穩(wěn)定域葉瓣圖建模方法。有部分學(xué)者對三軸數(shù)控機床上球頭銑刀顫振穩(wěn)定域葉瓣圖建模方法進行了研究，其中具有代表性的有：李忠群等人提出的基于龍格庫塔的時域求解法；dingh等人提出的全離散時域求解法；limingzhen等人提出的完全離散時域數(shù)值求解法。但是不論是龍格庫塔法，全離散法，還是改進的完全離散法本質(zhì)上都屬于差分類方法，它們的優(yōu)點是易于實現(xiàn)，但差分引起的誤差很難從根本上消除。同時，五軸數(shù)控機床和三軸數(shù)控機床的加工方式不同，刀具姿態(tài)不同，故球頭銑刀在三軸數(shù)控機床的建模方法并不能適用于其在五軸數(shù)控機床上；因此如何在更先進的五軸數(shù)控機床上對球頭銑刀進行高精度的顫振穩(wěn)定域葉瓣圖構(gòu)建將是一個亟待解決的問題。

技術(shù)實現(xiàn)要素：

為了解決現(xiàn)有技術(shù)存在的問題，本發(fā)明提供一種五軸數(shù)控機床加工中基于精細積分的球頭銑刀顫振穩(wěn)定域葉瓣圖建模方法。

本發(fā)明的技術(shù)方案為：

一種五軸數(shù)控機床加工中基于精細積分的球頭銑刀顫振穩(wěn)定域葉瓣圖建模方法，包括以下步驟：

步驟一，建立球頭銑刀刀具-工件動力學(xué)方程；

步驟二，求解球頭銑刀刀齒上的動態(tài)切削力ftx(t)和fty(t)；

步驟三，五軸數(shù)控機床平面加工中球頭銑刀與工件的接觸區(qū)域半解析建模；

步驟四，精細積分法對刀具-工件動力學(xué)方程時域數(shù)值求解；

步驟五，葉瓣圖構(gòu)建。

本發(fā)明的有益效果為：運用精細積分法對銑削系統(tǒng)二階動力學(xué)方程進行高精度時域數(shù)值求解，克服了傳統(tǒng)數(shù)值求解法不能同時兼顧計算精度、計算效率和穩(wěn)定性的弊端；精準識別出球頭銑刀存在前傾角和側(cè)傾角時與工件接觸區(qū)域邊界的投影方程，在單齒切削周期內(nèi)球頭銑刀視為圓弧切削的基礎(chǔ)上，通過球頭銑刀與工件的接觸區(qū)域投影邊界方程和切削刃不同時刻的投影方程的關(guān)系，確定出時域數(shù)值方程中所需要的瞬時參與切削的刀刃數(shù)目和參與切削刀刃的實際切削部位，利用floquet定理獲得了不同轉(zhuǎn)速下臨界切削深度，構(gòu)建出了球頭銑刀存在前傾角和側(cè)傾角時的顫振穩(wěn)定域葉瓣圖。

附圖說明

圖1是球頭銑刀銑削過程中再生顫振發(fā)生機理示意圖。

圖2是球頭銑刀銑削過程中再生顫振機理模型。

圖3是球頭銑刀銑削動力學(xué)模型。

圖4是球頭銑刀螺旋切削刃幾何模型。

圖5是球頭銑刀-工件接觸區(qū)域示意圖。

圖6是球頭銑刀-工件接觸區(qū)域邊界組成示意圖。

圖7是球頭銑刀三維坐標(biāo)x0-y0-z0建立示意圖。

圖8是推導(dǎo)出刀具前傾角為α，側(cè)傾角為零時的刀位點和刀觸點示意圖。

圖9是球頭銑刀三維坐標(biāo)x1-y1-z1和x-y-z建立過程示意圖。

圖10是球頭銑刀三維坐標(biāo)xp-yp-zp建立過程示意圖。

圖11是球頭銑刀三維坐標(biāo)xf-yf-zf建立過程示意圖。

圖12是球頭銑刀三維坐標(biāo)xm-ym-zm建立過程示意圖。

圖13是三維坐標(biāo)xm-ym-zm示意圖。

圖14是b號投影線求解示意圖。

圖15是m點示意圖。

圖16是n點示意圖。

圖17是d號線投影求解中沿y1軸等間距取值示意圖。

圖18是d號線投影方程求解過程示意圖。

圖19是接觸區(qū)域在xm-ym坐標(biāo)系投影示意圖。

圖20是第一條切削刃投影0-t5時刻位置。

圖21是tp＝0時刻，切削刃參與切削情況。

圖22是tp＜t1時刻，切削刃參與切削情況。

圖23是t1＜tp＜t2時刻，切削刃參與切削情況。

圖24是t2＜tp＜t3時刻，切削刃參與切削情況。

圖25是t3＜tp＜t4時刻，切削刃參與切削情況。

圖26是t4＜tp＜t5時刻，切削刃參與切削情況。

圖27是t5＜tp＜t6時刻，切削刃參與切削情況。

圖28是t6＜tp＜t7時刻，切削刃參與切削情況。

圖29是軸向切削深度小于l_jg時的接觸區(qū)域在xm-ym坐標(biāo)系下投影的邊界組成。

圖30是葉瓣圖構(gòu)建流程圖。

具體實施方式

以下結(jié)合技術(shù)方案和附圖詳細敘述本發(fā)明的具體實施方式。

五軸數(shù)控機床上加工中基于精細積分的球頭銑刀顫振穩(wěn)定域葉瓣圖建模方法，包括以下步驟：

步驟一，建立球頭銑刀刀具-工件動力學(xué)方程

將球頭銑刀刀具-工件系統(tǒng)簡化為如圖3所示的二自由度系統(tǒng)，僅考慮進給方向x和法向y方向的刀具振動因素，建立如下所示的動力學(xué)方程：

其中，mtx，ξx，ωnx分別是為刀具系統(tǒng)x方向的模態(tài)質(zhì)量，阻尼系數(shù)，固有頻率；mty，ξy，ωny分別是為刀具系統(tǒng)y方向的模態(tài)質(zhì)量，阻尼系數(shù)，固有頻率；ftx(t)和fty(t)分別為在x，y方向上作用在銑刀上的動態(tài)切削力。

步驟二，求解球頭銑刀刀齒上的動態(tài)切削力ftx(t)和fty(t)

2.1建立球頭銑刀刀刃幾何模型

如圖4所示，建立球頭銑刀切削刃曲線幾何模型，由于每個切削微元與其所對應(yīng)的軸向角存在一一映射關(guān)系，故切削微元坐標(biāo)值可以表示為關(guān)于軸向角的函數(shù)，式(2)為第j刀刃上第i個切削微元坐標(biāo)表達式。

其中，r為球頭銑刀半徑；μ切削刃螺旋角；t為切削過程中的時間(s)；k為第j刀刃上第i個切削微元的軸向接觸角，在一個切削刃上所能取值范圍為[0，π/2]；ψji(k)為第j刀刃上第i個切削微元徑向滯后角；φ10(t)為第一個切削刃端點處轉(zhuǎn)動的角度，n為刀具轉(zhuǎn)速(r/min)；φji(t)為第j刃上的第i個微元處瞬時徑向接觸角；nf為切削刃數(shù)目；xji(t)，yji(t)，zji(t)表示球頭銑刀第j刀刃上第i個切削微元在所建立的坐標(biāo)系下的坐標(biāo)值。

2.2球頭銑刀瞬時動態(tài)切削力計算

球頭銑刀第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為切向力系數(shù)，krc為徑向力系數(shù)，kac為軸向力系數(shù)；db＝r·dk，r為球頭銑刀半徑。

2.3瞬時動態(tài)切削厚度計算

刀具瞬時銑削厚度由兩部分組成，一部分是瞬時靜態(tài)切削厚度，另一部分是瞬時動態(tài)切削厚度，瞬時靜態(tài)切削厚度與顫振無關(guān)，故將其忽略。

僅考慮刀具x和y方向的振動，球頭銑刀第j刀刃上第i個切削微元(軸向角為k)的瞬時動態(tài)銑削厚度為：

其中，x(t)-x(t-t)，y(t)-y(t-t)表示當(dāng)前時刻t和前一個刀齒切削(t-t)時刻在x和y方向的動態(tài)振動矢量；t為時滯周期，在高速切削條件下，時滯周期等于單齒切削周期，即時滯周期nf為銑刀刀齒數(shù)；n為刀具轉(zhuǎn)速(r/min)。

2.4刀具瞬時動態(tài)切削力計算

通過坐標(biāo)變換，獲得第j刀刃上第i個切削微元x，y方向動態(tài)切削力：

由式(3)，式(4)和式(5)得到如下結(jié)果：

其中，axx,ji(t)、axy,ji(t)、ayx,ji(t)、ayy,ji(t)如下所示。

確定每個切削刃在t時刻的參與切削片段的數(shù)目和每個參與切削的片段所對應(yīng)的最大軸向角和最小軸向角，通過公式(6)就得到球頭銑刀上的動態(tài)切削力。

其中，

步驟三、五軸數(shù)控機床平面加工過程中球頭銑刀與工件接觸區(qū)域的半解析建模

球頭銑刀和工件的接觸區(qū)域是指在加工過程中刀具切入工件的區(qū)域。高速銑削條件下，在單齒切削周期內(nèi)，球頭銑刀可視為圓弧切削，因此銑刀-工件接觸區(qū)域蘊含著式(6)中所需要的單齒切削周期內(nèi)銑刀瞬時參與切削的刀刃數(shù)目和參與切削刀刃的實際切削部位等信息。

3.1進行刀具路徑規(guī)劃，設(shè)置加工參數(shù)

球頭銑刀在五軸數(shù)控機床上進行平面銑削，球頭銑刀球頭半徑為r，前傾角為α度，側(cè)傾角為β度，相鄰切削刀軌間距為l_xl，軸向切削深度為l_jg。

3.2確定前傾角為α度，側(cè)傾角為β度時球頭銑刀與工件接觸區(qū)域邊界組成，將接觸區(qū)域邊界方程的求解問題轉(zhuǎn)化為這些邊界在垂直于刀具軸線平面的投影方程的求解問題；

如圖5所示，前傾角為α，側(cè)傾角為β度時球頭銑刀與工件的接觸區(qū)域主要由圖6中所示的a，b，c，d四條線組成，其中，a號線為銑刀球頭與工件加工表面的交線，b號線為銑刀球頭與工件過渡表面的交線，c號線為本次走刀與上一次走刀在已加工表面留下的加工殘余而形成的線，d號線為銑刀球頭與上一次走刀留下的加工痕跡的交線。另外e號線為上一次走刀在工件加工表面留下的加工痕跡。

在切削靜力學(xué)和切削動力學(xué)模型中，確定刀刃某一時刻參加切削的上下邊界時都是在二維平面上進行的，因此根據(jù)球頭銑刀球頭部分沿垂直刀軸的方向具有單調(diào)分布的特點，將a，b，c，d四條線求解問題轉(zhuǎn)化為它們在垂直于刀具軸線平面的投影方程的求解問題。

3.3由刀具前傾角和側(cè)傾角為零度時的三維直角坐標(biāo)系x0-y0-z0推導(dǎo)出刀具前傾角為α，側(cè)傾角為零度時的三維直角坐標(biāo)系x-y-z

如圖7所示，建立刀具前傾角和側(cè)傾角都為零度時的三維直角坐標(biāo)系x0-y0-z0，其中刀具頂點為原點o0，刀軸線為z0軸，x0軸與刀具進給方向相同。下面將由x0-y0-z0推導(dǎo)出刀具前傾角為α，側(cè)傾角為零度時的三維直角坐標(biāo)系x-y-z，由于側(cè)傾角都為零度，故只需要在x0-z0平面討論即可。

如圖8所示，在x0-z0平面內(nèi)，將z0軸以o0為原點，順時針傾斜α度(從y0負方向向原點看，α值為正數(shù)，順時針；α值為負數(shù)，逆時針，此處以α為正進行闡述)，傾斜之后的軸線就是刀具前傾角為α度，側(cè)傾角為零度時的刀具軸線，該直線方程z0＝tan(90-α)×x0與直線z0＝r的交點cir_0，即為刀具前傾角為α度，側(cè)傾角為零度時的球心，經(jīng)計算cir_0的坐標(biāo)為xcir_o＝r/tan(90-α)，zcir_0＝r。

以cir_0為原點，r為半徑建立圓的方程(x0-xcir_0)²+(z0-ycir_0)²＝r²，此圓方程與x0軸的切點即為前傾角為α度，側(cè)傾角為零度時的球頭銑刀的刀觸點dc_0，經(jīng)計算dc_0的坐標(biāo)xdc_0＝r/tan(90-α)，zdc_0＝0。

(x0-xcir_0)²+(z0-ycir_0)²＝r²方程與z0＝tan(90-α)×x0方程兩個交點中y0值較小的點即為前傾角為α度，側(cè)傾角為0度時的球頭銑刀的刀位點dw_0。經(jīng)計算xdw_0＝r/tan(90-α)-r×sinα，zdw_0＝r-r×cosα。

如圖9所示，以dc_0為原點，建立三維直角坐標(biāo)系x1-y1-z1，其中，x1-z1平面與x0-z0平面相重合，x1與x0相互重合且方向相同，y1與y0平行且方向相同。

如圖9所示，以dw_0為原點，方程z0＝tan(90-α)×x0為z軸建立三維直角坐標(biāo)系x-y-z，其中x-z平面與x0-z0平面相重合，y與y0平行且方向相同，z軸正方向為遠離工件。x-y-z坐標(biāo)系下，銑刀球頭輪廓方程x²+y²+(z-r)²＝r²。

3.4通過x-y-z三維坐標(biāo)系，建立刀具前傾角為α，側(cè)傾角為β度時的xm-ym-zm坐標(biāo)系，獲取x-y平面上點與xm-ym平面上點的關(guān)系

如圖10所示，將x-y-z坐標(biāo)系中的x-y平面沿著z軸正方向平移r距離，使坐標(biāo)系的原點o與銑刀球心重合，獲得新的坐標(biāo)系xp-yp-zp。其中xp-yp-zp坐標(biāo)系中的xp軸與x軸相互平行，且方向相同；yp軸與y軸相互平行，且方向相同；zp軸與z軸相互重合，且方向相同。

由xp-yp-zp坐標(biāo)系獲得方式可知，xp-yp平面中的點(xp,yp)與x-y平面中的點(x,y)的關(guān)系為:

xp＝x，yp＝y(tǒng)(7)

如圖11所示，以xp軸為軸線，將yp-zp平面逆時針旋轉(zhuǎn)β度(xp軸正半軸向原點看，β值為正數(shù)，順時針；β值為負數(shù)，逆時針，以下以β為負值進行闡述)，旋轉(zhuǎn)之后，yp軸變?yōu)閥f軸，zp軸變?yōu)閦f軸，形成一個新的坐標(biāo)系xf-yf-zf，其中xf與xp軸重合且方向相同。

由xf-yf-zf坐標(biāo)系獲得方式可知，xp-yp-zp坐標(biāo)系中點(xp,yp,zp)與xf-yf-zf坐標(biāo)系中的點(xf,yf,zf)關(guān)系為：

令zf＝zp＝0，則可得：

xf＝xp，yf＝y(tǒng)pcosβ(8)

如圖12所示，將坐標(biāo)系xf-yf-zf中的xf-yf平面沿著zf負方向平移r距離，得到xm-ym-zm坐標(biāo)系，其中，xm軸與xf軸相互平行，且方向相同；ym軸與yf軸相互平行，且方向相同；zm軸與zf軸相互重合，且方向相同。如圖13所示，xm-ym-zm坐標(biāo)系的原點就是前傾角為α度，側(cè)傾角為β度時球頭銑刀的頂點，zm軸就是前傾角為α度，側(cè)傾角為β度時球頭銑刀的軸線，a，b，c，d四條邊界的求解問題轉(zhuǎn)化為它們在xm-ym平面的投影方程的求解問題。

由xm-ym-zm坐標(biāo)系獲得方式，可知xm-ym平面中的點(xm,ym)與xf-yf平面中的點(xf,yf)的關(guān)系為：

xm＝xf，ym＝y(tǒng)f(9)

由公式(7)，(8)和(9)可得x-y平面上點與xm-ym平面上點的關(guān)系：

xm＝x，ym＝y(tǒng)cosβ(10)

經(jīng)過以上推導(dǎo)可知，前傾角為α度，側(cè)傾角為零度時和前傾角為α度，側(cè)傾角為β度時，銑刀球頭的位置并沒有變化，刀觸點沒有變化，只是刀軸線變化了，刀位點發(fā)生了變化，也即目標(biāo)坐標(biāo)系由x-y-z變成了xm-ym-zm，這主要是由銑刀球頭的幾何特性所形成的。

由式(10)可知xm-ym平面上點與x-y平面上點的關(guān)系，因此，a，b，c，d四條線邊界曲線在xm-ym坐標(biāo)系下的投影方程的求解問題可以轉(zhuǎn)化為這四條線邊界線在x-y坐標(biāo)系下的投影方程的求解問題。

3.5確定a，b，c，d，e四條線在x-y坐標(biāo)系下的投影方程

3.5.1a號線在x-y坐標(biāo)系下投影方程求解過程

a號線為銑刀球頭與工件加工表面的交線。如圖8所示，工件加工表面在x-y-z坐標(biāo)系下的方程需要通過x0-y0-z0坐標(biāo)系與x-y-z坐標(biāo)系關(guān)系的獲得。

如圖9所示，在x0-z0坐標(biāo)系下，在工件加工表面任取a(xa0,za0)、b(xb0,zb0)兩點，由于工件加工面為平面，故za0＝zb0＝l_jg。通過x–z坐標(biāo)系與x0-z0坐標(biāo)系的關(guān)系，可獲得a、b兩點在x–z坐標(biāo)系的值(xa,za)，(xb,zb)，具體獲取方式如下：

獲取a、b在x–z坐標(biāo)系下坐標(biāo)之后，計算出經(jīng)過a、b兩點在x–z坐標(biāo)系下的直線方程z＝kx+b，其中此方程在x–z坐標(biāo)系下為直線方程，在x-y-z下可表示為工件加工表面方程。

在x-y-z下，將工件加工表面方程z＝kx+b與銑刀球頭方程x²+y²+(z-r)²＝r²方程聯(lián)立，即可得到a號線在x-y二維坐標(biāo)系下投影方程x²+y²+(kx+b-r)²＝r²。

3.5.2b號線投影方程求解過程

b號線為球頭銑刀垂直于進給方向的輪廓圓截面與工件過渡表面的交線，如圖14所示，此圓投影到x-y二維坐標(biāo)系下為一橢圓，橢圓方程為

3.5.3c號線投影方程求解過程

c號線為本次走刀與上一次走刀在工件已加工表面留下的加工殘余而形成的線，如圖15所示，在y1-z1坐標(biāo)系下，獲得銑刀球頭投影方程y1²+(z1-r)²＝r²，將其沿著y1負方向平移l_xl，得到方程(y1+l_xl)²+(z1-r)²＝r²，即為與本次刀觸點相對應(yīng)的上一次加工時的銑刀球頭輪廓方程。兩圓交點m所對的y1值為-l_xl/2，即為c號線在x1-y1二維坐標(biāo)系下的投影，由于x-y-z坐標(biāo)系與x1-y1-z1坐標(biāo)系的x-z平面與x1-z1平面相重合，y軸與y1軸相互投影方程平行，故c號線在x-y二維坐標(biāo)系下的投影方程為y＝-l_xl/2。

3.5.4d號線求解過程

d號線投影方程求解過程同a號和e號線投影方程有關(guān)，在求解d號線投影方程之前，先確定e號線投影方程。

1)e號線投影方程求解過程

e號線為上一次走刀在工件加工表面留下的加工痕跡。如圖15所示，在y1-z1坐標(biāo)系下，將z1＝l_jg帶入到與本次刀位點相對應(yīng)的上一次加工時的銑刀球頭輪廓方程(y1+l_xl)²+(z1-r)²＝r²中，得到的最大的y1即為e號線在x1-y1坐標(biāo)系下的投影，經(jīng)計算該值為由于x-y-z坐標(biāo)系與x1-y1-z1坐標(biāo)系中的x-z平面與x1-z1平面相重合，y軸與y1軸相互平行，因此，e號線在x-y二維坐標(biāo)系下的投影方程為

2)d號線投影方程求解過程

圖16中的d號線上的n點在x-y坐標(biāo)系下的投影可由a號線和e號線在x-y坐標(biāo)系下投影方程聯(lián)立得到，經(jīng)計算n點x值滿足方程y值為

利用分層思想可知，d號線上的其它點可以認為是不同軸向切削深度所對應(yīng)的a號線與e號的交點，因此將n點坐標(biāo)中的l_jg用變量k_v替換，可得單純從數(shù)學(xué)角度觀察n點坐標(biāo)，可知：在其它值一定的情況下，n點的坐標(biāo)值只與變量k_v值有關(guān)，改變k_v值，就能得到d號線上其它點的投影坐標(biāo)。

如圖17示，k_v最大取值為l_jg，最小取值為e號線在y1-z1坐標(biāo)系投影所對應(yīng)的z1坐標(biāo)，即m點所對用的z1值為了精準表示d號線投影方程，決定在x-y平面上沿y軸等間距在d號線投影上取7個點的坐標(biāo)值，具體如下所示：

如圖17所示，d號線在x1-y1平面投影所對應(yīng)的最大y1值為最小y1值為從y1的最大值和最小值中等間距取7值，通過方程(x大于零的解)和獲取7組所對應(yīng)的x，y值，即可得到d號線在x-y平面上沿y軸等間距分布的7個點，對這7個坐標(biāo)值進行牛頓插值即獲取如圖18所示的在x-y坐標(biāo)系下的六次多項式方程。

3.6確定a，b，c，d在xm-ym坐標(biāo)系下的投影方程

利用式(10)x＝xm，y＝y(tǒng)m/cosβ分別替換掉a，b，c，d在x-y坐標(biāo)系下投影方程中的x，y，獲得的新的方程就是a，b，c，d在xm-ym坐標(biāo)系下的投影方程，如圖19示，它們共同圍成的區(qū)域即為前傾角為α，側(cè)傾角為β時球頭銑刀與工件接觸區(qū)域在xm-ym坐標(biāo)系下的投影。

步驟四、精細積分法對刀具-工件動力學(xué)方程時域數(shù)值求解

4.1精細積分法對刀具-工件動力學(xué)方程時域數(shù)值分解

由式(6)可將式(1)表示為如下形式：

將式(11)表示為如下所示的哈密頓系統(tǒng)：

其中，

將式(12)中a(t)v(t)-a(t)v(t-t)用f(t)來表示，則對于非齊次方程(12)，由常微分理論可知，一般解為：

將時滯周期t均分為m份，即t＝m·τ，在[tp,tp+1]中，將f(t)表示為如下形式：

f(t)＝r0+r1(t-tp)(14)

其中，r0＝f(tp)＝a(tp)v(tp)-a(tp)v(tp-m·τ)

由式(13)和式(14)可將v(tp+1)表示

v(tp+1)＝t·[v(tp)+a0^-1(r0+a0^-1r1)]-a0^-1(r0+a0^-1r1+r1·τ)(15)

其中，將ε均分為λ＝2²⁰份，則

由于足夠小，故用taylor展開式去近似表示為

由式(16)和式(17)可知，

因此，用n次矩陣相稱就能得到t中的非單位矩陣部分，即

將r0和r1進一步分別表示為

r0＝apvp-apvp-m(20)

將(20)和(21)帶入到式(15)中，可得：

其中，mm＝ta0^-1-a0^-1，nn＝ta0^-2-a0^-2-a0^-1τ

若(i-nn/τ·ap+1)可逆，則式(22)可表示為

其中，

4.2axx,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時刻，按照公式(2)將球頭銑刀全部切削刃投影到第三部分所建立的xm-ym坐標(biāo)系下，其中n為刀具轉(zhuǎn)速，j從1到nf，k從0到π/2。

由銑刀刀齒數(shù)和構(gòu)建出的刀具-工件的接觸區(qū)域可知在時滯周期內(nèi)的任意時刻，只有一條切削刃參與切削，設(shè)此切削刃為第一條切削刃。如圖20所示，該切削刃在單齒切削周期內(nèi)到達與b號線相切位置；b投影線與a投影線交點位置；與d號線相切位置；a投影線與d投影線交點位置；d投影線與c投影線交點位置；c投影線與b投影線交點位置；時滯周期切削結(jié)束的時刻分別記為t1，t2，t3，t4，t5，t6，t7。

當(dāng)tp＝0時，如圖21所示，第一條切削刃沒有參與切削，當(dāng)tp在0-t1時間內(nèi)，如圖22所示，第一條切削刃沒有參與切削，此時kmax,1和kmin,1值都為零；當(dāng)tp在t1-t2時間內(nèi)，如圖23所示，第一條切削刃與b號線投影線相交，兩個交點所對應(yīng)的軸向角的最大值即為tp時刻該切削刃在式(24)中所對應(yīng)的kmax，1，最小值為tp時刻該切削刃在式(24)中所對應(yīng)的kmin,1。同理，可以確定tp分別在t2-t3，t3-t4，t4-t5，t5-t6時間段內(nèi)，第一條切削刃所對應(yīng)kmax,1和kmin,1值，具體如圖24-圖27所示；當(dāng)tp在t6-t7時間內(nèi)，如圖28所示，第一條切削刃沒有參與切削，此時kmax,1和kmin,1值都為零。

確定出各個切削刃在tp參與切削的最大軸向角和最小軸向角之后，通過式(24)可得到axx,p，axy,p，ayx,p，ayy,p值。

通過以上步驟可以得到tp+1時刻時axx,p+1，axy,p+1，ayx,p+1，ayy,p+1值。

步驟五、葉瓣圖構(gòu)建

5.1切削穩(wěn)定性判定方法

建立系數(shù)矩陣cp，該矩陣滿足離散映射：

vp+1＝cpvp(25)

vp是個(2m+4)維的向量：

矩陣cp為(2m+4)維矩陣：

pk為4×4矩陣等于公式(23)中的(i-nn/τ·ap+1)^-1(t+mm·ap-nn/τ·ap)，rk1為4×2矩陣，等于式公式(23)中的-(i-nn/τ·ap+1)^-1nn/τ·ap+1的前兩列，rk2為4×2矩陣，等于式(23)中的(i-nn/τ·ap+1)^-1(nn/τ·ap-mm·ap)的前兩列。

通過使用一系列離散cp(p＝0,1,2…,m-1)，構(gòu)建時滯周期t內(nèi)的過渡矩陣φ，亦即：

vp＝φv0(28)

式中，φ定義為：φ＝cm-1cm-2…c1c0。

由floquet理論可知，傳遞函數(shù)φ所有特征值模的最大值小于1，等于1和大于1，分別表示在該刀具轉(zhuǎn)速n和軸向切深l_jg下，切削處于穩(wěn)定狀態(tài)，臨界臨界狀態(tài)和不穩(wěn)定狀態(tài)。

5.2葉瓣圖構(gòu)建

主軸轉(zhuǎn)速不變，在[0,l_jg]范圍內(nèi)不斷改變軸向切深，通過5.1中的方法，獲得該轉(zhuǎn)速下所對應(yīng)的臨界軸向切深。改變軸向切深后的接觸區(qū)域投影邊界中只有a號線投影方程變化，b號線和c號線投影方程仍與l_jg時相同，具體如圖29所示。

改變主軸轉(zhuǎn)速，獲得不同主軸轉(zhuǎn)速在[0,l_jg]范圍內(nèi)所對應(yīng)的臨界軸向切深。最終，構(gòu)建出臨界軸向切深隨主軸轉(zhuǎn)速變化的函數(shù)關(guān)系，即為顫振穩(wěn)定域葉瓣圖，整個流程如圖30所示。