杜宏洋，陶 濤，侯瑞生，顏宗卓，姜歌東，梅雪松

(1.西安交通大學(xué) 機械工程學(xué)院，西安 710049；2.機械制造系統(tǒng)工程國家重點實驗室(西安交通大學(xué))，西安 710049；3.智能機器人陜西省重點實驗室(西安交通大學(xué))，西安 710049；4. 河北工程大學(xué) 機械與裝備工程學(xué)院，河北 邯鄲 056038)

機床誤差主要包括幾何誤差、熱誤差、切削力誤差和控制誤差等。 對于精密加工，熱誤差可以占到總誤差的40%～70%[1-2]。 主軸作為機床的關(guān)鍵零部件，也是機床的主要熱源之一。 研究減小主軸熱誤差的方法對提高機床加工精度具有十分重要的意義。 研究表明，熱誤差補償是一種經(jīng)濟、有效的提高機床精度的方法，而建立高精度、高魯棒性的熱誤差模型是熱誤差補償技術(shù)中十分關(guān)鍵的環(huán)節(jié)，很大程度上決定了最終的補償效果。 廣大學(xué)者針對主軸熱誤差建模進行了大量卓有成效的研究，提出了多種建模算法，包括多元線性回歸[3-5]、神經(jīng)網(wǎng)絡(luò)[6-9]、支持向量機[10-12]和灰色系統(tǒng)[13-14]等。 但是，模型的精度和魯棒性仍需進一步提高以滿足工廠的實際使用需求。 Yang等[15]指出，前述建模方法魯棒性差的主要原因為機床熱變形是整體溫度場的作用結(jié)果，利用少量離散的測溫點建立熱變形的預(yù)測模型時會丟失溫度信息，產(chǎn)生圖1所示的“偽滯后”現(xiàn)象，即距離熱源較近位置的溫度變化超前于熱變形，而距離熱源較遠位置的溫度變化滯后于熱變形。 因此，當選取的溫度測點位置不合適時，熱變形ΔL(τ)和溫度t(τ)之間的關(guān)系曲線呈現(xiàn)較強的非線性，且在升降溫時不重合，采用ΔL(τ)=f1(t1(τ))或ΔL(τ)=f2(t2(τ))的建模方式難以同時描述升降溫曲線。 楊建國等[16]針對熱變形偽滯后問題指出，對于主軸單端受熱情況，在靠近熱源大概x=0.4L處，主軸熱變形ΔL(τ)和溫度t(τ)呈近似的線性關(guān)系且升降溫曲線基本重合。 但是，該結(jié)論在實際應(yīng)用時存在主軸上無法布置溫度傳感器，主軸箱上與主軸0.4L處等效的溫度點不易尋找且受結(jié)構(gòu)限制也可能無法布置溫度傳感器的問題。 目前，多采用時間序列建模算法解決機床熱變形偽滯后問題，將溫度和熱變形的歷史數(shù)據(jù)作為模型輸入，彌補以離散測溫點替代整體溫度場進行建模時造成的信息缺失[17-19]。 此外，動態(tài)神經(jīng)網(wǎng)絡(luò)因包含延遲或反饋環(huán)節(jié)，也常被用來解決熱變形偽滯后問題[20]。 時間序列算法和動態(tài)神經(jīng)網(wǎng)絡(luò)的應(yīng)用對熱誤差模型的精度和魯棒性的提高起到了一定的效果，但是也引入了新的問題。 如在應(yīng)用時間序列建模時沒有統(tǒng)一的定階方法，并且模型系數(shù)完全由試驗數(shù)據(jù)確定，不具備物理意義，而動態(tài)神經(jīng)網(wǎng)絡(luò)需要設(shè)計更符合偽滯后效應(yīng)的反饋單元結(jié)構(gòu)[21]。 Fraser等[22]指出，熱誤差模型的數(shù)學(xué)形式不具備物理意義也是影響模型精度和魯棒性的重要因素。 Hou等[23]和顏宗卓等[24]針對該問題，提出了將理論建模與經(jīng)驗建模相結(jié)合的方法：首先，將主軸箱熱傳導(dǎo)簡化為半無限大平板傳熱，計算得到其溫度場和熱變形的解析表達式；然后，通過一定的數(shù)學(xué)簡化得到主軸箱熱變形的模型形式；最后，利用實驗獲取模型系數(shù)。 這種建模方法提高了模型的物理支持，但是推導(dǎo)過程中的數(shù)學(xué)簡化在一定程度上降低了模型的精度。

為了降低偽滯后效應(yīng)對熱誤差模型精度和魯棒性的影響，同時為了提高模型的物理意義，本文將主軸簡化為一維桿件，從傳熱學(xué)的角度推導(dǎo)出一維桿在變化熱源條件下熱變形的一階自回歸模型。 通過理論推導(dǎo)確定了自回歸模型階數(shù)為一階，并根據(jù)模型系數(shù)的表達式確定了系數(shù)與主軸物理特性、自回歸模型時間間隔以及熱源條件的關(guān)系。 整個推導(dǎo)過程不存在數(shù)學(xué)簡化，僅僅是進行形式變換，因此模型精度更高。

1 主軸熱變形計算及一階自回歸模型推導(dǎo)

1.1 主軸溫度場和熱變形理論求解

將主軸簡化為圖1(a)所示的一維桿結(jié)構(gòu)，在左端施加熱源，通過圓柱面進行對流散熱，其導(dǎo)熱微分方程和定解條件為：

(1)

式中：t(x,τ)為一維桿溫度分布函數(shù)(℃)，x為一維桿軸向坐標(m)，τ為時間(s)，h為表面?zhèn)鳠嵯禂?shù)(W/(m2·K))，λ為導(dǎo)熱系數(shù)(W/(m·K))，R為一維桿半徑(m)，tam為環(huán)境溫度(℃)，ρ為一維桿密度(kg/m3)，c為比熱容(J/(kg·℃))，t(x,0)為初始時刻一維桿的溫度場分布(℃)，f(x)為任意函數(shù)，L為一維桿長度(m)，q為熱流密度(W/m2)。

引入過余溫度θ(x,τ):

θ(x,τ)=t(x,τ)-tam

(2)

并且設(shè)：

a=λ/(ρc),b=2h/(ρcR),θ0(x)=f(x)-tam

(3)

將式(2)和式(3)代入式(1)可得：

(4)

對式(4)進行求解可得溫度場分布：

(5)

式中：

(6)

(7)

從式(5)中可以看出，當初始過余溫度θ0(x)=0時，一維桿位置x在τ時刻的溫升與熱流密度q成正比。 此推論在后續(xù)將被用于估計不同轉(zhuǎn)速下主軸受熱的熱流密度之間的關(guān)系。 一維桿的熱變形為

(8)

式中：ΔL(τ)為一維桿的軸向熱膨脹量(m)，α為熱膨脹系數(shù)(1/K)。

1.2 主軸熱變形一階自回歸模型推導(dǎo)

(9)

2)在mΔτ時刻一維桿左端熱流密度變?yōu)閝2，此時視為第二階段的初始時刻，初始溫度為

(10)

因此，第二階段熱變形：

(11)

由上述推導(dǎo)可知，一維桿熱變形可表述為

ΔL(nΔτ)=C1ΔL((n-1)Δτ)+C2(q)

(12)

C1=e-bΔτ

(13)

C2(q)=-qαa/(λb)(e-bΔτ-1)

(14)

由式(13)和式(14)中可看出，C1與一維桿物理特性及時間間隔Δτ相關(guān)，而一維桿物理特性通常可以認為是固定不變的，因此在時間間隔選定后，C1可視為常數(shù)。C2(q)與一維桿物理特性、時間間隔及熱流密度相關(guān)，在一維桿物理特性和時間間隔確定的條件下，系數(shù)C2與熱流密度q成正比。

主軸在實際使用時為三維傳熱，而且結(jié)構(gòu)比一維桿要復(fù)雜很多，散熱系數(shù)等邊界條件也不易精確確定。 因此本節(jié)的推導(dǎo)只是為主軸熱誤差提供一種建模的數(shù)學(xué)形式，并為不同轉(zhuǎn)速(熱流密度)下模型系數(shù)的修正提供一定的理論依據(jù)，而模型的實際系數(shù)需要根據(jù)實驗測得的數(shù)據(jù)進行確定。

2 一階自回歸模型的有限元仿真

對一階自回歸模型進行有限元仿真。 取有限元單元類型為brick 20node 226，一維桿R=0.1 m，長度L=0.8 m，導(dǎo)熱系數(shù)λ=50 W/(m·K)，密度ρ=7 850 kg/m3，比熱容c=460 J/(kg·℃)，彈性模量E=206 GPa，泊松比為0.3，熱膨脹系數(shù)為12×10-6/K，表面?zhèn)鳠嵯禂?shù)h=55 W/(m2·K)，設(shè)定環(huán)境溫度、桿的初始溫度以及參考溫度均為20 ℃。

為一維桿施加熱流密度時以工廠典型的零件加工過程為參考，即上午加工4.0 h，中午停機1.5 h，下午加工4.0 h，最后停機14.5 h完成一天的生產(chǎn)過程。 另外，機床在加工不同的零件時，加工工藝的不同導(dǎo)致主軸受熱的熱流密度不同。 因此，進行兩組計算，分別取升溫時的熱流密度為2 000 W/m2和3 000 W/m2，降溫時的熱流密度為0。 設(shè)自回歸模型時間間隔Δτ=30 s，則計算可得C1=e-bΔτ≈0.990 9，C2(2000)≈0.198 5，C2(3 000)≈0.297 7，C2(0)=0。

利用得到的一階自回歸模型系數(shù)計算一維桿熱變形，并與有限元仿真數(shù)據(jù)進行對比，結(jié)果如圖2和圖3所示。 可以看出,前述一階自回歸建模方法可有效描述一維桿在變化熱源條件下的熱變形規(guī)律，但是也存在誤差：一方面是由于在進行有限元在計算時并不是單純的一維傳熱，還存在徑向的熱傳導(dǎo)與變形；另一方面是數(shù)值計算方法本身存在的誤差。

圖2 有限元與自回歸模型計算熱變形對比(q=2 000 W/m2)

圖3 有限元與自回歸模型計算熱變形對比(q=3 000W/m2)

3 一階自回歸模型實驗驗證

3.1 T65車床Z向熱誤差一階自回歸建模

在實際使用時機械主軸采用雙端軸承支撐結(jié)構(gòu)，即為雙端受熱。 楊建國等[16]經(jīng)研究指出，一維桿在單端熱源條件下的熱變形規(guī)律可以推廣至雙端熱源。 因此，第1、2節(jié)的分析結(jié)果可應(yīng)用至實際主軸的熱變形。 本研究以海德曼T65車床為對象，進行主軸熱伸長一階自回歸建模方法的驗證。 溫度測量采用PT100溫度傳感器，量程為-20～100 ℃，精度為A級；熱變形測量采用合肥置信WD502A型電渦流傳感器，量程為2 mm，非線性誤差為0.5%，分辨率為0.2 μm。 溫度點布置見圖4(a)。 實驗時，將渦流傳感器調(diào)整至量程中部以減少非線性誤差，并采用連續(xù)啟停的方式，以60 s為一個周期。 在一個周期內(nèi)主軸運轉(zhuǎn)55 s后準停在某一固定位置5 s用于位移數(shù)據(jù)的采集。 這種測量方式可以減少主軸轉(zhuǎn)動時的跳動帶來的誤差以及檢棒自身的制造誤差。

(a)溫度傳感器布置

(b)機床X和Z向熱變形測量

本研究進行了3組實驗，以第1組實驗數(shù)據(jù)進行建模，并以另外兩組實驗數(shù)據(jù)進行模型魯棒性驗證。 3組實驗中主軸均運行12 h然后停機降溫12 h，轉(zhuǎn)速分別為2 000 、1 000 和1 500 r/min，測量結(jié)果如圖5所示。

(a)溫度曲線

(b)熱變形曲線

由圖5可看出，在轉(zhuǎn)速2 000 r/min的條件下，主軸軸承位置溫升達到10 ℃以上，會使主軸箱在豎直方向以及主軸在軸向產(chǎn)生較大的熱變形，從而影響機床的加工精度。 床身整體溫升較小，在靠近主軸箱處溫升稍大，但是因為床身尺寸較大，其熱變形對加工精度的影響不可忽略。 由熱變形曲線可知，機床在X向熱變形較小，基本控制在10 μm以內(nèi)，Z向熱變形較大，當主軸轉(zhuǎn)速為2 000 r/min時，Z向熱變形達到30 μm。 本文依據(jù)一維桿軸向熱變形規(guī)律推導(dǎo)得到一階自回歸模型，而且機床X向熱誤差基本滿足加工要求，因此本文主要針對機床Z向熱變形進行模型驗證及補償。

張老師啊，楊校長使勁咽了口唾沫說，李書記今天來學(xué)校找我談過了。他的意思是，先讓你忙一段家里的事，教學(xué)的工作先放一放。不是，不是那意思，你別誤會。我是想你如果忙起來，顧得了這頭，顧不了那頭。要不你先把家里的事忙完，再盡早回來上課。

由機床結(jié)構(gòu)分析可知，機床Z向熱變形主要受床身和主軸的熱膨脹影響，其中床身熱膨脹使傳感器架和檢棒之間的距離變大，而主軸熱膨脹使傳感器架和檢棒之間的距離變小。 在機床升溫初始階段，主軸溫升較快，其熱變形占主導(dǎo)地位；升溫一段時間后，主軸溫升變緩，床身的熱變形占主導(dǎo)地位。 這就是機床Z向熱變形曲線先迅速變小，一段時間后變平緩甚至增大的原因。 因此機床Z向熱誤差為

ΔLZ=ΔLbed-ΔLsp

(15)

式中：ΔLbed為床身Z向熱變形，ΔLsp為主軸Z向熱變形。

床身的熱膨脹變形為

ΔLbed=α*L*ΔT3

(16)

式中：L為主軸箱后部至X軸絲杠位置的床身長度，約為1.2 m，如圖4(a)中所示，α為鑄鐵的熱膨脹系數(shù)，取10×10-6/℃，ΔT3為床身中部位置的溫升，近似視為床身的平均溫升。

綜上所述，基于一階自回歸方法得到的模型為

(17)

圖6 機床Z向熱誤差建模(n=2 000 r/min)

ΔLzmlr=-2.544 4-1.411 8ΔT1+0.354 3ΔT2+26.417 6ΔT3-20.741 0ΔT4

(18)

從圖6的曲線形式可初步判斷出自回歸模型的殘差不滿足白噪聲的特性。 為進一步判斷，本文對自回歸模型的殘差進行LB檢驗，構(gòu)造的統(tǒng)計量為

(19)

表1 自回歸模型殘差的LB檢驗

由圖6可知，一階自回歸建模方法與多元線性回歸建模方法對于建模數(shù)據(jù)具有相近的擬合精度，擬合殘差基本控制在±5 μm。 但是，多元線性回歸模型依據(jù)殘差平方和最小標準對數(shù)據(jù)進行擬合，模型系數(shù)缺乏物理意義。 從式(18)中可以看出多元線性回歸模型床身的系數(shù)偏大，該模型對床身溫度的波動較為敏感。

3.2 一階自回歸模型魯棒性驗證

為驗證一階自回歸模型的魯棒性，用前述建立的模型預(yù)測主軸轉(zhuǎn)速為1 000、1 500 r/min時機床Z向熱變形，并與多元線性回歸模型進行對比。

表2 不同轉(zhuǎn)速下主軸箱后端在前100 min內(nèi)的溫升

由表2可得：

(20)

(21)

(22)

(23)

利用求得的自回歸模型預(yù)測機床在轉(zhuǎn)速為1 000 、1 500 r/min時Z向熱變形，并與式(18)的多元線性回歸模型的結(jié)果進行對比，如圖7所示。

從圖7中可以看出，在主軸轉(zhuǎn)速為1 000 r/min時，一階自回歸模型預(yù)測殘差為4 μm，預(yù)測精度為69%；多元線性回歸模型的預(yù)測殘差為6 μm，預(yù)測精度為54%。 在主軸轉(zhuǎn)速為1 500 r/min時，一階自回歸模型預(yù)測殘差為5 μm，預(yù)測精度為70%；多元線性回歸模型預(yù)測殘差為9 μm，預(yù)測精度為47%。 因此，一階自回歸模型對于非建模工況具有更好的精度保持性。

(a)n=1 000 r/min

(b)n=1 500 r/min

3.3 機床Z向熱誤差補償

利用前述分析結(jié)果對T65機床進行熱誤差補償。 本實驗所用T65機床采用西門子828D數(shù)控系統(tǒng)，具有專用的熱誤差補償接口，西門子系統(tǒng)將熱誤差補償值分為與位置無關(guān)的補償值和與位置相關(guān)的補償值，其補償原理為

ΔKZ=K0(T)+tanβ(T)×(PZ-P0)

(24)

式中：ΔKZ為位置PZ上的溫度補償值，K0(T)為與位置無關(guān)的補償值，PZ為軸的實際位置，P0為軸的參考位置，tanβ(T)為與位置相關(guān)的溫度補償?shù)南禂?shù)。

本文只考慮機床Z向與位置無關(guān)的熱誤差。 利用自主開發(fā)的溫度補償模塊采集溫度數(shù)據(jù)，并將式(17)建立的一階自回歸熱誤差模型寫入模塊內(nèi)部計算熱誤差補償量，通過將模塊并口與西門子PP 72/48D 2/2A PN外設(shè)模塊的X111口連接，將溫度和補償量數(shù)據(jù)傳輸至PLC，編寫PLC程序?qū)⒀a償量寫至Z軸的SD43900參數(shù)中。 為實現(xiàn)熱誤差補償，需將參數(shù)MD32750設(shè)定為1，即與位置無關(guān)的溫度補償功能生效，通過“垂度補償+溫度補償”參數(shù)可實時查看生效的補償值。

實施熱誤差補償后，機床主軸轉(zhuǎn)速為2 000 r/min 時Z向熱變形測量結(jié)果見圖 8。 從圖8中可以看出，機床Z向熱誤差控制在10 μm內(nèi)，滿足工廠的實際使用需求，從而驗證了一階自回歸模型的有效性。

(a)溫度曲線

(b)熱變形曲線

4 結(jié) 論

本文提出了變熱源條件下主軸軸向熱變形的一階自回歸建模方法，并通過有限元仿真及實驗驗證了模型的有效性，主要結(jié)論如下：

1)一階自回歸建模形式符合主軸在變熱源條件下的物理變形規(guī)律，建模魯棒性比多元線性回歸更高，并且不受偽滯后效應(yīng)的影響。

2)自回歸模型系數(shù)C1是與主軸物理特性及時間間隔Δτ相關(guān)的常數(shù)；C2(q)是與主軸物理特性、時間間隔及熱流密度相關(guān)的系數(shù)，當主軸轉(zhuǎn)速不同時，C2(q)需要依據(jù)熱流密度進行修正。

本文提出的主軸熱變形一階自回歸建模方法適用于機床環(huán)境溫度基本恒定的情況。 當車間環(huán)境溫度波動較大時，機床溫度場分布和熱變形規(guī)律會發(fā)生較大改變，從而影響熱誤差模型精度。 如何使模型適應(yīng)環(huán)境溫度的變化還需進一步研究。