基于有限元及貝葉斯網(wǎng)絡(luò)的數(shù)控滾齒機(jī)熱誤差建模研究

樊隆祥，周啟武，彭東林，鄭 永，徐 是

(重慶理工大學(xué) 機(jī)械檢測技術(shù)與裝備教育部工程中心， 重慶 400054)

0 引言

機(jī)床熱誤差占已知加工誤差的40%～70%[1]，嚴(yán)重影響機(jī)床加工精度。對于數(shù)控滾齒機(jī)，熱誤差是一個主要誤差源[2]，有必要建立有效的補(bǔ)償方式。建立一個準(zhǔn)確的數(shù)控滾齒機(jī)熱誤差模型是一個重要環(huán)節(jié)，好的建模結(jié)果使得補(bǔ)償更為有效，更利于提升滾齒加工精度。

就數(shù)控機(jī)床熱誤差模型而言，李彬等[3]利用遺傳算法自適應(yīng)的全局優(yōu)化搜索特性和小波神經(jīng)網(wǎng)絡(luò)的時頻局部特性建立了數(shù)控機(jī)床的熱誤差補(bǔ)償模型，提高了熱誤差補(bǔ)償?shù)木燃靶剩籐iu等[4]提出一種改進(jìn)的灰狼優(yōu)化算法(IGWO)和廣義回歸神經(jīng)網(wǎng)絡(luò)相結(jié)合(GRNN)的滾齒機(jī)熱誤差模型，將模糊聚類和平均影響值相結(jié)合，從而減少溫度變量間的耦合。本文以改造的數(shù)控滾齒機(jī)作為對象，進(jìn)行熱-結(jié)構(gòu)仿真分析，基于仿真數(shù)據(jù)，選擇多個對機(jī)床主軸熱誤差影響大的部件的溫度數(shù)據(jù)，對溫度變量進(jìn)行模糊聚類分析，然后通過與主軸熱位移相關(guān)性比較篩選得到滾齒機(jī)熱誤差建模的熱關(guān)鍵點(diǎn)，同時減少了溫度傳感器數(shù)量，最后基于貝葉斯網(wǎng)絡(luò)建立數(shù)控滾齒機(jī)熱誤差的數(shù)學(xué)模型。

1 數(shù)控滾齒機(jī)模型

該數(shù)控滾齒機(jī)主要由立柱、直驅(qū)滾刀軸系統(tǒng)、工作臺、進(jìn)給系統(tǒng)、床身等組成，其結(jié)構(gòu)如圖1所示。三維模型主要為刀架部分模型，對其余結(jié)構(gòu)做了簡化，如圖2所示。

圖1 滾齒機(jī)組成 圖2 滾齒機(jī)三維模型

2 數(shù)控滾齒機(jī)熱源和邊界條件分析

2.1 數(shù)控滾齒機(jī)熱源分析計算

數(shù)控滾齒機(jī)床熱源主要來自滾刀切削熱、電機(jī)發(fā)熱及軸承的發(fā)熱[5]。

2.1.1切削熱

切削過程中，切削熱大部分由切屑帶走，取5%的切削熱量傳至滾刀[6]，切削熱的計算公式：

Q=Pν

(1)

式中：Q為熱功率;P為切削力;ν為主軸轉(zhuǎn)速。

滾刀最大力矩計算[7]：

(2)

P=2Mmax/d

(3)

式中：m為法向模數(shù)；S為軸向進(jìn)給量；T為吃刀深度；ν0為主軸轉(zhuǎn)速；z為齒輪齒數(shù)；K材為齒輪材料修正系數(shù)，根據(jù)加工材料選取相應(yīng)的系數(shù)；K硬為工件硬度修正系數(shù)；K螺為螺旋角修正系數(shù)；d為滾刀外徑。

2.1.2軸承熱

軸承的摩擦發(fā)熱是滾齒機(jī)的一個主要熱源，本結(jié)構(gòu)中滑動軸承是主要熱源，其摩擦熱計算為[8]

Q1=fpv

(4)

式中：Q1為摩擦升熱；f為摩擦因數(shù)；p為軸承載荷；v為軸頸圓周速度。

該數(shù)控滾齒機(jī)機(jī)床中滾動軸承為推力球軸承，其發(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為軸承型號和負(fù)載的相關(guān)常數(shù)；ν為潤滑劑運(yùn)動黏度；p1為軸承摩擦的力矩計算載荷；dm為軸承中徑。

2.1.3主軸直驅(qū)電機(jī)發(fā)熱

主軸電機(jī)的發(fā)熱主要為定子銅耗和鐵芯損耗，其發(fā)熱表示為[10]：

Qd=P1(1-η)

(9)

式中：Qd為電機(jī)生熱；P1為電機(jī)輸入功率；η為電機(jī)效率。

2.2 邊界條件分析計算

滾齒加工中，刀架主軸、電機(jī)轉(zhuǎn)子轉(zhuǎn)動與空氣及滾刀的切削液冷卻均為強(qiáng)迫對流換熱。強(qiáng)迫對流換熱系數(shù)公式[6]：

(10)

(11)

(12)

式中：α為對流換熱系數(shù)；Nu、Re、Pr為努塞爾特數(shù)、雷諾數(shù)、普朗特數(shù)；λ為流體導(dǎo)熱系數(shù)；ds為主軸當(dāng)量直徑；ω為主軸角速度；νf為流體運(yùn)動黏度。

電機(jī)水冷為管內(nèi)強(qiáng)迫對流，對流換熱公式[11]：

(13)

(14)

式中：h1為換熱系數(shù)；D為幾何特征定型尺寸；u為流體特征速度；ν為水的運(yùn)動黏度。

3 數(shù)控滾齒機(jī)熱-結(jié)構(gòu)仿真

采用以下參數(shù)作為仿真條件，分析滾齒機(jī)受熱與變形。齒輪參數(shù)：m=6，z=20，β=10°，材料為45鋼；滾刀參數(shù)：d=105 mm，z1=1，軸向進(jìn)給量S=2 mm/r，n=160 r/min，吃刀深度T=1 mm。用于機(jī)床仿真的主要材料參數(shù)如表1所示，表2為根據(jù)加工條件進(jìn)行的熱源強(qiáng)度值。

表1 機(jī)床仿真主要零件材料參數(shù)

表2 機(jī)床熱源強(qiáng)度

環(huán)境溫度設(shè)定為20 ℃，分析熱邊界條件，表3為刀架部分的強(qiáng)迫對流換熱系數(shù)。

表3 機(jī)床換熱系數(shù)

利用有限元計算軟件完成穩(wěn)態(tài)熱-結(jié)構(gòu)仿真。圖3為模型的四面體網(wǎng)格劃分和三維坐標(biāo)的標(biāo)識，研究主要圍繞機(jī)床刀架主軸熱誤差展開，因此對主軸部分做了網(wǎng)格細(xì)化，主軸細(xì)化程度對主軸熱變形計算結(jié)果無明顯差別，模型的單元數(shù)為397 046，節(jié)點(diǎn)數(shù)為632 875。

以計算出的熱功率和對流系數(shù)作為熱分析條件進(jìn)行穩(wěn)態(tài)熱仿真，得到如圖4所示的機(jī)床達(dá)到熱平衡狀態(tài)時的溫度分布圖。結(jié)果顯示，水冷雖然帶走了電機(jī)定子大部分的熱量，但仍有一部分滯留在電機(jī)內(nèi)部，電機(jī)轉(zhuǎn)子溫升明顯，溫度達(dá)到41.6 ℃；滑動軸承部分的溫度上升也較多，尾座軸承達(dá)到33 ℃；滾刀冷卻效果較好，溫度變化較小。

以穩(wěn)態(tài)熱分析結(jié)果作為結(jié)構(gòu)仿真分析的條件，并設(shè)置機(jī)床底座的底面為固定約束，得到圖5—8所示的機(jī)床熱位移分布情況，其結(jié)果的趨勢與溫度分布趨勢相當(dāng)，靠近軸承和電機(jī)內(nèi)部的區(qū)域因溫度較高導(dǎo)致熱位移較大，最大的變形發(fā)生在電機(jī)轉(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向熱位移

為得到機(jī)床溫升和主軸熱位移變化情況，進(jìn)行滾齒機(jī)瞬態(tài)熱-結(jié)構(gòu)仿真，提取機(jī)床熱源附近7個重要部件溫度(分別為T1-T7，其中滑動軸承Ⅱ位于主軸中部)和滾刀主軸熱位移的數(shù)據(jù)如圖9、圖10所示。隨時間推移，滾齒機(jī)各部分溫度持續(xù)上升，滾刀主軸熱位移也隨溫度的上升而增大，與穩(wěn)態(tài)仿真結(jié)果分析一樣，主軸的熱位移主要為X和Y2個方向。

圖9 滾齒機(jī)各部分溫升曲線

圖10 滾刀主軸熱位移曲線

4 貝葉斯網(wǎng)絡(luò)的滾齒機(jī)熱誤差建模

4.1 基于模糊聚類的溫度優(yōu)化

機(jī)床熱-結(jié)構(gòu)的仿真數(shù)據(jù)一定程度上反映了機(jī)床的誤差情況，因此提取瞬態(tài)仿真的重要溫度變量數(shù)據(jù)來進(jìn)行溫度優(yōu)化，從而獲得機(jī)床熱關(guān)鍵點(diǎ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ù)學(xué)平方法求出相似系數(shù)矩陣R的傳遞閉包R，運(yùn)算如下：

R→R2→R4→…

(16)

表4所示為溫度變量分類情況及與熱位移相關(guān)系數(shù)大小。依據(jù)溫度變量Ti與熱位移E的相關(guān)性系數(shù)大小得到機(jī)床熱關(guān)鍵點(diǎn)，同時電機(jī)主軸的相關(guān)性相對較小，因此以滾齒機(jī)刀架、刀架尾座軸承、滑動軸承Ⅱ作為滾齒機(jī)熱誤差建模的熱關(guān)鍵點(diǎn)。

表4 溫度變量分類

4.2 基于貝葉斯網(wǎng)絡(luò)的滾齒機(jī)熱誤差的構(gòu)造理論

貝葉斯網(wǎng)絡(luò)是一種模擬人在推理過程中對于因果關(guān)系的不確定性的處理模型。一組變量X={X1,X2,…,Xn}的貝葉斯網(wǎng)絡(luò)包含由組中所有變量節(jié)點(diǎn)共同組成的網(wǎng)絡(luò)S，如圖11所示。

圖11 網(wǎng)絡(luò)結(jié)構(gòu)S

網(wǎng)絡(luò)結(jié)構(gòu)S是由各個節(jié)點(diǎn)之間共同構(gòu)成的有向無環(huán)圖，節(jié)點(diǎn)之間采用單箭頭進(jìn)行連接則表示2個節(jié)點(diǎn)變量之間的依賴關(guān)系，箭頭一端為“因”，一端為“果”，節(jié)點(diǎn)間的不確定性通過概率P量化表示[14]。簡言之，把某個研究系統(tǒng)中涉及的隨機(jī)變量，根據(jù)是否條件獨(dú)立繪制在一個有向圖中，就形成了貝葉斯網(wǎng)絡(luò)。

圖11中節(jié)點(diǎn)X1影響到節(jié)點(diǎn)X2，用從X1指向X2的箭頭建立有向弧(X1，X1)，權(quán)值用條件概率P(X2|X1)表示，同時，對于隨機(jī)變量X的聯(lián)合概率分布由各自的局部條件概率分布相乘得到：

(17)

貝葉斯網(wǎng)絡(luò)的熱誤差建模過程主要分為3步：貝葉斯網(wǎng)絡(luò)先驗網(wǎng)絡(luò)構(gòu)造、樣本數(shù)據(jù)的網(wǎng)絡(luò)學(xué)習(xí)以及貝葉斯網(wǎng)絡(luò)推理的熱誤差預(yù)測，如圖12所示[15]。

圖12 貝葉斯網(wǎng)絡(luò)的機(jī)床熱誤差建模過程

根據(jù)貝葉斯網(wǎng)絡(luò)，由刀架尾座軸承、刀架、滑動軸承Ⅱ作為熱關(guān)鍵點(diǎn)與滾刀主軸熱誤差一起構(gòu)成網(wǎng)絡(luò)的節(jié)點(diǎn)集合，為簡化網(wǎng)絡(luò)和降低建模的難度，假設(shè)溫度節(jié)點(diǎn)間是相互獨(dú)立的，構(gòu)造的3個溫度點(diǎn)下的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)如圖13所示。

圖13 3個溫度點(diǎn)下的貝葉斯網(wǎng)絡(luò)

4.3 機(jī)床熱誤差實驗及誤差預(yù)測模型

實驗采集機(jī)床熱關(guān)鍵點(diǎn)的溫升和滾刀主軸熱誤差的數(shù)據(jù)，以滾刀主軸X方向的熱誤差作為實驗對象進(jìn)行空載實驗，實驗裝置如圖14所示。滾刀主軸轉(zhuǎn)速160 r/min，環(huán)境溫度20 ℃，溫度測量采用DS18B20溫度傳感器。

根據(jù)4.2節(jié)貝葉斯網(wǎng)絡(luò)的理論建立熱誤差模型。首先對實驗數(shù)據(jù)變量域進(jìn)行確定，數(shù)據(jù)的特征統(tǒng)計如表5所示，以尾座滑動軸承溫度T1和熱誤差E為例。ΔT1為各溫度變量T1相鄰采樣點(diǎn)的差值(單位：℃)，ΔE為熱誤差數(shù)據(jù)E相鄰采樣點(diǎn)的差值(單位：μm)。

圖14 數(shù)控滾齒機(jī)溫度及熱誤差采集實驗裝置

表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獨(dú)立分布概率(整個數(shù)據(jù)中該變量落在不同狀態(tài)域的概率)，如表7所示。

表7 變量的獨(dú)立分布概率

得出Δ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ù)進(jìn)行熱誤差結(jié)果預(yù)測。數(shù)控滾齒機(jī)熱關(guān)鍵點(diǎn)溫度數(shù)據(jù)如圖15所示，根據(jù)貝葉斯網(wǎng)絡(luò)建立的滾齒機(jī)熱誤差模型，其熱誤差實驗數(shù)據(jù)與預(yù)測模型的結(jié)果比較如圖16所示。熱誤差預(yù)測模型的變化趨勢與實際熱誤差的變化趨勢較為一致，同時在7 h內(nèi)二者殘差最大為3.5 μm，預(yù)測結(jié)果較為準(zhǔn)確。

圖15 滾齒機(jī)及熱關(guān)鍵點(diǎn)溫度曲線

圖16 滾齒機(jī)熱誤差實驗數(shù)據(jù)與貝葉斯網(wǎng)絡(luò)預(yù)測值曲線

5 結(jié)論

1) 基于傳熱理論完成了數(shù)控滾齒機(jī)整機(jī)的熱特性分析，明確了機(jī)床的主要熱源，進(jìn)行機(jī)床熱-結(jié)構(gòu)仿真，分析了機(jī)床溫度及熱位移情況。

2) 采用模糊聚類對溫度變量的優(yōu)化，消除了溫度變量間的耦合，使關(guān)鍵溫度測點(diǎn)從7個減小為3個，減少了實驗所需的傳感器，降低了建模的難度。

3) 采用貝葉斯網(wǎng)絡(luò)建立的數(shù)控滾齒機(jī)熱誤差模型對熱誤差進(jìn)行預(yù)測，并把模型預(yù)測誤差數(shù)據(jù)與實測誤差數(shù)據(jù)進(jìn)行比較，殘差最大為3.5 μm，表明所提熱誤差模型是有效的，預(yù)測精度較高。