高速主軸系統熱特性分析與實驗

馬 馳，楊 軍，趙 亮，梅雪松，施 虎，王新孟

（西安交通大學 機械制造系統工程國家重點實驗室，陜西 西安710049）

高速與高精度加工是現代制造業的發展方向，而高速主軸的應用對實現高速與高精度加工至關重要.但是，為了改善傳動性能，高速主軸引入內置電機與高速軸承，導致主軸受內部與外部多種熱源的影響，加之主軸部件材料性能各異，且結構復雜，導致熱源位置、強度、膨脹系數及主軸結構的相互耦合產生了復雜的熱特性，導致主軸溫度場分布不均勻引起了熱誤差，并具有時變非線性非穩態特性；大量研究表明，熱誤差是導致機床精度降低的主要因素之一，其中熱誤差約占機床總誤差的70%［1-2］，且機床越精密，熱誤差占機床總誤差的比重越大，可見，熱誤差是引起數控機床精度降低的主要誤差源.因此，研究高速主軸系統熱特性能夠有效避免實際加工中由于溫度場分布不均導致機床加工精度降低的問題.

提高機床加工精度是學者們的畢生追求，熱誤差建模與補償對提高機床加工精度至關重要，通過模糊聚類分析，楊建國等［3］應用回歸分析方法建立了主軸熱變形與溫度場間的線性函數映射關系，YANG 等［4］利用神經網絡建立主軸熱誤差與溫度間的非線性函數映射關系.楊軍等［5］利用時間序列模型的預測精度高、計算速度快等特點建立了熱誤差建模.進給系統熱特性研究也是提高機床加工精度的一個研究方向，楊軍等［6］利用激光干涉儀，紅外熱像儀及溫度位移測量系統，研究了不同進給速度下的機床熱特性.近年來，利用有限差分（FDM）與有限元（FEM）理論分析機床溫度場分布規律及熱變形變化的方法成為研究熱點，Bossmann等［7］使用有限差分法研究了主軸系統生熱機理與散熱機制，Zhao等［8］建立了主軸熱特性分析模型，并研究了主軸溫度場分布及熱變形變化規律，Creighton等

［9］及Kyung等［10］研究了主軸-軸承系統在高速運轉時的動態與熱態特性，Holkup等［11］與Kolar等

［12］建立了預測主軸溫度場分布、熱伸長、軸承剛度、動態特性與接觸載荷的熱力學模型.上述模型都將熱載荷與對流換熱系數等作為重要的邊界條件，但沒有意識到結合面與連續體傳熱間的區別，因此，忽略了接觸熱導對溫度場與熱變形等仿真結果的影響，有些學者將各種不同結合面間接觸熱導設置為經驗常數，導致熱特性仿真精度較低.事實上，結合面間接觸熱導與接觸材料的表面粗糙形貌、接觸壓力、微凸體變形形式、材料屬性及溫度等因素是非線性函數關系，因此，在進行熱特性仿真時，忽略結合面間接觸熱導或將設置為經驗常數是不合理的.另一方面，主軸系統中有大量的結合面，如軸承與軸、軸承與軸承座、主軸箱與端蓋等，當兩粗糙表面接觸時，由于表面粗糙度的存在產生了非完美接觸，即使結合面間壓力很大，兩粗糙表面的非完美接觸也會使得熱量傳遞在結合面處受阻，而產生明顯的溫度降，尤其是靠近熱源處的接觸熱導對主軸熱特性仿真結果影響很大.

本文以高速主軸為研究對象，建立了三維有限元模型進行高速主軸的瞬態熱-結構耦合分析.建立基于幾何-力學-熱分析的綜合接觸熱導預測模型以考慮結合面間接觸熱導對主軸溫度場與熱變形仿真結果的影響，基于分形理論建立了通用的接觸熱導預測模型，考慮了粗糙表面形貌與微凸體3種接觸變形對接觸熱導的影響，有效避免了基于傳統統計學理論與實驗測量法的精度不足及通用性弱等缺點；運用擬靜力學分析方法求解角接觸球軸承生熱功率；應用數值傳熱學基本理論求解主軸部件的對流換熱系數，將上述邊界條件施加到主軸熱-結構分析模型中以提高仿真結果的精度，并進行熱特性實驗驗證該有限元模型的正確性，為高速主軸熱平衡設計提供有益幫助.

1 主軸系統

本文的建模對象為高速電主軸，實體模型如圖1所示.主軸最高轉速nmax＝30 000r／min，主軸兩端安裝有一對型號為FAG B7205-C-T-P4S的角接觸球軸承，其安裝方式為背對背，使用預緊螺母對軸承進行壓緊，預緊力為600N，主軸冷卻系統的工作機理為通過感應出水口溫度控制系統溫度變化，當出水口溫度高于零界溫度時，冷卻系統開始工作，2個配重盤安裝在主軸兩端以調節主軸回轉誤差，4個溫度傳感器安裝在主軸系統內部以測量前軸承、后軸承、電機定子與冷卻水套的溫度.

圖1 高速主軸實體模型Fig.1 Solid model of high-speed spindle

2 主軸熱特性分析模型

為了得到溫度場與熱變形的準確數值解，應對實體模型進行合理簡化，且對實體模型進行恰當的網格劃分，仿真結果的可信度依賴于邊界條件如接觸熱導、熱載荷及對流換熱系數等設置是否合理.對該有限元模型做如下假設：1）忽略主軸系統與環境間的相互熱輻射；2）假設冷卻水套內流體流速足夠大，且恒定，流體的速度波動可以忽略；3）假設環境溫度恒定為15 ℃.

2.1 有限元模型

對螺栓孔、鍵槽與內外倒角等對主軸溫度場與熱變形仿真結果影響很小的結構進行合理簡化，考慮到軸承轉速很高，可將軸承簡化為一空心圓柱體.使用具有20節點的正六面體網格對實體模型進行剖分，為了保證仿真結果的可信度，并保證在整個熱平衡過程中主軸溫度場分布和熱變形的變化與主軸轉速嚴格同步，對溫度梯度較大的區域網格劃分較細，如軸承與電機處，該有限元模型共有298 147個實體單元與1165 521個節點，如圖2所示.

在沒有熱量傳遞或者熱量傳遞不明顯的位置，接觸熱導對結合面間熱量傳遞的阻礙作用并不明顯，而靠近熱源處的接觸熱導對結合面間傳熱影響很大，如軸承內圈與軸，軸承外圈與軸承座處的接觸熱導，這些位置處的接觸熱導直接影響由軸承傳遞進入軸與軸承座熱量的多少，而由軸承傳遞進入軸與軸承座的熱量對主軸的溫度場分布與熱變形變化影響很大，在高速主軸系統中，軸承是主要熱源，由于軸承內圈與軸，軸承外圈與軸承座處的接觸熱導靠近熱源，其余位置如主軸箱與端蓋處的接觸熱導遠離熱源，因此，除軸承處的接觸熱導外，其余位置處的接觸熱導對主軸系統的熱特性的影響可以忽略.主軸系統部件的材料如表1 所示.表中，ρ 與c分別表示主軸部件材料的密度與比熱容，E 為彈性模量，ν 為泊松比，α 為熱伸長系數，λ 為熱導率.瞬態熱分析的控制方程可表示為

式中：T 為熱力學溫度，T ＝T（x，y，z，t）為隨時間與位置分布的溫度場，λx、λy與λz分別表示x，y 與z方向的熱導率，Q 為內部熱源.

圖2 有限元模型Fig.2 Finite element model

表1 主軸系統材料Tab.1 Materials of spindle system

基于胡克定律計算熱變形，如方程（2）所示

式中：ε為應變，ΔT 為溫升.

2.2 接觸熱導模型

由于沒有意識到連續體與結合面間傳熱的區別，因此，傳統主軸熱特性分析模型［7-12］都忽略了接觸熱導對溫度場與熱變形仿真結果的影響，一些學者將接觸熱導設置為經驗常數，導致了溫度場與熱變形的仿真精度較低；另一方面，從微觀角度看任何表面都是粗糙的，即粗糙表面是由大量不同尺度的微凸體組成，因此，兩粗糙表面間的接觸發生在有限離散點上，當熱流通過非完美接觸表面時，熱流會在這些離散點處發生收縮，即當熱量從一個物體傳遞到另一個物體時，在結合面處受到阻礙，接觸熱導的存在導致了熱量不能順利通過結合面，這樣在結合面兩側會產生明顯的溫度降.現有文獻中關于接觸熱導的研究方法大多是基于統計學理論［13-14］與實驗測量法［14-15］.然而，基于統計學理論的研究方法［13-14］對采樣長度與測量儀器的分辨率很敏感，而實驗直接法［14-15］也有其缺點，如通用性不強，適應性差，并且實驗結果具有一定的不確定度，且當材料、表面粗糙度、加工方法與表面接觸壓力等任一因素改變時，就要重新設計實驗；并且，結合面間接觸熱導的預測是一個跨學科問題，包括，即幾何、力學與熱3個子問題，且接觸熱導分析的核心是幾何及力學問題，因此，在進行主軸熱特性分析中，提出基于幾何-力學-熱分析綜合的接觸熱導模型以考慮其對主軸熱特性仿真結果的影響.

2.2.1 分形表面表征 分形幾何學中的Weierstrass-Mandelbrot函數［16］（W-M）被用來表征具有無序、非靜態隨機與自仿射等特性的粗糙表面形貌，本文利用WM 函數對軸承內外圈粗糙形貌進行表征：

式中：x 與z（x）分別為粗糙表面的測量位置與輪廓高度函數，D 與G 分別為粗糙表面輪廓分形維數及尺度常數，常數γ 為大于1 的非整正數，n1與粗糙表面的最低截斷頻率有關.

Majumdar與Tien［17］給出了W-M 函數的功率譜函數：

式中：w 為頻率.

W-M 函數也滿足如下關系

可從方程（4-5）看出W-M 函數是自相似函數，因此，其對采樣長度與測量儀器的分辨率不敏感，且分形參數G 與D 是尺度獨立的參數.對方程（4）兩端取對數，可得

由方程（6）可看出，在雙對數lg（ω）-lgS（ω）坐標系中，W-M 函數的功率譜函數近似為一條直線，該直線斜率為（5-2D），其在縱坐標軸上的截距為2（D-1）lg G-lg 2lnγ.如圖3所示為D＝1.6、G＝9.46×-6m、L＝1μm 且γ ＝1.5 時，利用W-M函數仿真的三維分形粗糙表面形貌及其等高面.

2.2.2 微凸體變形特性 為了確定結合面間的接觸參數，有必要對微凸體的3種變形形式：彈性、彈塑性與塑性變形進行研究.Mikic［18］提出了結合面間接觸熱導的預測模型，該模型假設所有微凸體均發生彈性變形，Yovanovich［19］提出了假設所有微凸體均發生塑性變形的接觸熱導的預測模型.事實上，在載荷的作用下，微凸體會發生彈性、塑性與彈塑性3種變形.因此，Mikic［18］與Yovanovich［19］提出的模型所做的假設并不合理.

圖3 三維分形表面形貌Fig.3 Three-dimensional fractal surface morphology

當微凸體所受最大赫茲接觸壓力達0.6 H 時，微凸體開始發生塑性變形，與微凸體發生塑性變形開始時刻對應的零界接觸面積ac1［20］可表示為

式中：H 為兩接觸材料中較小的微觀硬度，1／E ＝（1-ν）／E1＋（1-ν22）／E2，其中E1，E2，ν1，ν2分別表示兩接觸材料的彈性模量與泊松比.

當微凸體所受最大接觸壓力達到H 時，微凸體發生完全塑性變形，與微凸體發生完全塑性變形時刻對應的零界接觸面積ac2可表示為

接觸面積處于a 與a＋da 間的微凸體數量n（a）［17］可表示為

式中：aL為微凸體最大接觸面積.

所有微凸體的總接觸面積Ar可表示為

當D ≠1.5時，整個接觸面上的載荷F 可表示為

當D＝1.5時，整個接觸面上的載荷F 可表示為

2.2.3 接觸熱導模型 處于接觸狀態的微凸體會發生3種變形形式：彈性、彈塑性與塑性變形，處于3種變形形式的單個微凸體形成的接觸熱導hce、hcep與hcp可分別表示為

結合面間總接觸熱導Hc可表示為

當D ≠1.5時，結合面間總接觸熱導Hc可表示為

當D ＝1.5時，結合面間總接觸熱導Hc可表示為

［18］與Yovanovich［19］的統計學模型的預測精度高，因此，基于幾何-力學-熱分析的綜合接觸熱導模型更加合理.

圖4 接觸熱導預測值與現存模型及實驗數據的比較Fig.4 Comparison of predictions by present model with existing models and experimental data

2.3 熱源生熱功率

2.3.1 電機生熱功率 精確計算電機銅損、鐵損及其他附加損耗非常困難，因此，本文采用效率分析法求解電機生熱功率.電機風損為

式中：d 與L 分別為轉子外徑及長度，h 為定轉子間的間隙厚度，μair為空氣的動力黏度，f 為轉子轉動頻率.

電機效率［21］為

式中：Pbearing為軸承功率損耗，其計算方法將在2.3.2部分詳細討論，Pout為輸出機械能，當主軸無負載運行時，Pout＝0，Pin為電機輸入功率.

電機生熱功率Pmotor可由方程（19）與（20）得到

電機生熱功率在定子與轉子間的分配由電機轉差率fslip與同步頻率fsync共同決定［22］

式中：Protor與Pstator分別為轉子與定子生熱功率.

2.3.2 軸承生熱功率 Harris［23］給出了角接觸球軸承生熱功率P 的表達式

式中：n與M 分別表示軸承轉速與其摩擦力矩.軸承摩擦力矩M 由負荷項Ml與速度項Mv組成

1）軸承外載荷引起的摩擦力矩Ml

外載荷引起的摩擦力矩Ml［23］可表示為

式中：f1為與軸承結構有關的系數，dm為軸承中徑，Fβ由 軸 承 受 力 狀 態 決 定，且 Fβ＝max…（0.9Fa／tanα-0.1Fr，Fr），Fa與Fr分別為軸承軸向與徑向載荷，α為軸承受載后的接觸角.

2）潤滑劑黏性摩擦引起的摩擦力矩Mv

由潤滑劑黏性摩擦引起的摩擦力矩Mv［24］可表示為

式中：vo為潤滑劑的運動黏度，fo與軸承類型及潤滑方式相關.

2.4 對流換熱

對流換熱系數αh

［25］可表示為

式中：hgap為水力半徑，λf為流體導熱系數，Nu 為努塞爾數.努塞爾數可根據不同的對流換熱模型進行計算，如表2所示.

表2 努塞爾數的計算Tab.2 Calculation of Nusselt number

表中：Pr 為普朗特數，Re為雷諾數，lt為冷卻水道長度，w 為水流平均速度，η為流體動力黏度，dk為回路直徑，ds為主軸箱直徑，g 為重力加速度，ξ為摩擦因數，ξ＝1／［0.79ln（Re）-1.64］2.

3 實驗系統和測量原理

3.1 熱特性實驗

實驗系統如圖5所示，主軸系統最高轉速nmax＝30 000r／min，冷卻方式為智能溫度控制，系統內部溫度控制在0～70 ℃，以防止其溫度過高而突然失效，最大熱伸長40μm，X 與Z 方向最大徑向熱變形分別為50與20μm.測量設備及功能如下：溫度位移測量系統是基于美國NI-SCXI1600架構的，利用精密磁吸式溫度傳感器PT100測量主軸系統溫度，其安裝位置如下：前軸承（T1），后軸承（T2），電機（T3），冷卻水套（T4），環境（T5），主軸殼體（T6）.應用高精密電渦流傳感器測量主軸熱漂移，其安裝位置如圖5與6所示，且溫度位移測量為同步測量.

3.2 分形參數辨識

利用OLYMPUS LEXT OLS4000型激光共聚焦顯微鏡對軸承部件粗糙形貌進行測量，如圖7所示.

圖5 熱特性測試實驗現場Fig.5 Thermal characteristic experiment

圖6 位移傳感器安裝方式Fig.6 Spindle three-spot installation diagram

如圖8所示為軸承外圈表面三維粗糙形貌，外圈表面沿X 方向輪廓高度如圖9所示，其中R 表示為粗糙表面輪廓高度，其值作為分形參數的辨識數據源，S 為采樣長度.

圖7 OLYMPUS LEXT OLS4000型激光共聚焦顯微鏡Fig.7 Olympus LEXT OLS4000laser microscope

圖9 軸承外圈X 方向輪廓高度Fig.9 Profile Height of Xdirection

利用如圖9所示軸承外圈表面沿X 方向輪廓的數據，并結合方程（4）求得粗糙表面輪廓的功率譜函數，如圖10所示，并應用最小二乘法對粗糙表面輪廓的功率譜函數進行擬合，可得外圈的分形參數為G＝1.66，D＝7.8×10-14m，同理，可計算軸承內圈表面的分形參數.

4 仿真結果及實驗對比分析

4.1 溫度場分析

圖10 軸承滾道功率譜與其最小二乘擬合曲線Fig.10 Power spectrum and least-squares fitting curves of Xdirection morphology

圖11 主軸轉速分布圖Fig.11 Spindle speed distribution

圖12 主軸系統溫度場云圖Fig.12 Nephogram of temperature field distribution when thermal equilibrium

為了模擬主軸實際加工狀態，設計主軸轉速分布如圖11所示，達到熱平衡狀態時，主軸系統的溫度場分布如圖12所示.主軸系統內部溫度場不均勻且梯度明顯，最高溫度θ＝55.37 ℃（轉子），t為時間，電機定子生熱功率（1 635 W）是主軸系統所有部件中最大的，但其溫度卻不是最高的，這是因為定子與冷卻系統間對流換熱明顯，顯著降低了定子的溫度.冷卻水套帶走了定子大部分熱量，因此，冷卻水套溫度升高了，但是，冷卻水套與主軸箱接觸，有相當一部分熱量由冷卻水套傳遞進入主軸箱，主軸箱溫度也相應升高.

前軸承溫度最高可達41.2℃，且前軸承較后軸承更加靠近電機轉子，因此，由轉子經軸進入前軸承的熱量多于進入后軸承的熱量，因此，前軸承溫度高于后軸承.軸承與軸承座，軸承與軸間溫度梯度明顯，這是因為這些位置處存在接觸熱導使得軸承熱量擴散到軸承座與軸更加困難，但是，在主軸運轉過程中，軸承仍在不斷產生熱量，因此，由于接觸熱導的阻力作用使得軸承產生的熱量在其自身上發生了累積，導致其局部溫度過高.將前軸承、后軸承、電機定子與冷卻水套溫度的仿真值與實驗值進行比較，如圖13所示，可以看出，考慮結合面間接觸熱導的有限元模型能準確預測測點溫度，且測點溫度在初期迅速上升，然后逐漸飽和到最終溫度.

圖13 溫度測量值與仿真值比較Fig.13 Comparison of the temperature histories，simulated，experimental

為了驗證結合面間接觸熱導對仿真結果的影響，對忽略接觸熱導的有限元模型進行仿真，并將前軸承與后軸承溫度分別與測量值和考慮接觸熱導的有限元模型進行比較，如圖14所示，可以看出，考慮接觸熱導的有限元模型的仿真精度遠高于不考慮接觸熱導的有限元模型.且考慮接觸熱導時，前后軸承溫度分別高于不考慮接觸熱導時的前后軸承的溫度.原因如下：結合面間接觸熱導的存在使得由一個物體向另一個物體的傳熱受到阻礙，即考慮接觸熱導時，由于軸承內圈與軸，軸承外圈與軸承座間接觸熱導的存在，使得軸承（熱源）產生的熱量傳遞到軸與軸承座變得更加困難，即軸承結合面間接觸熱導的阻礙作用限制了熱量從軸承傳向軸與軸承座，但是，在主軸持續運轉過程中，軸承（熱源）仍然在不斷的產生熱量，因此，熱量在軸承上產生了積累，使得軸承局部溫度過高，故當考慮接觸熱導時，軸承溫度比不考慮接觸熱導時的軸承溫度高，如圖14所示；而對于軸與軸承座，考慮接觸熱導時，由于接觸熱導對結合面間傳熱的阻礙作用，使得由軸承傳遞到軸與軸承座上的熱量比不考慮接觸熱導時的熱量少，因此，當考慮接觸熱導時，軸與軸承座溫度比不考慮接觸熱導時的溫度低.并且結合面間的接觸熱導越小，結合面兩側溫差越大.

圖14 溫度測量值與仿真值比較Fig.14 Comparison of temperature experimental data with the simulation results

4.2 熱變形分析

如圖15所示為主軸系統熱平衡時的熱變形云圖（放大250倍），可以看出，主軸系統最大變形量為60.8μm（Y＋端），主軸X、Y 與Z 方向的最大變形量分別為58.0、10.5與8.0μm，前軸承溫度高于后軸承導致主軸左側熱變形較主軸右側對應位置的熱變形大.

圖15 主軸系統熱變形場云圖Fig.15 Nephogram of thermal deformation of spindle system

為了研究主軸系統熱變形對加工精度的影響，對3個方向的熱變形仿真值與測量值進行比較，如圖16所示，可以看出X 方向熱變形最明顯，其對加工精度有明顯影響，Y 方向熱變形也較大，即主軸會沿軸向伸長，對加工精度也有很大影響.

熱變形測量值、考慮結合面間接觸熱導與不考慮結合面間接觸熱導的有限元模型的熱變形的比較如圖17所示，可看出，考慮結合面間接觸熱導的熱變形的精度比不考慮結合面間接觸熱導的熱變形精度較高，并且考慮結合面接觸熱導的有限元模型X，Y 與Z3個方向熱變形均大于不考慮接觸熱導3個方向對應的熱變形.原因如下：當考慮接觸熱導時，由軸承傳遞到軸的熱量比不考慮接觸熱導時的熱量少，因此，對于軸而言，其熱變形也會相應減小，如圖17所示.

圖16 主軸X，Y 與Z 方向熱變形測量值與仿真值比較Fig.16 Comparison between thermal errors of measurement data with simulation results of X，Yand Z directions

圖17 主軸熱變形測量值與仿真值比較Fig.17 Comparison of thermal deformations of measurement data with the simulation results with and without TCC

因此，在考慮接觸熱導情況下，軸承產生的熱量在其自身產生了積累，且由軸承進入主軸與軸承座的熱量少于不考慮接觸熱導時的熱量，因此，考慮接觸熱導情況下軸承溫升的仿真值大于考慮接觸熱導的情況，不考慮接觸熱導情況下其熱變形值大于考慮接觸熱導的情況.

5 結 語

本文提出了高速主軸熱特性分析的系統化建模方法.建立了基于幾何-力學-熱的綜合接觸熱導預測模型，該模型考慮到粗糙表面形貌與微凸體接觸變形對結合面間接觸熱導的影響，且有效避免了統計學方法的不準確性與實驗測量法通用性不強的缺陷，分析了電機與軸承生熱功率求解方法，給出了主軸部件與流體間的對流換熱系數的計算方法，進行熱平衡實驗驗證建模方法的有效性，分析了主軸溫度分布規律與熱變形變化，結果表明，本文建立的主軸熱特性分析模型能夠準確預測主軸溫度場與熱變形，且仿真精度遠高于忽略接觸熱導的有限元模型.雖然本文模型是針對高速主軸建立的，但模型不失一般性，其分析方法可為機床整機或進給系統熱特性仿真提供參考.結合面間接觸熱導模型適用于所有固定結合面間傳熱分析.

（

）：

［1］MAYR J，JEDRZEJEWSKI J，UHLMANN E，et al.Thermal issues in machine tools［J］.CIRP Annals-Manufacturing Technology，2012，61（2）：771-791.

［2］ABELE E，ALTINTAS Y，BRECHER C.Machine tool spindle units［J］.CIRP Annals-Manufacturing Technology，2010，59（2）：781-802.

［3］楊軍，梅雪松，趙亮，等.基于模糊聚類測點優化與向量機的坐標鏜床熱誤差建模［J］.上海交通大學學報，2014，8：1175-1182.YANG Jun，MEI Xue-song，ZHAO Liang，et al.Thermal error modeling of jig-borer based on fuzzy clustering and support vector machine［J］.Journal of Shanghai Jiaotong University，2014，8：1175-182.

［4］YANG J，SHI H，FENG B，et al.Applying neural network based on fuzzy cluster pre-processing to thermal error modeling for coordinate boring machine［J］.Procedia CIRP，2014，17：698-703.

［5］楊軍，梅雪松，馮斌，等.時序分析在電主軸熱誤差建模中的應用［J］.計算機集成制造系統，2014（1）：10-18.YANG Jun，MEI Xue-song，FENG Bin，et al.Application of time series in the thermal error modeling of highspeed spindle system［J］.Computer integrated manufac-turing system，2014（1）：10-18.

［6］楊軍，施虎，梅雪松，等.雙驅伺服進給系統熱誤差的試驗測量與預測模型構建［J］.西安交通大學學報，2013，47（11）：53-59.YANG Jun，SHI Hu，MEI Xue-song，et al.Prediction model construction and experimental measurement of the thermal error of dual servo feed drive system［J］.Journal of Xi'an Jiaotong University，2013，47（11）：53-59.

［7］BOSSMANNS B，TU J F A thermal model for high speed motorized spindles［J］.International Journal of Machine Tools and Manufacture，1999，39（9）：1345-1366.

［8］ZHAO H，YANG J，Shen J.Simulation of thermal behavior of a CNC machine tool spindle［J］.International Journal of Machine Tools and Manufacture，2007，47（6）：1003-1010.

［9］CREIGHTON E，HONEGGER A，TULSIAN A，et al.Analysis of thermal errors in a high-speed micromilling spindle［J］.International Journal of Machine Tools and Manufacture，2010，50（4）：386-393.

［10］KYUNG CHOI J，GIL LEE D.Thermal characteristics of the spindle bearing system with a gear located on the bearing span［J］.International Journal of Machine Tools and Manufacture，1998，38（9）：1017-1030.

［11］HOLKUP T，CAO H，KOLáR P，et al.Thermo-mechanical model of spindles［J］.CIRP Annals-Manufacturing Technology，2010，59（1）：365-368.

［12］KOLAR P，HOLKUP T.Prediction of machine tool spindle's dynamics based on a thermo-mechanical model［J］.MM Science Journal，2010：166-171.

［13］GREENWOOD J A，WILLIAMSON J B P.Contact of nominally flat surfaces［J］.Proceedings of the Royal Society of London.Series A.Mathematical and Physical Sciences，1966，295（1442）：300-319.

［14］COOPER M G，MIKIC B B，YOVANOVICH M M.Thermal contact conductance［J］.International Journal of heat and mass transfer，1969，12（3）：279-300.

［15］MCWAID T H，MARSCHALL E.A comparison of elastic and plastic contact models for the prediction of thermal contact conductance ［J］. W?rme-und Stoffübertragung，1993，28（8）：441-448.

［16］BERRY M V，LEWIS Z V.On the Weierstrass-Mandelbrot fractal function［J］.Proceedings of the Royal Society of London.A.Mathematical and Physical Sciences，1980，370（1743）：459-484.

［17］MAJUMDAR A，TIEN C L.Fractal network model for contact conductance［J］.Journal of Heat Transfer，1991，113（3）：516-525.

［18］MIKIC B B.Thermal contact conductance；theoretical considerations［J］.International Journal of Heat and Mass Transfer，1974，17（2）：205-214.

［19］YOVANOVICH M M.Thermal contact correlations［J］.AIAA Paper，1982，81：83-95.

［20］MAJUMDAR A，BHUSHAN B.Fractal model of elastic-plastic contact between rough surfaces［J］.Journal of Tribology，1991，113（1）：1-11.

［21］KYLANDER G.Thermal modelling of small cage induction motors［M］.Goteborg：Chalmers University of Technology，1995：1-113.

［22］AVALLONE E A，BAUMEISTER T，SADEGH A.Marks'Standard Handbook For Mechanical Engineers：Standard Handbook for Mechanical Engineers ［M］.New York：Mcgraw-Hill Professional，2006：1-386.

［23］HARRIS T A.Rolling bearing analysis［M］.New York：Wiley，1984：1-1013.

［24］JONES A B.A general theory for elastically constrained ball and radial roller bearings under arbitrary load and speed conditions［J］.Journal of Fluids Engineering，1960，82（2）：309-320.

［25］XU M，JIANG S，CAI Y.An improved thermal model for machine tool bearings［J］.International Journal of Machine Tools and Manufacture，2007，47（1）：53-62.