非線性銑削熱傳導逆問題的界面熱流辨識*

馮 勇 汪木蘭 王保升

(①南京工程學院機械工程學院，江蘇南京 211167;②先進制造技術江蘇省高校重點建設試驗室，江蘇 南京 211167)

高速切削工件時，由于切削速度快、時間短，熱量在刀具、工件及空氣中的傳遞形式主要為熱傳導。若在切削加工問題中已知工件的界面熱流密度或者界面溫度分布，就可以應用熱傳導方程來研究工件的溫度分布狀況，這一問題稱為熱傳導正問題(DHCP:direct heat conduction problem)。但通常工件的界面熱流或者界面溫度分布無法直接通過在工件表面測量來獲得。在這種情況下，需要由工件內部某些位置的溫度值來反演出界面熱流密度，這一問題稱為熱傳導逆問題(IHCP:inverse heat conduction problem)［1］。

國內外學者對高速切削加工熱傳導問題開展了大量研究。Radulescu 和Kapoor 利用分離變量法建立了銑削中刀具溫度的三維模型，并用解析法計算了熱流量［2］;Lazoglu 利用有限差分法求解了切削熱傳導模型，計算了銑削加工過程中刀具和切屑的溫度場［3］，M.R.Golbahar 通過對二維熱傳導逆問題的求解，對功能梯度材料的時變熱流進行了分析預測［4］。山東大學張士軍、劉戰(zhàn)強采用熱傳導方法對涂層刀具的導熱系數(shù)、熱流量、切削溫度場進行計算和有限元仿真［5］;安徽理工大學黃紹服、李君等人，采用熱傳導方法建立了高速銑削薄壁件溫度場的數(shù)學模型［6］;重慶大學周憶、梁錫昌等人采用熱源法和熱傳導法對超高速銑削加工中切削溫度場進行建模［7］。而結合實測溫度數(shù)據(jù)，考慮切削過程中溫度疊加的非線性效應，采用優(yōu)化算法對銑削加工中界面熱流密度進行準確辨識的研究卻極為少見。本文在建立銑削加工過程中工件的熱傳導數(shù)學模型基礎上，采用PSO 算法，通過與工件內部的測量溫度值的比較來確定被測工件的動態(tài)界面熱流密度。由于該方法只需工件內部進行溫度測量，不破壞被測工件的表面完整性，使得所計算出的界面條件比較準確，應用于熱傳導正問題的求解可以更準確得出工件的溫度分布狀況。

1 銑削熱傳導問題的數(shù)學描述

銑削加工是一個斷續(xù)加工過程，刀具中每個切削刃都經歷一個切入和切出的過程，切入工件時所產生的熱量較小，將要切出時所產生的熱量最大，切出后熱源強度降為零。銑削加工中的移動熱源形式可假設為三角形分布(圖1 所示)，熱流密度可表示為［8］:

n=0，1，2，...，t 為加工周期;t1為刀具切入時間。

采用熱傳導法求解工件在切削區(qū)的溫度場分布問題可抽象成二維傳熱問題，所建立的銑削加工過程如圖2 所示(圖中n 為主軸轉速;aw為徑向切削深度;R 為刀具半徑)。此時，工件固定不動，刀具以主軸轉速n 旋轉并以進給速度vf沿X 方向前進，忽略材料塑性變形所產生的熱量，并假定工件的比熱、密度以及熱傳導系數(shù)為常數(shù)，對于干切削過程，銑削面與空氣的對流換熱系數(shù)比較小，工件表面可看成絕熱邊界。根據(jù)熱力學第一定律，可得銑削加工的導熱微分方程為［9］:

式中:T 為工件上任一點的溫度，℃;k 為工件的熱傳導系數(shù)，W/(m·℃);ρ 為工件的密度，kg/m3;c 為比熱容，J/(kg·℃);vf為刀具沿x 方向的進給速度，m/s;τ為切削時間，s。給定邊界條件、初始條件分別如式(3)所示:

式中:q 為熱流密度，W/(m2℃);R 為刀具半徑，m;T0為初始時刻工件溫度，℃。若給定q 及T0的值來求解上述問題，稱之為熱傳導正問題。若通過某些點的溫度測量值來求得T，則稱之為熱傳導逆問題。

圖1 銑削加工中的移動熱源形式

圖2 銑削加工工作示意圖

求解熱傳導逆問題一般采用數(shù)值計算，近年來研究發(fā)現(xiàn)人工智能中的一些算法可以較好地應用于該問題的求解，常見的方法有啟發(fā)式算法和非啟發(fā)式算法。啟發(fā)式算法具有收斂速度快的優(yōu)點，但是不能保證全局收斂，使得問題陷入局部最小值，得不到最優(yōu)解。非啟發(fā)式算法包括蟻群算法、模擬退火算法、PSO 算法等，這些算法通過設置一定的約束條件，既可在求解空間內找到全局最優(yōu)解又可保證求解精度［10］。

2 基于PSO 算法的熱流密度辨識方法

2.1 辨識原理

采用PSO 算法開展銑削加工的動態(tài)熱流密度辨識研究的流程圖如圖3 所示，其基本分析步驟如下:

(1)測出工件內部某一位置處的切削溫度T(t)隨時間變化曲線，并設定界面熱流的初值q0。

(2)應用熱傳導正問題求解方法得出工件任一位置的溫度分布情況及工件內部A 位置處的溫度場T0(t)。

(3)由于熱流密度初值并不準確，因此T(t)與T0(t)存在著一定誤差，實際問題中，為提高求解精度往往采用測量多個點的溫度來評定綜合偏差，其表達式如式(4)所示。

(4)按照PSO 算法的求解思路，通過使不同時刻的綜合偏差達到預定的極小值，即可求得動態(tài)的熱流密度變化曲線。

圖3 辨識方法流程圖

2.2 算例計算

為驗證上述算法的可行性，并進一步研究不同參數(shù)及慣性權重值對尋優(yōu)效果的影響，對式(5)所給定的算例開展熱流密度辨識研究。

該算例中熱流密度q 未知，其他參數(shù)為已知:熱傳導系數(shù)k=47.7 W/(m ·K)，熱擴散率α=0.082 34 cm2/s，對流換熱系數(shù)h=20 W/(m2℃)，x1=0.016 m，δ=0.03 m。將上述熱傳導逆問題求解轉化為關于綜合偏差的優(yōu)化問題，選取測量點位置的坐標值分別為:x=1、3、5、7、9 mm，y=2、4 mm，總計10 個位置點，粒子的適應度方程如式(6)所示。求解過程中設定粒子群的總數(shù)為50，待求熱流密度的搜索范圍為(100 000，1 000 000)，速度系數(shù)取為0.1。隨機參數(shù)r1和r2為［0，1］間的隨機數(shù)，學習因子c1和c2取為2。最大搜索次數(shù)為100;最優(yōu)適應值小于100;初始權重為0.9，最終權重為0.4。

表1 t=1 s 粒子的最優(yōu)適應值的變化

首先，假定t=1 s 時刻的真實熱流密度q=335 200 W/(m2℃)，通過算例所給出的熱傳導模型，計算t=1 s 時刻10 個被測點的溫度值，…，作為真實參考溫度;其次，任選一個熱流密度值代入到熱傳導數(shù)學模型中，求得被測點溫度的計算值為T1，T2，…，T10;最后，按照上述熱流密度辨識方法的步驟對該問題進行優(yōu)化運算。計算過程中，粒子最優(yōu)位置的適應值變化情況如下表1 所示，從表中可以看出，隨著搜索次數(shù)的增加，粒子群最優(yōu)位置的適應值逐漸下降，熱流密度的計算值與實際值之間的誤差越來越小，在粒子的迭代次數(shù)超過60 以后，粒子的最優(yōu)適應值滿足要求。

圖4 不同測點個數(shù)n 計算得的熱流相對誤差曲線

圖4 為選取測量點個數(shù)為10 和20 時，熱流的相對偏差變化曲線，從中可以看出當選測點個數(shù)n=20時的相對誤差比n=10 時小，迭代效果更好且粒子的收斂能力也很好。因此，在條件允許的情況下，選擇更多的測量點能夠減小迭代過程中粒子最優(yōu)位置與目標位置間的偏差量。

3 銑削溫度測量

圖5 測溫實驗原理

3.1 實驗原理

待測點溫度測量裝置原理如圖5 所示，系統(tǒng)的坐標原點位于圖示工件的右頂端線的中心位置，X、Z 軸方向如圖所示，Y 軸由右手法則確定。將工件剖分成兩塊，工件剖面處粘貼絕緣片(云母片)，將標準熱電偶置于剖分后的工件之間，并使其略低于工件頂面1.0 mm 左右，然后將兩塊工件在工作臺上對齊夾緊。切削過程中，當熱電偶節(jié)點感受到被測點溫度時毫伏計輸出相應的電壓值，通過比對電壓值與溫度的關系曲線就可以獲得工件上任意位置點的溫度變化曲線［11-12］。實驗中選用鎳鉻-鎳硅組成的K 型標準熱電偶，該熱電偶具有測量范圍廣(0～1 300 ℃)、反應靈敏、精度高等特點。

3.2 實驗系統(tǒng)構建

3.2.1 試驗條件

(1)機床:立式加工中心，型號MCV850。

(2)工件材料及刀具:工件采用AISI1045 鋼;試件的規(guī)格為80 mm×80 mm×30 mm;刀具采用硬質合金雙刃銑刀，直徑為12 mm。

(3)測試裝置:由熱電偶回路、信號放大裝置、采集卡、計算機所組成，應用MFC AppWizrd 建立基于對話框的數(shù)據(jù)采集界面，采樣頻率為6 000 Hz，最終測試裝置如圖6 所示。K 型標準熱電偶選用直徑1 mm。

(4)環(huán)境溫度:10 ℃。

3.2.2 實驗方案

主軸轉速為3 000 r/min，軸向切削深度為1.5 mm，進給速度為100 mm/min;銑削方式:順銑干式切削;環(huán)境溫度:10 ℃;被檢測點位置的橫坐標分別為4、5、6、8、9、10、12、15、20、25，縱坐標分別為0、5，共計20個坐標位置。

圖6 溫度采集系統(tǒng)

4 銑削熱流密度辨識

采用前述算例分析中的計算參數(shù)，在所建立的銑削熱傳導模型基礎上，按照PSO 算法的求解思路，將該熱傳導逆問題的求解轉化為誤差的優(yōu)化問題。計算分析所選用工件和刀具材料的物理參數(shù)如表2 所示。

表2 工件和刀具的材料參數(shù)

為研究銑削加工過程中的動態(tài)溫度分布情況，必須計算不同時刻的界面熱流值。為此，結合20 個被測點的動態(tài)溫度變化曲線，對不同離散時刻的熱流值開展求解。表3 所示為不同時刻粒子的最優(yōu)位置適應值，從表中可以看出，粒子群最優(yōu)位置的適應值最大為65.6385，滿足偏差分析要求。將上述離散時刻的熱流密度值按照時間變化繪制成熱流曲線，如圖11 所示。由曲線可以看出，銑削加工中不同時刻的熱流呈明顯的非線性變化，全局最大熱流密度與最小熱流密度值分別為2.856 ×106和2.823 ×106W/(m2·℃)，分別出現(xiàn)在5 s 和10 s 左右。此外，從整個熱流變化曲線來看，在加工的前期(0～10s)，界面熱流先發(fā)生了急劇的增大，而后發(fā)生了驟降;在加工的中期(10～35 s)，界面熱流趨于穩(wěn)定變化階段;而在加工的后期(35～48 s)，界面熱流又發(fā)生了類似于加工前期的變化情形。上述變化趨勢表明界面熱流受到物理邊界的影響較大。

表3 典型時刻粒子的全局最優(yōu)適應度

圖7 界面熱流變化曲線

5 結語

本文基于高速銑削機理分析，建立了熱傳導逆問題求解的優(yōu)化模型，闡述了應用PSO 算法求解熱傳導逆問題的求解流程，并通過所述熱流密度辨識方法，結合實驗數(shù)據(jù)對AISI 1045 鋼在銑削過程中的熱流密度變化情況進行了分析。分析結果基本滿足了銑削加工熱傳導問題求解的誤差要求，為準確研究銑削加工的溫度分布提供了理論依據(jù)。

(1)單次切削過程中，界面熱流分布呈非線性變化，全局最大熱流與最小熱流值分別為2.856 ×106和2.823 ×106W/(m2·℃)。

(2)整個優(yōu)化過程中，粒子群最優(yōu)位置的適應值最大為65.638 5，界面熱流的計算值與實際值之間的偏差僅為2.635%，滿足了銑削加工熱傳導問題求解的誤差要求。

(3)從整個加工過程的界面熱流變化曲線分析來看，界面熱流的變化存在3 個典型階段，其中前期、后期兩個階段的界面熱流變化比較顯著，受邊界條件的影響較大。

［1］錢煒祺，周宇，何開鋒，等.非線性熱傳導逆問題的表面熱流辨識方法［J］.空氣動力學報，2012，30(2):145 -150.

［2］Radulescu R，Kapoor S G.An analytical model for prediction of tool temperature fields during continuous and interrupted cutting［J］.Journal of Engineering for industry，1994，116:135 -143.

［3］Ismail Lazoglu，Yusuf Altintas.Prediction of tool and chip temperature in continuous and interrupted machining ［J］.International Journal of Machine Tools ＆ Manufacture，2002，42(9):1011 -1022.

［4］Haghighi M R G，Eghtesad M，Malekzadeh P，et al.Two -dimensional inverse heat transfer analysis of functionally graded materials in estimating time -dependent surface heat flux［J］.Numerical Heat Transfer，2008，54:744 -762.

［5］張士軍.涂層刀具切削溫度及其測試技術研究［D］.山東:山東大學，2009.

［6］黃紹服，李君，何慶，等.薄壁件高速銑削溫度場研究［J］.安徽理工大學學報:自然科學版，2006，26(1):14 -17.

［7］周憶，梁錫昌，陳斌.超高速銑削加工的溫度場計算及生產應用［J］.中國機械工程，2003，14(13):1158 -1161.

［8］范晶晶.高速加工切削溫度軟測量技術及其試驗系統(tǒng)［D］.常州:常州大學，2013.

［9］陳明，袁人煒，凡孝勇，等.三維有限元分析在高速銑削溫度研究中的應用［J］.機械工程學報，2002，38(7):155 -161.

［10］錢煒祺，何開鋒，桂業(yè)偉，等.非穩(wěn)態(tài)表面熱流反演算法研究［J］.空氣動力學學報，2010，28(2):155 -161.

［11］馮勇，汪木蘭，王保升.高速切削熱及溫度預測研究進展［J］.機械設計與制造，2012，255(5):266 -268.

［12］Feng Yong，Wang Mulan，Wang Baosheng，et al.Influence of the process variables on the temperature distribution in AISI 1045 turning proceedings［A］.Advanced Materials Research，2012，578:1364-1368.