利用Kriging和Monte Carlo進行球頭銑刀銑削力概率分布預測

孫洪哲，孫志禮

（1.遼寧工業大學 機械工程與自動化學院，遼寧 錦州 121001；2.東北大學 機械工程與自動化學院，遼寧 沈陽 110819）

1 引言

球頭銑刀作為一種加工復雜曲面時的重要銑削刀具，為能夠在加工前確定合理的加工工藝，確定合理銑削參數，對球頭銑刀銑削力概率分布的研究越發重要。因此可以提高生產效率，降低加工成本等具有十分重要的現實意義。在加工過程中的銑削參數中，銑削力占據著重要的位置。銑削力的大小能夠影響銑削時的機床功率消耗并決定銑削時產生熱量的多少、刀具磨損的程度等。銑削力還能夠影響加工精度和加工質量，因此需要建立準確合理的銑削力模型。文獻[1]在斜角切削理論的基礎上考慮了刀刃上的耕犁力，進行銑削力的預測。文獻[2]利用球頭銑刀加工彎管內表面，并用預測結果與實際結果對比驗證預測的準確性。文獻[3]利用Z-Map法對求解了參與切削的切削刃片段，建模過程中考慮了刀具柔性以及刀具偏心的影響，可以很好地符合實際工況。文獻[4]基于BP神經網絡進行銑削力預測，并將預測結果與實驗結果進行比較，證明了BP神經網絡方法適用于對銑削力進行預測。

以上銑削力的預測方法中通常將銑削參數視為精確無誤差的確定性變量。然而，實際工程問題中由于機床制造過程中加工和裝配誤差、工作環境偶然因素等影響，可以精確描述的信息因素極少，大部分因素是非精確的、不完全的。確定性模型中不考慮不確定因素對銑削力的影響，基于此設計的銑削方案，在工程應用中可能會產生一定偏差，造成加工誤差過大、刀具磨損嚴重等不利結果。

以雙刃球頭銑刀等高加工剛性工件為主要研究對象，對銑削力進行建模，考慮隨機因素對銑削力的影響，提出一種基于文獻[5-6]的球頭銑刀銑削力概率分布計算方法，與解析方法進行比較驗證了提出方法的合理性，為后續分析球頭銑刀銑削精度打下基礎。

2 球頭銑刀銑削力建模

銑削力模型在刀具坐標系下建立的，如圖1所示。坐標原點為刀尖位置；Xc軸的正方向為銑刀進給方向在水平面內的投影方向；Zc軸正方向為沿銑刀軸線向上；Yc軸與Xc軸、Zc軸形成右手系。

根據幾何原理，螺旋導程不變的球頭銑刀切削刃的幾何形狀能夠簡單地看作是恒定螺旋角的立銑刀切削刃在球面上所產生的投影，如圖2所示。

圖1 刀具坐標系Fig.1 The Cutter Coordinate System

圖2 恒定導程球面螺旋線切削刃Fig.2 Constant Lead Spiral Cutting Edge for Ball End Mill

由于球頭銑刀銑削刃曲線十分復雜，刀刃上各點所處的銑削狀態不同。根據文獻[7]，將球頭銑刀的銑削刃分割成眾多微小的連續單元，并通過積分求得刀具整體銑削力。選用文獻[9]提出的微元銑削力表達式，如式（1）所示。

式中：dFt、dFr、dFa—切向、徑向和副法向的微元銑削力，N；Kte、Kre、Kae—微元切削刃切向、徑向和副法向的耕犁力系數，N/mm；Kts、Krs、Kas—微元切削刃切向、徑向和副法向的剪切力系數，N/mm2；dS—微元切削刃長度，mm；db—微元切屑寬度，mm；fn—微元切屑厚度，mm。

根據圖2的幾何關系，可以求得微元切削刃長度：

式中：R—刀具半徑；k—從Zc軸負方向向上度量的微元切削刃的軸向位置角；α—公稱螺旋角；微元切屑寬度可表示為db=Rdk。

微元切屑厚度模型選用Lee和Altinas提出的描述球頭銑刀水平進給的微元切屑厚度計算模型，即fn=fcsinksinθj

式中：θj＝ψj-φ（k），（klow＜k＜kup）。fn—微元切屑厚度，mm；fc—每齒進給量，mm/z；θj—第j號切削刃上微元點的水平位置角；ψj—第j號切削刃的位置角；kup、klow—切觸段上、下限。

將微元銑削力表達式投影到刀具坐標系下并對微元銑削力進行積分，再將各切削刃的銑削力求和，即可得到整體銑削力，如式（2）所示。

3 銑削力系數識別與刀刃切觸段確定

銑削力系數識別是銑削力建模的核心內容，是保證銑削力模型準確性的重要組成部分。根據文獻[8]，首先將銑削力模型整理成銑削力與銑削力系數之間的顯式線性函數，然后通過平面槽切實驗測取平均銑削力，再應用最小二乘法估算銑削力系數。

將耕犁力系數表示為常數，將剪切力系數表示為軸向位置角的3次多項式：

式中：Kis0、Kis1、Kis2、Kis3（i=t、r、a）—待定的常數系數。將耕犁力系數與剪切力系數代入式（1），并進行積分，得到單條切削刃銑削力的表達式：

根據文獻[8]球頭銑刀加工時，刀具和工件接觸區域的邊界形狀通常是不規則的空間曲線。利用解析法求出邊界曲線方程進而確定該接觸區域的方法十分復雜的。因此，采用基于離散仿真的Z-Map方法，具體方法可參考文獻[9]。銑削系數的識別方法由平面槽切實驗的方法來確定[9]。選取文獻[9]中HGT BLM0606硬質合金球頭銑刀對40Cr合金鋼工件進行平面槽切的實驗。測量平均銑削力并利用最小二乘法計算出銑削力系數得到銑削力的解析表達式。

4 球頭銑刀銑削力預測

Kriging模型包含確定性和隨機兩個部分，確定性部分一般采用最小二乘多項式擬合，隨機部分為高斯過程。與其它代理模型例如響應面法，神經網絡法相比，Kriging模型的優勢是在樣本點較少的情況下可以得到較精確的預測結果[5]。Monte Carlo模擬是一種基于概率統計的隨機模擬方法，能夠應用于多個機械產品領域，例如數控車床與數控銑床領域。Monte Carlo方法的主要理論基礎是概率論中的大數定理和中心極限定理[6]。應用Monte Carlo方法，對銑刀及加工參數等隨機參數進行大量的抽樣，并帶入到模型中進行計算得到樣本的響應值，即銑刀的切削力，并對響應值進行統計。在實際生產加工中，銑刀的形狀和加工參數具有隨機性，這些參數一般服從正態分布。由于球頭銑刀的銑削力是各參數的隱函數，球頭銑刀銑削力的解析形式也是一種對銑刀的銑削力的近似顯式表達式。為了驗證解析表達式的準確新，因此利用Kriging方法近似得出銑刀的銑削力其的另一個顯示表達式，并比較兩個表達式的結果。由于Kriging模型具有較快的收斂速率，因此采用銑刀的銑削力解析形式計算出平面槽切實驗測量50組實驗結果，得到50組X，Y，Z方向的平均銑削力，并用這些數據對Kriging模型進行訓練，選擇文獻[10]中的學習函數擬合出近似的球頭銑刀銑削力表達式。訓練過程，如圖3所示。

圖3 Kriging模型訓練過程Fig.3 The Cutter Engagement in Slot and Semi-Slot Machining

以球頭銑刀公稱螺旋角、銑刀直徑、銑削深度、每齒進給量為隨機變量，其分布情況，如表1所示。以銑削力為研究對象，根據分布形式，利用解析方法計算出50組X，Y，Z方向的銑削力的響應值，并用這50組數據作為設計樣本空間對Kriging模型進行訓練。Kriging模型中的參數θ采用極大似然方法估計。學習函數選用文獻[10]中的U函數，其表達式，如式（6）所示。

式中：μG（x）—Kriging 預測值；σG（x）—Kriging 方差。

表1 隨機變量表Tab.1 Random Parameters Table

根據表1中隨機參數的分布，基于Monte Carlo模擬，對各隨機變量進行抽樣（105次），將其代入前面建立的球頭銑刀銑削力計算模型中計算；同時將其帶入訓練好的Kriging預測模型中進行預測。將球頭銑刀銑削力計算模型中的結果和Kriging模型預測的結果用MATLAB軟件處理并分別繪制概率密度曲線。

X、Y、Z軸方向的銑削力的概率密度函數圖像，如圖4～圖6所示。其中實線為用計算模型預測的銑削力概率密度分布曲線，虛線為用Kriging模型預測得到的銑削力概率密度分布曲線。圖中可以看出，兩種方法的各方向的銑削力基本服從正態分布。從各方向分力的概率密度曲線可以看出Kriging模型預測結果與計算模型預測結果基本一致，可以表明提出的方法適用于預測球頭銑刀銑削力的分布模型，Kriging模型是計算模型的一個驗證，其表達式簡單，方便計算，為后續分析球頭銑刀銑削精度進行可靠性分析等打下基礎。

圖4 X軸方向銑削力概率密度曲線Fig.4 PDF of X Milling Force

圖5 Y軸方向銑削力概率密度曲線Fig.5 PDF of Y Milling Force

圖6 Z軸方向銑削力概率密度曲線Fig.6 PDF of Z Milling Force

5 結論

基于Kriging和Monte Carlo提出一種球頭銑刀銑削力概率分布預測和分析方法。與傳統方法相比，提出的方法考慮了銑刀的形狀和加工參數等隨機參數對銑刀銑削力的影響，更加符合工程實際狀況。

該方法發揮了Kriging方法需要樣本少，擬合收斂快的優勢，可以在快速得到銑削力的概率分布特性。上述特性對提高加工精度、研究刀具使用壽命等具有重要參考價值及指導意義，為后續分析球頭銑刀銑削精度提供基礎。并且提出的方法具有深入研究的價值，適合應用于其他領域。