樊隆祥,周啟武,彭東林,鄭 永,徐 是
(重慶理工大學 機械檢測技術(shù)與裝備教育部工程中心, 重慶 400054)
機床熱誤差占已知加工誤差的40%~70%[1],嚴重影響機床加工精度。對于數(shù)控滾齒機,熱誤差是一個主要誤差源[2],有必要建立有效的補償方式。建立一個準確的數(shù)控滾齒機熱誤差模型是一個重要環(huán)節(jié),好的建模結(jié)果使得補償更為有效,更利于提升滾齒加工精度。
就數(shù)控機床熱誤差模型而言,李彬等[3]利用遺傳算法自適應(yīng)的全局優(yōu)化搜索特性和小波神經(jīng)網(wǎng)絡(luò)的時頻局部特性建立了數(shù)控機床的熱誤差補償模型,提高了熱誤差補償?shù)木燃靶?;Liu等[4]提出一種改進的灰狼優(yōu)化算法(IGWO)和廣義回歸神經(jīng)網(wǎng)絡(luò)相結(jié)合(GRNN)的滾齒機熱誤差模型,將模糊聚類和平均影響值相結(jié)合,從而減少溫度變量間的耦合。本文以改造的數(shù)控滾齒機作為對象,進行熱-結(jié)構(gòu)仿真分析,基于仿真數(shù)據(jù),選擇多個對機床主軸熱誤差影響大的部件的溫度數(shù)據(jù),對溫度變量進行模糊聚類分析,然后通過與主軸熱位移相關(guān)性比較篩選得到滾齒機熱誤差建模的熱關(guān)鍵點,同時減少了溫度傳感器數(shù)量,最后基于貝葉斯網(wǎng)絡(luò)建立數(shù)控滾齒機熱誤差的數(shù)學模型。
該數(shù)控滾齒機主要由立柱、直驅(qū)滾刀軸系統(tǒng)、工作臺、進給系統(tǒng)、床身等組成,其結(jié)構(gòu)如圖1所示。三維模型主要為刀架部分模型,對其余結(jié)構(gòu)做了簡化,如圖2所示。
圖1 滾齒機組成 圖2 滾齒機三維模型
數(shù)控滾齒機床熱源主要來自滾刀切削熱、電機發(fā)熱及軸承的發(fā)熱[5]。
2.1.1切削熱
切削過程中,切削熱大部分由切屑帶走,取5%的切削熱量傳至滾刀[6],切削熱的計算公式:
Q=Pν
(1)
式中:Q為熱功率;P為切削力;ν為主軸轉(zhuǎn)速。
滾刀最大力矩計算[7]:
(2)
P=2Mmax/d
(3)
式中:m為法向模數(shù);S為軸向進給量;T為吃刀深度;ν0為主軸轉(zhuǎn)速;z為齒輪齒數(shù);K材為齒輪材料修正系數(shù),根據(jù)加工材料選取相應(yīng)的系數(shù);K硬為工件硬度修正系數(shù);K螺為螺旋角修正系數(shù);d為滾刀外徑。
2.1.2軸承熱
軸承的摩擦發(fā)熱是滾齒機的一個主要熱源,本結(jié)構(gòu)中滑動軸承是主要熱源,其摩擦熱計算為[8]
Q1=fpv
(4)
式中:Q1為摩擦升熱;f為摩擦因數(shù);p為軸承載荷;v為軸頸圓周速度。
該數(shù)控滾齒機機床中滾動軸承為推力球軸承,其發(fā)熱功率計算[9]:
q=1.047×10-4M·n
(5)
M=M0+M1
(6)
式中:q為功率;M為軸承摩擦力矩;n為主軸轉(zhuǎn)速;M0為載荷摩擦力矩;M1為粘性摩擦力矩。
(7)
M1=f1p1dm
(8)
式中:f0為軸承型號與潤滑方式的相關(guān)常數(shù);f1為軸承型號和負載的相關(guān)常數(shù);ν為潤滑劑運動黏度;p1為軸承摩擦的力矩計算載荷;dm為軸承中徑。
2.1.3主軸直驅(qū)電機發(fā)熱
主軸電機的發(fā)熱主要為定子銅耗和鐵芯損耗,其發(fā)熱表示為[10]:
Qd=P1(1-η)
(9)
式中:Qd為電機生熱;P1為電機輸入功率;η為電機效率。
滾齒加工中,刀架主軸、電機轉(zhuǎn)子轉(zhuǎn)動與空氣及滾刀的切削液冷卻均為強迫對流換熱。強迫對流換熱系數(shù)公式[6]:
(10)
(11)
(12)
式中:α為對流換熱系數(shù);Nu、Re、Pr為努塞爾特數(shù)、雷諾數(shù)、普朗特數(shù);λ為流體導(dǎo)熱系數(shù);ds為主軸當量直徑;ω為主軸角速度;νf為流體運動黏度。
電機水冷為管內(nèi)強迫對流,對流換熱公式[11]:
(13)
(14)
式中:h1為換熱系數(shù);D為幾何特征定型尺寸;u為流體特征速度;ν為水的運動黏度。
采用以下參數(shù)作為仿真條件,分析滾齒機受熱與變形。齒輪參數(shù):m=6,z=20,β=10°,材料為45鋼;滾刀參數(shù):d=105 mm,z1=1,軸向進給量S=2 mm/r,n=160 r/min,吃刀深度T=1 mm。用于機床仿真的主要材料參數(shù)如表1所示,表2為根據(jù)加工條件進行的熱源強度值。
表1 機床仿真主要零件材料參數(shù)
表2 機床熱源強度
環(huán)境溫度設(shè)定為20 ℃,分析熱邊界條件,表3為刀架部分的強迫對流換熱系數(shù)。
表3 機床換熱系數(shù)
利用有限元計算軟件完成穩(wěn)態(tài)熱-結(jié)構(gòu)仿真。圖3為模型的四面體網(wǎng)格劃分和三維坐標的標識,研究主要圍繞機床刀架主軸熱誤差展開,因此對主軸部分做了網(wǎng)格細化,主軸細化程度對主軸熱變形計算結(jié)果無明顯差別,模型的單元數(shù)為397 046,節(jié)點數(shù)為632 875。
以計算出的熱功率和對流系數(shù)作為熱分析條件進行穩(wěn)態(tài)熱仿真,得到如圖4所示的機床達到熱平衡狀態(tài)時的溫度分布圖。結(jié)果顯示,水冷雖然帶走了電機定子大部分的熱量,但仍有一部分滯留在電機內(nèi)部,電機轉(zhuǎn)子溫升明顯,溫度達到41.6 ℃;滑動軸承部分的溫度上升也較多,尾座軸承達到33 ℃;滾刀冷卻效果較好,溫度變化較小。
以穩(wěn)態(tài)熱分析結(jié)果作為結(jié)構(gòu)仿真分析的條件,并設(shè)置機床底座的底面為固定約束,得到圖5—8所示的機床熱位移分布情況,其結(jié)果的趨勢與溫度分布趨勢相當,靠近軸承和電機內(nèi)部的區(qū)域因溫度較高導(dǎo)致熱位移較大,最大的變形發(fā)生在電機轉(zhuǎn)子部位,約92 μm;滾刀主軸部分的變形較小,約35 μm,同時圖6—8的結(jié)果顯示滾刀主軸部分主要表現(xiàn)為X和Y2個方向的熱位移,Z方向則相對較小。
圖3 網(wǎng)格劃分 圖4 穩(wěn)態(tài)溫度場
圖5 穩(wěn)態(tài)熱位移 圖6 X向熱位移
圖7 Y向熱位移 圖8 Z向熱位移
為得到機床溫升和主軸熱位移變化情況,進行滾齒機瞬態(tài)熱-結(jié)構(gòu)仿真,提取機床熱源附近7個重要部件溫度(分別為T1-T7,其中滑動軸承Ⅱ位于主軸中部)和滾刀主軸熱位移的數(shù)據(jù)如圖9、圖10所示。隨時間推移,滾齒機各部分溫度持續(xù)上升,滾刀主軸熱位移也隨溫度的上升而增大,與穩(wěn)態(tài)仿真結(jié)果分析一樣,主軸的熱位移主要為X和Y2個方向。
圖9 滾齒機各部分溫升曲線
圖10 滾刀主軸熱位移曲線
機床熱-結(jié)構(gòu)的仿真數(shù)據(jù)一定程度上反映了機床的誤差情況,因此提取瞬態(tài)仿真的重要溫度變量數(shù)據(jù)來進行溫度優(yōu)化,從而獲得機床熱關(guān)鍵點,方法采用模糊聚類[12],首先構(gòu)建模糊等價矩陣,計算各溫度變量間的相似系數(shù)rij組合成n階相似系數(shù)矩陣(rij)n×n,記為R。相似系數(shù)rij采用皮爾遜相關(guān)系數(shù)計算方式。相關(guān)系數(shù)rxy公式[13]:
(15)
式中:x、y為變量。
采用數(shù)學平方法求出相似系數(shù)矩陣R的傳遞閉包R,運算如下:
R→R2→R4→…
(16)
表4所示為溫度變量分類情況及與熱位移相關(guān)系數(shù)大小。依據(jù)溫度變量Ti與熱位移E的相關(guān)性系數(shù)大小得到機床熱關(guān)鍵點,同時電機主軸的相關(guān)性相對較小,因此以滾齒機刀架、刀架尾座軸承、滑動軸承Ⅱ作為滾齒機熱誤差建模的熱關(guān)鍵點。
表4 溫度變量分類
貝葉斯網(wǎng)絡(luò)是一種模擬人在推理過程中對于因果關(guān)系的不確定性的處理模型。一組變量X={X1,X2,…,Xn}的貝葉斯網(wǎng)絡(luò)包含由組中所有變量節(jié)點共同組成的網(wǎng)絡(luò)S,如圖11所示。
圖11 網(wǎng)絡(luò)結(jié)構(gòu)S
網(wǎng)絡(luò)結(jié)構(gòu)S是由各個節(jié)點之間共同構(gòu)成的有向無環(huán)圖,節(jié)點之間采用單箭頭進行連接則表示2個節(jié)點變量之間的依賴關(guān)系,箭頭一端為“因”,一端為“果”,節(jié)點間的不確定性通過概率P量化表示[14]。簡言之,把某個研究系統(tǒng)中涉及的隨機變量,根據(jù)是否條件獨立繪制在一個有向圖中,就形成了貝葉斯網(wǎng)絡(luò)。
圖11中節(jié)點X1影響到節(jié)點X2,用從X1指向X2的箭頭建立有向弧(X1,X1),權(quán)值用條件概率P(X2|X1)表示,同時,對于隨機變量X的聯(lián)合概率分布由各自的局部條件概率分布相乘得到:
(17)
貝葉斯網(wǎng)絡(luò)的熱誤差建模過程主要分為3步:貝葉斯網(wǎng)絡(luò)先驗網(wǎng)絡(luò)構(gòu)造、樣本數(shù)據(jù)的網(wǎng)絡(luò)學習以及貝葉斯網(wǎng)絡(luò)推理的熱誤差預(yù)測,如圖12所示[15]。
圖12 貝葉斯網(wǎng)絡(luò)的機床熱誤差建模過程
根據(jù)貝葉斯網(wǎng)絡(luò),由刀架尾座軸承、刀架、滑動軸承Ⅱ作為熱關(guān)鍵點與滾刀主軸熱誤差一起構(gòu)成網(wǎng)絡(luò)的節(jié)點集合,為簡化網(wǎng)絡(luò)和降低建模的難度,假設(shè)溫度節(jié)點間是相互獨立的,構(gòu)造的3個溫度點下的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)如圖13所示。
圖13 3個溫度點下的貝葉斯網(wǎng)絡(luò)
實驗采集機床熱關(guān)鍵點的溫升和滾刀主軸熱誤差的數(shù)據(jù),以滾刀主軸X方向的熱誤差作為實驗對象進行空載實驗,實驗裝置如圖14所示。滾刀主軸轉(zhuǎn)速160 r/min,環(huán)境溫度20 ℃,溫度測量采用DS18B20溫度傳感器。
根據(jù)4.2節(jié)貝葉斯網(wǎng)絡(luò)的理論建立熱誤差模型。首先對實驗數(shù)據(jù)變量域進行確定,數(shù)據(jù)的特征統(tǒng)計如表5所示,以尾座滑動軸承溫度T1和熱誤差E為例。ΔT1為各溫度變量T1相鄰采樣點的差值(單位:℃),ΔE為熱誤差數(shù)據(jù)E相鄰采樣點的差值(單位:μm)。
圖14 數(shù)控滾齒機溫度及熱誤差采集實驗裝置
表5 數(shù)據(jù)特征統(tǒng)計
根據(jù)表5,將ΔT1和ΔE的值域均勻的劃為10個區(qū)間,每個區(qū)間對應(yīng)于一個狀態(tài)域,完成數(shù)據(jù)的離散化,如表6所示。
表6 變量狀態(tài)域劃分
確定變量對應(yīng)的變量域后,得出ΔT1和ΔE獨立分布概率(整個數(shù)據(jù)中該變量落在不同狀態(tài)域的概率),如表7所示。
表7 變量的獨立分布概率
得出ΔT1和ΔE之間的條件概率P(ΔT1|ΔE),如表8所示。
表8 條件概率
計算max{Pk(ΔEj|ΔT1m,ΔT4n,ΔT5p)},即變量ΔTi(i=1,4,5)的狀態(tài)域確定的條件下,變量ΔEj最大概率處于1~10的某個狀態(tài)域(k為變量數(shù)據(jù)組號;j,m,n,p為狀態(tài)域編號),取此時j所對應(yīng)的誤差變量的狀態(tài)域作為本組ΔEk的預(yù)測區(qū)間,取狀態(tài)域的中值作為ΔEk的預(yù)測值,得出:
Ek+1=Ek+ΔEk
(18)
式中:Ek+1為第k+1組的熱誤差預(yù)測值;ΔEk為第k+1組和第k組熱誤差差值的預(yù)測值。
建立貝葉斯網(wǎng)絡(luò)熱誤差預(yù)測模型后,導(dǎo)入實驗數(shù)據(jù)進行熱誤差結(jié)果預(yù)測。數(shù)控滾齒機熱關(guān)鍵點溫度數(shù)據(jù)如圖15所示,根據(jù)貝葉斯網(wǎng)絡(luò)建立的滾齒機熱誤差模型,其熱誤差實驗數(shù)據(jù)與預(yù)測模型的結(jié)果比較如圖16所示。熱誤差預(yù)測模型的變化趨勢與實際熱誤差的變化趨勢較為一致,同時在7 h內(nèi)二者殘差最大為3.5 μm,預(yù)測結(jié)果較為準確。
圖15 滾齒機及熱關(guān)鍵點溫度曲線
圖16 滾齒機熱誤差實驗數(shù)據(jù)與貝葉斯網(wǎng)絡(luò)預(yù)測值曲線
1) 基于傳熱理論完成了數(shù)控滾齒機整機的熱特性分析,明確了機床的主要熱源,進行機床熱-結(jié)構(gòu)仿真,分析了機床溫度及熱位移情況。
2) 采用模糊聚類對溫度變量的優(yōu)化,消除了溫度變量間的耦合,使關(guān)鍵溫度測點從7個減小為3個,減少了實驗所需的傳感器,降低了建模的難度。
3) 采用貝葉斯網(wǎng)絡(luò)建立的數(shù)控滾齒機熱誤差模型對熱誤差進行預(yù)測,并把模型預(yù)測誤差數(shù)據(jù)與實測誤差數(shù)據(jù)進行比較,殘差最大為3.5 μm,表明所提熱誤差模型是有效的,預(yù)測精度較高。