劉子建　皮皓杰　熊敏

摘 要：為了降低加工過程的熱誤差，提高數控機床加工精度，基于時序相關分析理論與數值計算方法，建立了一種以溫度場分布及加工參數為輸入的新型機床主軸熱誤差建模方法.所建模型由熱誤差模型、主軸動壓軸承熱特性模型以及主軸熱傳遞模型三部分組成.該方法首先根據時序相關理論建立熱誤差與溫度測點之間的相關模型，再通過灰色相關理論完成關鍵溫度測點位置與數量的優(yōu)化，同時，基于數值計算與熱傳導理論，建立了動壓主軸系統(tǒng)熱特性模型.以一臺大型龍門導軌磨床為實驗對象，建立了磨床主軸箱熱誤差預測模型.實驗結果表明，所建立的熱誤差模型具有良好的熱誤差辨識性能.

關鍵詞：時序相關分析；雷諾方程；有限差分法； 熱誤差

中圖分類號：TH161 文獻標識碼：A

機床熱誤差是精密加工誤差的主要來源之一，對其控制和補償是提高機床加工精度的關鍵技術＼[1＼].機床熱誤差補償涉及熱誤差檢測、熱誤差預測模型和熱誤差實時補償三類問題，其中熱誤差預測模型的建立最為關鍵＼[2＼].

機床熱誤差控制的最新技術是基于熱誤差預測模型的誤差實時補償技術（RTEC）＼[3＼].由于機床熱誤差可以看作是關于機床熱分布、位置以及運動方向的確定性函數＼[4＼]，故而，可以建立以機床溫度場分布和運動方向為自變量的熱誤差模型.同時，機床溫度場分布特性可以用若干溫度測點處的數據來表征＼[5＼].

本文運用時間序列分析理論擬合熱誤差與相關溫度測點數據之間的關系，采用數值計算方法建立某高精度磨床主軸熱特性模型，提高了此類熱誤差模型的預測精度和魯棒性.

1 機床熱誤差模型的理論基礎

1.1 時序分析建模方法

時序分析是對觀測所得的有序隨機數據進行分析和處理的一種方法，適用于輸入不可測或不確定測量數據的描述，常被用于機床熱誤差建模.

根據時序分析理論＼[6＼]，對一組測點進行第i次測量的熱誤差χi可以分解為兩部分，如式（1）和式（2）所示，一部分是確定性部分，完全取決于各溫度測點處溫度序列Tr，i的線性組合fi；另一部分是隨機性部分，完全獨立于Tr，i，由殘量εi確定（如儀器的測量誤差）.同時，假定εi是零均值的獨立序列，即白噪聲序列，因而，對于不同i值，εi是相互獨立的.

χi=fi+εi.（1）

fi=β0+β1T1，i+β2T2，i+…+βRTR，i.（2）

式中：i=1，2，…，I，βo為常數，βi為系數，r為溫度測點序號，R為溫度測點總數.

基于熱誤差χi和各溫度測點處溫度Tr，i的實測樣本序列，采用最小二乘法獲取系數β0和βi的估計值0和i，如式（3）所示：

=TTT-1TTχ.（3）

式中，χ=χ1，χ2，…，χIT，

=o，1，2，…，RT，

T=1T1，1T2，1…TR，11T1，2T2，2…TR，21T1，IT2，I…TR，I.

將系數估計值代入式（2）和（1），得到殘量序列εi，如式（4）所示：

εi=χi-o+1T1，i+2T2，i，…，RTR，i.（4）

采用一階自回歸模型描述殘差，如式（5）所示：

εi=φ1εi-1+ai.（5）

采用最小二乘估計求解上述模型的參數，如式（6）所示.同時，ai是均值為0的正態(tài)分布序列：

1=∑Ii=2εiεi-1/∑Ii=2ε2i.（6）

綜上所述，由一組測點確定的熱誤差計算模型如式（7）所示：

χi=0+1T1，i+2T2，i+…+RTR，i+10+1T1，i-1+2T2，i-1+…+RTR，i-1+ai.（7）

1.2 機床溫度測點優(yōu)化方法

機床各部分溫度的分布特性是熱誤差產生的內在原因，因而，溫度測點數量和位置的選取將直接影響式（7）計算誤差的準確性.基于已有的測試數據，采用灰色關聯度分析方法進行溫度測點位置和數量的優(yōu)化選取，可以將測點數目控制在合理范圍內的同時，提高式（7）的計算精度.

灰色關聯度分析是按照系統(tǒng)中各特征參量序列之間的相似程度進行系統(tǒng)分析的一種方法[7-8]，保證各特征參量間的等效性和同序性，是進行關聯度分析的前提.

采用極差變換法的一種改進形式，對熱誤差與各溫度測點的實測結果進行無量綱化處理，如式（8）所示：

x（k）=max x（k）-x（k）max x（k）-min x（k），

k=1，2，…，n.（8）

設對應于熱誤差實測值的x0為參考序列，對應于各溫度測點的溫度實測值的xi為比較序列，則x0對xi在第k點的灰關聯系數如式（9）所示：

ξ0i（k）=min imin kΔ0i（k）+ρmax imax kΔ0i（k）Δ0i（k）+ρmax imax kΔ0i（k）.（9）

其中，Δ0i（k）為第k點x0與xi差值的絕對值，ρ為分辨系數，min imin kΔ0i（k）為Δ0i（k）的兩級最小值，max imax kΔ0i（k）為Δ0i（k）的兩級最大值.兩序列之間的灰關聯度為：

γ0i=1m∑mk=1ξ0i（k）.（10）

式中：γ0i為比較序列xi與參考序列x0的關聯度，m為比較序列的數據個數.

將各個比較序列對同一參考序列的關聯度按大小順序排列起來，即組成了反映各個比較序列對同一參考序列的“主次”關系的關聯序.依據關聯序即可對磨床主軸各溫度測點位置分布及測點數量進行優(yōu)化處理.

1.3 動壓軸承的熱特性建模方法

精密磨床運行時，主軸軸承的摩擦熱是熱誤差的主要來源，通常也是主要的溫度測點.在機床的實際運行中，高速旋轉和結構復雜的特點，難以用傳感器對軸承處的溫度進行準確測量.通過建立主軸的熱特性描述模型，計算不同工況下主軸軸承處的溫度參數用于熱誤差建模，對于提高模型的精度和適用性都有主要意義.

建立主軸動壓滑動軸承油膜熱特性計算模型的思路是通過對雷諾方程、溫粘方程等的聯立求解，得到動壓軸承達到平衡狀態(tài)時的平均溫度.計算流程如圖1所示.

由于剛性徑向滑動軸承的油膜壓強分布滿足雷諾方程，基于相似理論，將坐標原點置于軸承軸線上，且位于軸承寬度的1/2處，對雷諾方程進行無量綱化后得到＼[9＼]：

φ3φ+DB2γh-3γ=6φ.（11）

其中， =ψ2p/6μVtB，ψ為相對間隙，μ為潤滑油動力粘度，p為油膜壓強，Vt為軸頸上任意點的切向速度，φ=x/R，x為軸頸上任意點的周向坐標，R為軸瓦半徑，為油膜無量綱厚度，=1+εcos φ，ε為偏心率，γ=2y/B，y為軸向坐標，B為軸承寬度，D為軸瓦直徑.

式（11）基于如下假設：

1）潤滑油重量及慣性力忽略不計；

2）油膜為牛頓流體，且為層流；

3）潤滑油不可壓縮；

4）潤滑油為定常流體；

5）載荷方向為垂直于軸向.

對于任意動壓軸承油膜，若求解區(qū)域為φ1≤φ≤φ2，0≤γ≤1，則將求解區(qū)域網格化，且令i表示沿φ方向的列數，j表示沿γ方向的行數.根據有限差分原理，將式（11）化為有限差分方程求解.

1.3.1 雷諾方程的有限差分化

采用半步長插入點上壓強的中差商來表示任意節(jié)點上壓強的一階導數，相鄰半步長插入點上壓強的一階導數的中差商來表示壓強的二階導數，則式（11）可表示為：

求解上述方程時，先確定承載區(qū)四周邊界的壓強，進而求解整個求解區(qū)域上的壓強分布.

1.3.2 油膜偏位角的修正

由于已知載荷方向，為了確保油膜壓力分布求解的準確性，需要進行油膜偏位角的校核，從而保證油膜承載力的合力與載荷方向相反.

已知載荷方向為垂直于軸向，則軸承的偏位角應保證油與膜合力方向相反，采用迭代法求解穩(wěn)定狀態(tài)下軸承的偏位角.修正方程為：

1.3.3 溫粘方程的求解

油膜溫升的求解需要考慮溫度對潤滑油粘度的影響.為了克服采用能量方程求解粘度過程的復雜且耗時的缺點，采用潤滑油標準粘度η50（t=50℃）來計算任意潤滑油溫度下的粘度＼[10＼]：

η=ηxexp 400lg η50ηx1t+95-0.001.（17）

其中，t為油膜溫度，ηx=0.18 mPa·s.

將修正后的粘度值代入雷諾方程進行迭代，直到收斂，此時，所得壓強分布即為精確解.

1.3.4 油膜溫升的計算

已知軸頸表面的摩擦阻力為：

Ff=vωR2Bψ∑ζa=1∫10∫φa2φa11+311φdφd.（18）

軸頸表面的摩擦系數為f=f/，其中， =x2+z2.

基于流量連續(xù)原理和差分原理，軸承泄油量為：

Q=∑ζa=12R3ωψBΔφΔγ∑m-103i，n-1-i，n.（19）

基于熱平衡公式，假設由摩擦功轉化的熱量全部由潤滑油及熱擴散帶走，則熱平衡方程式為＼[11＼]：

Δt=f·FfψDB/cρQψUDB+πksψU.（20）

則油膜溫度為：

t=t0+f·FfψDB/cρQψUDB+πksψU.（21）

其中，t0為潤滑油初始溫度； c為潤滑油比熱容；ρ為潤滑油密度；Q為軸承泄油量，ks為散熱系數.

1.4 主軸熱傳遞模型

已知主軸運轉過程中，動壓軸承油膜溫度如式（21）所示，下面討論求解主軸熱特性分布的另一個關鍵條件，即主軸熱傳遞模型.

采用圓柱坐標系描述主軸溫度分布，則不含內熱源的非穩(wěn)態(tài)導熱微分方程為＼[12＼]：

2Tr2+1rTr+2Tz2=1αTt.（22）

其中，T為溫度，r為半徑，z為軸向坐標，t為時間，α為主軸散熱系數.

式（22）還基于如下假設：

1）主軸材質均勻且各向同性；

2）主軸材料的特性系數為常數；

3）物體內無溫度躍變或物體相變；

4）主軸溫度不隨角度φ的變化而變化.

采用分離變量法，假定式（22）有如下變量分離形式的解：

2 熱誤差模型精度分析

下面以某高精度龍門導軌磨床主軸箱系統(tǒng)（如圖2所示）為研究對象，依據上述研究結果進行測試和建模，并對熱誤差模型精度進行分析驗證.

首先根據經驗和初步分析確定12個測試點：主軸后端皮帶輪附近布置2個測點（1～2號測點）；主軸箱體上布置4個測點（3～6號測點）；主軸前段布置3個測點（7～9號測點）；電機上布置1個測點（10號測點）；安裝座上布置2個測點（11～12號測點），如圖2所示.進行溫度測試采樣.測試條件為，磨頭主軸以1 040 r/min和1 500 r/min的恒定轉速各空載運行2.5 h，每間隔5 min用非接觸式紅外溫度測量儀測量12個測點的溫度各一次，與此同時，用千分表測量主軸在X，Y，Z三個方向上對應的熱漂移，總共獲取2組溫度和熱誤差數據.

2.1 主軸箱系統(tǒng)溫度測點優(yōu)化

采用灰色關聯度分析法，將所測得的主軸在X，Y，Z三個方向上的熱誤差作為參考序列，所測得的12組溫度測點數據作為比較序列.利用式（8）對各數據序列進行處理，通過式（9）和式（10）求得各比較序列與參考序列間的灰關聯度，如圖3所示.

經過關聯度大小的對比，排除關聯度小的溫度測點后，確定影響主軸三個方向熱漂移的關鍵溫度測點總數由12個降低為8個，分別是：1號、3號、4號、5號、6號、7號、11號及12號等8個溫度測點.其中，1號與7號測點位于主軸軸承處，需采用2.3節(jié)與2.4節(jié)建立的熱特性模型對該兩個測點處的溫度進行預測.

2.2 主軸熱特性建模與精度分析

實驗中所采用的磨床的主軸箱動壓軸承及其潤滑油的參數包括：寬徑比0.8，初始偏位角90°，偏心率0.5，軸頸半徑0.120 m，軸瓦半徑0.167 m，內切圓半徑0.121 m，軸承張角147°，初始潤滑油動力粘度0.02 Pa·s，初始潤滑油溫度28 ℃，潤滑油比熱1 888 J/（kg·℃），潤滑油密度867 kg/m3，軸承散熱系數80 W/m2·℃，主軸轉速為1 470 r/min或1 500 r/min.

采用主軸熱特性模型和主軸熱傳遞模型，對1號與7號測點處溫度時間序列進行預測，同時，將預測值序列i與兩個測點處的溫度實測值序列xi進行對比，并計算其相對誤差δ，結果如表1所示.相對誤差的計算公式如式（27）所示：

δ=xi-ixi.（27）

由表1可知：預測值序列與實測值序列的相對誤差較小，相對誤差平均值為14.2%.其中，1號測點相對誤差最大值13.3%，7號測點相對誤差36.4%.

因此，主軸熱特性辨識模型的預測數據與實測數據基本吻合，另外考慮到1號和7號測點位于主軸外露的軸頸處，測試所得溫度與軸承處的溫度本身存在一定的誤差，所以，本文建立的模型能有效地描述主軸的熱特性.

2.3 熱誤差模型實驗驗證

綜合主軸熱特性辨識數據和溫度測點數據，利用式（7）進行熱誤差模型的驗證.

首先，基于第1組數據中各溫度測點數據序列和熱誤差數據序列，建立主軸X，Y，Z三個方向上的熱誤差預測模型.然后，選取未使用過的第2組數據中的溫度數據，代入熱誤差預測模型中，輸出X，Y，Z方向的熱誤差預測值；最后，將預測值與熱誤差實測值比較.同時，預測過程中，第1號和7號測點的溫度數據采用主軸箱熱特性模型的輸出值.模型的熱誤差預測值與實測值，以及兩者之間的殘差如圖4，圖5和圖6所示.

3 結 論

1）基于空載實驗的精密導軌磨床熱誤差建模方法，研究的熱誤差不受其它加工因素的耦合影響，降低了熱誤差建模的復雜程度，對于切削力不大的精密磨削具有實用意義.

2）基于時序相關分析理論與數值計算理論的熱誤差模型能有效地對機床熱誤差進行辨識.

3）實驗驗證了該建模方法的有效性與可行性，為機床精密加工精度的提高與熱誤差實時補償技術提供了理論借鑒.

參考文獻

[1] 林偉青， 傅建中， 許亞洲， 等. 基于最小二乘向量機的數控機床熱誤差預測[J]. 浙江大學學報， 2008， 42（6）： 905-908.

LIN Weiqing， FU Jianzhong， XU Yazhou， et al. Thermal error prediction of numerical control machine tools based on least squares support vector machines [J]. Journal of Zhejiang University， 2008， 42（6）： 905-908.（In Chinese）

[2] 李永祥， 楊建國.灰色系統(tǒng)模型在機床熱誤差建模中的應用[J]. 中國機械工程， 2006， 17（23）： 2439-2442.

LI Yongxiang， YANG Jianguo. Application of grey system model to thermal error modeling on machine tools [J]. China Mechanical Engineering， 2006， 17（23）： 2439-2442.（In Chinese）

[3] YANG Hong， NI Jun. Dynamic neural network modeling for nonlinear， nonstationary machine tool thermally induced error[J]. International Journal of Machine Tool& Manufacture， 2005， 45（4）： 455-465.

[4] JOSEF M， JERZY J， ECKART U， et al. Thermal issues in machine tools [J]. Annals of CIRP， 2012， 61（2）： 771-791.

[5] WECKA M， MCKEOWNB P， BONSEC R， et al. Reduction and compensation of thermal errors in machine tools [J]. Annals of CIRP， 1995， 44（2）： 589-598.

[6] 楊叔子， 吳雅， 軒建平， 等.時間序列分析的工程應用[M]. 武漢：華中科技大學出版社， 2007.

YANG Shuzi， WU Ya， XUAN Jianping， et al. Engineering application of time series analysis[M]. Wuhan： Huazhong University of Science and Technology Press， 2007.（In Chinese）

[7] 劉思峰， 黨耀國， 方志耕.灰色系統(tǒng)理論及其應用[M]. 北京：科學出版社， 2005.

LIU Sifeng， DANG Yaoguo， FANG Zhigeng. Theory and application of grey system[M]. Beijing： Science Press， 2005.（In Chinese）

[8] 閆嘉鈺， 張宏韜， 劉國良，等. 基于灰色綜合關聯度的數控機床溫度測點分組優(yōu)化[J]. 湖南大學學報： 自然科學版，2008， 35（4）： 37-41.

YAN Jiayu， ZHANG Hongtao， LIU Guoliang， et al. Optimization of measuring points for machine tool thermal error modeling based on grouping of synthetic grey correlation method [J]. Journal of Hunan University： Natural Science， 2008， 35（4）： 17-41.（In Chinese）

[9] 楊沛然. 流體潤滑數值分析[M]. 北京：國防工業(yè)出版社， 2005.

YANG Peiran. Numerical analysis of fluid lubrication[M]. Beijing： National Defense Industry Press， 2005.（In Chinese）

[10]李中華. 數控機床主軸軸承部件溫度及熱力學建模[J]. 機械制造， 2010， 48（6）： 30-33.

LI Zhonghua. Temperature of spindle bearings of CNC machine tools and thermodynamic modeling[J]. Mechanical Manufacturing， 2010， 48（6）： 30-33.（In Chinese）

[11]孫軍， 桂長林. 內燃機曲軸強度研究的現狀、討論與展望[J]. 內燃機學報， 2002， 20（2）： 179-184.

SUN Jun， GUI Changlin. A review of crankshaft strength analysis for internal combustion engines [J]. Journal of Combustion Engine， 2002， 20（2）： 179-184.（In Chinese）

[12]俞昌銘. 熱傳導[M]. 北京：高等教育出版社， 1983.

YU Changming. Heat conduction [M]. Beijing： Higher Education Press， 1983.（In Chinese）