基于熱特性分析的幾何與熱復合的定位誤差建模及補償

馮文龍，南江紅，楊建國

(1. 上海航天電子技術研究所，上海 201109； 2. 上海交通大學 機械與動力工程學院，上海 200240)

0 引言

為能生產高精度的工業產品，減小加工過程中的各項誤差是一種最為有效的途徑。影響加工精度的誤差包括各部件的幾何誤差、熱誤差、切削力誤差、伺服控制誤差及插補誤差等[1]。其中，熱誤差的比重達到了40%～70%[2]。對于絲杠螺母系統，絲杠受內部加工熱源的影響，其溫度發生較大變化，從而受熱膨脹產生熱誤差。內部加工熱源包括絲杠螺母摩擦熱源、軸承熱源、變速箱熱源、電機熱源等[3]。在熱誤差的影響下，機床加工誤差往往會超出容許誤差10倍以上，故控制機床熱誤差越來越受到人們重視[4]。

減小熱誤差的方法一般可分為誤差補償法和誤差防止法[5-7]。誤差防止法是通過優化機床結構設計和提高各部件加工精度的方法來提高機床的整體性能。但因存在多種物理限制，此種方法對機床精度的提升非常有限，且應用此方法消耗的成本較大[8]。誤差補償法是人為制造一項誤差用于補償，該誤差與原始誤差相比，大小相等，方向相反，可抵消原始誤差[9]。基于現有先進的誤差識別技術的支持，機床原始誤差能被激光干涉儀等測量設備測得。相較于誤差防止法，誤差補償不乏是一種經濟有效的精度提升方案。文獻[10]中提出了基于正交多項式的熱誤差建模方法，補償后，在既定運動軌跡下，幾何誤差及其熱誤差減小了80%以上，但其使用單個溫度點來表征整個絲杠的熱狀態具有先天性缺陷，一旦絲杠的運動范圍改變，該熱誤差模型將不再適用于變化后的加工環境，誤差補償效果也將大打折扣。

本文考慮復雜的加工狀況，即絲杠并非處于熱均勻狀態。把絲桿軸細分成多個細小的熱均勻單元，并將每個細小熱均勻單元應用到熱均勻絲杠定位誤差模型，計算可得任意時刻每個細小單元的幾何與熱復合的定位誤差值，繼而把這些細小單元的定位誤差模型疊加，建立整個絲杠軸的幾何與熱復合的定位誤差模型。通過實驗驗證，結果表明：補償后，絲杠螺母傳動系統的熱誤差明顯減小，機床加工精度顯著提高。

1 絲杠熱特性分析

為研究絲杠的熱誤差，利用立式加工中心(VMC850E)展開實驗。

以立式加工中心的Y軸為例,絲杠螺母傳動系統機構及其傳動機理如圖1所示，在傳動過程中，系統根據需要移動的位移量，依據其絲杠和螺母的螺距值，可精確計算得到伺服電機需要旋轉的角度值。然而絲杠本身存在螺距誤差，且當外界環境變化時，其螺距誤差也發生變化，是一項時變誤差。雖然該誤差導致的絲杠轉動圈數相同，但其帶來的位移變化卻并不相同。故對絲杠螺母系統定位誤差的研究，可轉化為對絲杠螺距誤差的研究，即對絲杠伸長量的研究。絲杠軸正常工作時，會受到電機熱源、前后軸承熱源、環境溫度變化以及絲杠和螺母摩擦熱源的影響，從而引起絲杠本身的熱溫升，繼而發生熱伸長，導致絲杠螺距發生變化。可在電機外殼、前后軸承座、絲杠螺母以及床身等多處布置溫度傳感器，并實時監測這些位置的溫度，來評估其對絲杠伸長的影響。

在對絲杠展開研究之前，為便于理解，作如下假設：

1) 絲杠被認為是一個實心的圓柱體，其熱伸長或者冷縮現象均被認為是線性的。

2) 絲杠表面與內部(截面上)的溫度一致。

3) 以空氣熱傳導方式影響絲杠溫度的熱源，實驗證明其影響較小，如電機熱源、前后軸承熱源等,均不在本文考慮范圍之內。本文認為絲杠和螺母的摩擦熱源是導致絲杠伸長的主要因素。

4) 在同一進給速度下，認為絲杠和螺母摩擦產生的熱量是相同的。

5) 將絲杠等分成多個細小的單元，并假設這些單元是熱均勻的，即在任何時候，獨立單元的每個位置均處于同一溫度下，且任何一個獨立單元，只要其處在同一溫度下，該單元的熱狀態就認為是相同的，與歷史過程無關，即每個單元的溫度值與熱伸長量一一對應。

絲杠受熱源的影響，會發生熱脹冷縮現象，其長度變化量可表示為

ΔLscr=ΔTscr·L·α

(1)

式中：ΔLscr表示絲杠長度的變化量，單位為μm，正值表示伸長，負值表示縮短；ΔTscr表示絲杠溫度的變化量，單位為℃，正值表示升高，負值表示降低；L表示絲杠的總長度，單位為m；α表示絲杠的線性熱膨脹系數，即絲杠改變1 ℃長度的變化量，取值為11.7×10-3，單位為℃-1。

絲杠的熱伸長量等同于絲杠的定位熱誤差量，使用激光干涉儀，對絲杠的定位誤差進行測量，測量結果如圖2所示，每個時間點的定位誤差曲線形狀基本保持不變，每條曲線可由0 min這條曲線繞著起始點旋轉一定角度得到。故所測得的定位誤差數據可拆分成兩部分，一部分是固定的幾何定位誤差，只跟坐標位置信息相關，也即0 min這條曲線；另一部分是熱復合定位誤差，跟旋轉角度及坐標信息相關，將旋轉角度對應的斜率變化定義為此時絲杠的熱狀態參數。實測定位誤差采用誤差分離技術，表示為

δp(y,T)=δpG(y)+δpT(T)·y

(2)

式中：δp(y,T)表示跟Y軸坐標位置及溫度相關的定位誤差；δpG(y)為只跟Y軸位置相關的幾何定位誤差；δpT(T)表示熱狀態參數，即跟旋轉角相關的斜率變化值；y表示對應的Y軸坐標值。

圖2 Y軸幾何與熱復合的定位誤差Fig.2 Geometrical error and thermally induced positioning error of Y-axis

絲杠定位熱誤差等同于絲杠的熱伸長量，對比式(1)和式(2)，可得到熱狀態參數的表達式

δpT(T)=ΔTscr·α

(3)

由式(3)可知，熱狀態參數δpT(T)跟絲杠的溫升量呈線性關系，只要能準確地測量出絲杠的實時溫度，即可求解熱狀態參數，繼而獲得定位熱誤差，完成對其的建模過程。然而絲杠在工作條件下將一直處于旋轉狀態，無法在其表面直接吸附傳感器用以測量其溫度值，可通過測量絲杠螺母的實時溫度，來近似得到絲杠的溫度。在絲杠螺母摩擦的過程中，絲杠會受絲杠螺母摩擦熱源的影響，溫度升高，同時，絲杠也會向低于其溫度的外界空氣釋放出熱量，溫度降低。絲杠在溫升階段的熱狀態可由以下瞬態熱平衡方程表示

(4)

式中：Qnut為絲杠和螺母摩擦產生的熱量；hw為強制對流熱傳導系數；A為熱擴散面積；Tscr、Tamb分別為絲杠溫度和環境溫度；ρ為絲杠材料密度；c為絲杠材料比熱容；V為絲杠體積；tup為絲杠溫升的累積時間。

為得到絲杠瞬態溫度Tscr的表達式，對式(4)作解微分方程運算。絲杠瞬態溫度可表示為

(5)

而當絲杠停止轉動，絲杠和螺母停止摩擦，機床處于停止狀態時，絲杠不再受摩擦熱源的影響，而只是向低于其溫度的外界空氣散熱，此時絲杠處于瞬態熱平衡，熱平衡方程為

(6)

式中：hc表示自由對流熱傳導系數；tdown表示絲杠降溫累積時間。

對式(6)作解微分方程運算，得到降溫狀態時絲杠瞬態溫度表達式，即

(7)

絲杠被認為是一個實心的圓柱體，故其體積V和散熱面積A都可用絲杠直徑d來表示，綜合式(2)、(3)、(5)和式(7)，可得到溫升和冷卻階段絲杠定位誤差求解表達式：

(8)

由式(8)可知，絲杠的定位誤差與絲杠和螺母的摩擦熱量、絲杠直徑、散熱面積、強制對流系數、自由對流熱傳導系數、溫升時間和冷卻時間等參數相關。其中散熱面積和運動范圍相關，絲杠和螺母的摩擦熱源與進給速度相關。當加工參數發生變化時，絲杠溫度預測方程式也可通過調整模型系數，來適應環境的變化，以此保證絲杠溫度函數對絲杠溫度的準確預測。摩擦熱求解表達式為

(9)

式中：Qnut為絲杠螺母摩擦產生的發熱量，單位為W；M為絲杠螺母的摩擦力矩，單位為N·mm；n為絲杠的轉速，單位為r/min；a為系數；F為機床進給速度，單位為mm/min；S為絲杠的導程，單位為mm。

2 絲杠定位誤差建模

2.1 絲杠定位誤差模型理論

分析絲杠在工作環境中的熱狀態可得絲杠在溫升和冷卻階段的定位誤差模型，但其表達式中仍有很多未知的參數，不能直接用于定位誤差的計算，本節將求解這些未知系數。

以2 000 mm/min的進給速度來回移動Y軸工作臺，移動范圍為整個行程。絲杠溫度會伴隨著移動而上升，在Y軸螺母處布置溫度傳感器，用以測量和螺母接觸處絲杠的溫度。絲杠和螺母摩擦共同產生摩擦熱，且因絲杠和螺母緊密連接，故可認為在絲杠和螺母接觸處絲杠和螺母的溫度是相同的。與此同時，在靠近絲杠螺母傳動系統后軸承的床身處布置溫度傳感器，監控絲杠溫升量基準，因整個實驗過程耗時較長，很難保證整個過程中環境溫度不發生變化。而絲杠螺母傳動系統中，后軸承處受到絲杠螺母摩擦熱源影響較小，其溫度僅受環境溫度變化的影響，能有效測量出絲杠在環境溫度影響中的溫度變化，螺母處溫度傳感器值和此處溫度的差值即絲杠僅受絲杠螺母摩擦熱源影響時的溫度變化量。移動Y軸工作臺直到絲杠到達熱平衡狀態，繼而停止工作臺，進入絲杠冷卻階段，在這個過程中，實時測量并記錄螺母溫度和絲杠基準溫度。螺母位置和軸承端位置的溫度傳感器布置情況如圖3所示。

圖3 溫度傳感器安裝圖Fig.3 Temperature sensors installation on machine tool

圖4 溫升和冷卻階段絲杠溫度變化量Fig.4 Temperature variations during warming and cooling phases

指數擬合測得的溫度曲線，可得絲杠溫度的預測模型。測量得到的絲杠溫度數據擬合情況如圖4所示。擬合得到的絲杠溫度的預測模型可表示為

(10)

由式(10)可知，冷卻階段的指數項系數與溫升階段不同，究其原因：理論模型的前提是溫升到熱平衡后再降溫，其降溫的起點是熱平衡狀態，但在實際實驗中，溫升階段歷時140 min，看起來曲線變化趨勢已經逼近熱平衡狀態，但其終究還是沒有達到熱平衡，若在此點降溫，其相距熱平衡狀態必有一個時間間隔。對比式(10)和式(5)、(7)，可得

(11)

式中：絲杠密度ρ=7.85×103kg/m3；絲杠比熱容c=0.46×103J/(kg·℃)；絲杠直徑d=40 mm；絲杠長度L=850 mm；絲杠散熱面積A=0.11 m2。

根據上述已知參數，求解式(11)，可以得到：hw=13.4 W/(m2·℃),Qnut=16.5 W,hc=12.0 W/(m2·℃)。

把求得的參數代入式(8)，其中線性熱膨脹系數取值α=11.7×10-3μm/mm，可得絲杠定位誤差預測模型，其表達式為

δp(y,T)=

(12)

在實際加工過程中，在溫升階段，依據定位誤差模型，實時定位誤差可易得；在熱平衡后冷卻階段，實時定位誤差也可輕松獲得。而在實際的加工中，很少出現理想的溫升到熱平衡再冷卻的情況，一般都是溫升一段時間就開始冷卻，冷卻一段時間又開始溫升，溫升和冷卻階段隨機交替出現。故下文將研究上述溫升和冷卻階段隨機交替出現時定位誤差的計算方法。

如圖5所示，當絲杠從冷態開始移動，一直處于溫升狀態，熱狀態參數也隨著絲杠溫度的上升而變大，但在A點，絲杠停止運動，開始冷卻，此時的絲杠可看作一個均勻的整體，且絲杠各點溫度基本相同，找到與其熱狀態相同的一個點，且處于自由冷卻中，復制其變化趨勢就能得到絲杠中間過程冷卻的熱狀態參數變化趨勢。作一條平行于時間軸的直線，交冷卻階段曲線于B點，平移曲線BD至A、B點重合即可得到曲線AC，絲杠中間過程冷卻后的熱狀態參數變化趨勢如曲線AC所示。

圖5 中間過程冷卻情況下的熱狀態參數變化趨勢Fig.5 Thermal parameter variation when cooling in the midway of warming up

絲杠在溫升過程中，中途進行冷卻后，絲杠處于冷卻過程中，此時，絲杠重新開始移動，絲杠溫度必定又會上升，絲杠的熱狀態參數變化趨勢如圖6所示，絲杠在經歷第1次溫升段OA后，在A點冷卻，冷卻至F點后，再度溫升，如曲線段FG，直至熱平衡。過F點作時間軸的平行線，交第1次溫升曲線OH于E點，平移曲線EH直至E點和F點重合，即可獲得2次溫升曲線段FG。在絲杠溫升階段FG過程中，如果再次出現冷卻和冷卻后再溫升的情況，可通過類似的方法計算絲杠的熱狀態參數，即當絲杠冷卻階段和溫升階段隨機交替時絲杠熱狀態參數都可通過上述計算方法計算得到。

圖6 絲杠二次溫升對應的熱狀態參數變化趨勢Fig.6 Thermal parameter variation in double warming situation

2.2 分段指數模型理論

2.1節建立了幾何與熱復合的定位誤差模型，兼容溫升及冷卻交替出現情況。該模型利用絲杠全行程來回移動模擬加工運動，把絲杠等效為一個均勻的整體，繼而研究絲杠的熱狀態。但在實際的加工中，基本不存在全行程來回移動的情況，實際的加工狀況是圍繞工件，在全行程的某個區間內來回移動。實際的加工條件和2.1節中建立的誤差模型前提不相符，這也預示著建立的誤差模型不能直接應用于實際加工的誤差補償中。

在實際的加工中，加工工件不同，運行的G代碼程序不同。機床的運動軌跡是無規律及不可預測的，這勢必帶來一個嚴重問題，即絲杠各點受到摩擦熱源的影響將不再相同，絲杠各點的熱狀態將不再相等。使用一個熱狀態參數已經不能表征整根絲杠的熱狀態。本文采用對絲杠進行細分的方法，把絲杠等分成若干足夠小的單元，因其足夠小，其內部就不存在溫度梯度，近似可認為其絕對熱均勻，對每個等分單元應用2.1節建立的誤差模型，得到每個等分單元的幾何與熱復合的定位誤差模型，再將這些等分單元的誤差模型疊加，即可得到整根絲杠的誤差模型。

為驗證這一方法的可行性，實驗研究機床的Y向運動。為模擬實際工件的加工狀況，使絲杠各點處于不同的熱狀態。本實驗采用各區域依次熱機的方法，熱機軌跡為：

1) 螺母工作臺在位置坐標為200 mm和600 mm之間以進給速度2 000 mm/min來回移動10 min；

2) 螺母工作臺再在位置坐標為0 mm和400 mm之間以進給速度2 000 mm/min來回移動10 min；

3) 螺母工作臺在位置坐標為400 mm和800 mm之間以進給速度2 000 mm/min來回移動10 min。

運動過程中，實時記錄機床的位置坐標值，如圖7所示。絲杠在這30 min溫升階段的熱狀態參數模型將不同于上述全行程運動溫升過程，因其運動范圍只是全行程的一半，故其散熱面積是全行程溫升情況的一半，用A2=A/2表示。上述30 min溫升過程中，絲杠溫升和冷卻時熱狀態參數計算方法可表示為

圖7 溫升階段工作臺位置坐標Fig.7 Recorded location of nut during warming phase

經過實驗驗證，選定最小的熱均勻單元長度為40 mm，將Y軸等分為20份，對這20個細小熱均勻單元分別建模。建模結果如表1所示。

表1 20個單元的熱狀態參數值

20個獨立單元的熱狀態參數分別用δpT1,δpT2,δpT3,…,δpT20表示，整根絲杠的定位誤差可表示為

(14)

式中：δpG(y)表示Y軸的幾何定位誤差；p為絲杠的位置，單位為mm。

幾何定位誤差可通過最小二乘法建模，建模結果可表示為

幾何定位誤差和熱狀態參數建模結束后，依據式(14)，即可預測絲杠任意時刻的定位誤差。為評估定位誤差模型的模型精度，在上述30 min熱機過程中，每隔10 min測量一次絲杠的定位誤差，并對比測量誤差值和模型預測誤差值之間的差異，用以驗證模型的準確性。對比結果如圖8所示，結果表明：上文所建立的定位誤差模型能準確預測上述30 min實驗過程中的定位誤差，該誤差在5 μm以內，預測精度大于80%。

此定位誤差建模方法對運動軌跡沒有諸多限制，其適用于任何加工軌跡后的定位誤差預測，可廣泛應用于實際生產加工過程中，只需要同時對每個獨立細分單元做熱狀態分析，準確記錄每個單元的溫升和冷卻時間，并依據模型算法計算每個單元各個時刻的熱狀態，根據式(14)疊加，即可獲得絲杠定位誤差模型。用戶可根據自己的要求來細分絲杠，細分單元越小，模型精度越高，定位誤差建模流程圖如圖9所示。

圖8 定位誤差模型預測值與實測值對比圖Fig.8 Comparison of positioning errors between the measured values and the predicted values

圖9 定位誤差建模流程圖Fig.9 Algorithm flow chart of positioning error mode

3 補償結果分析

為驗證模型的魯棒性和有效性，在立式加工中心VMC850E上進行一系列驗證實驗。首先，機床空載運行，以2 000 mm/min的進給速度在位置為300 mm和500 mm之間來回移動30 min，分別測量機床Y軸補償前和補償后的定位誤差；繼而，機床停機10 min，自由冷卻，然后再分別測量機床Y軸補償前和補償后的定位誤差，對比這兩組補償前后的定位誤差量，評估補償的效果。對比結果如圖10所示，在300～500 mm段定位誤差在溫升后有明顯增加，冷卻后又有明顯下降。經過指數模型補償后，定位誤差從30 μm減小至8 μm，然而多項式模型補償后殘差更大，這是因為用單個溫度傳感器來表征整個絲杠的熱狀態對復雜的運動情況已不再適用。

圖10 補償結果對比Fig.10 Comparison of compensation results

接下來進行加工實驗，對6個試件進行定位孔加工，前2個試件分別在有補償和無補償的條件下進行，因從冷態開始加工，這2個試件的加工條件認為是初始狀態；再加工2個試件，這段過程被認為是溫升過程；最后加工2個試件，分別在有補償和無補償的條件下進行，這2個試件的加工條件認為是溫升后狀態。加工過程和試件如圖11所示。

對加工完成試件的4個孔距送交第三方檢測機構“上海市計量測試技術研究院”檢測，測試結果如表2所示，孔距從最大10 μm減小為3 μm，補償效果顯著。

圖11 試件加工過程Fig.11 Machining process of test pieces

mm

4 結束語

針對旋轉編碼器定位的絲杠螺母傳動系統，在溫升和冷卻階段分析絲杠軸的熱特性，建立絲杠軸瞬態溫度和時間的指數關系表達式，根據溫升或冷卻時間計算絲杠實時瞬態溫度。基于熱膨脹理論，根據絲杠軸溫度，建立熱均勻絲杠幾何與熱復合的定位誤差模型。考慮到復雜的加工狀況，絲杠并非處于熱均勻狀態，把絲桿軸細分成多個細小熱均勻單元，對每個細小熱均勻單元應用熱均勻絲杠定位誤差模型，可計算得到任意時刻每個細小單元的幾何與熱復合的定位誤差值，繼而把這些細小單元的定位誤差模型疊加，即可建立整個絲杠軸的幾何與熱復合的定位誤差模型。應用該誤差模型對多個試件進行加工，經過補償后，試件的孔距精度明顯提高。研究結果表明：應用該誤差補償方法，可有效地減小加工熱源影響時批量生產零件的尺寸偏差，相比于基于單溫度點的誤差補償方法，該補償方法更有效。不僅如此，因其是對絲杠螺母傳動系統進行建模，故其對機床結構約束較小，可適用于大部分現有機床結構。針對大批量的零件加工任務，可實現各加工機床之間組網，并上載顯示各機床的加工狀態，實現對各機床加工狀態的實時監控，此項工作可作為本文的后續研究方向。

[1] LEE D S, CHOI J Y, CHOI D H. ICA based thermal source extraction and thermal distortion compensation method for a machine tool[J]. International Journal of Machine Tools and Manufacture, 2003, 43(6): 589-597.

[2] RAMESH R, MANNAN M A, POO A N. Error compensation in machine tools: a review Part II: thermal errors[J]. International Journal of Machine Tools and Manufacture, 2000, 40(9): 1257-1284.

[3] XUZ Z, LIU X J, KIM H K, et al. Thermal error forecast and performance evaluation for an air-cooling ball screw system[J]. International Journal of Machine Tools and Manufacture, 2011, 51 (7/8): 605-611.

[4] RAMESH R, MANNAN M A, POO A N. Support vector machines model for classification of thermal error in machine tools[J]. The International Journal of Advanced Manufacturing Technology , 2002, 20(2): 114-120.

[5] NI J. CNC machine accuracy enhancement through real-time error compensation[J]. Journal of Manufacturing Science and Engineering, 1997, 119(4B): 717-725.

[6] 楊建國. 數控機床誤差綜合補償技術及應用[D]. 上海: 上海交通大學, 1998.

[7] 馮文龍. 大型數控機床多誤差元素建模及綜合補償[D]. 上海: 上海交通大學, 2016.

[8] WANG W, ZHANG Y, YANG J G, et al. Geometric and thermal error compensation for CNC milling machines based on Newton interpolation method[J]. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2013, 227(4): 771-778.

[9] HAN Z Y, JIN H Y, LIU Y L, et al. A Review of Geometric Error Modeling and Error Detection for CNC Machine Tool[J]. Applied Mechanics and Materials, 2013, 303/304/305/306: 627-631.

[10] FAN K G, YANG J G, YANG L Y. Orthogonal polynomials-based thermally induced spindle and geometric error modeling and compensation[J]. The International Journal of Advanced Manufacturing Technology, 2013, 65(9/10/11/12): 1791-1800.