0  引言
   
機床的熱特性問題一直阻礙加工精度的進一步提高，熱控制技術是研制高精度機床工作中急需解決的關鍵技術之一。由于機床在正常工作條件下，熱源(電機、軸承等)產生的熱量傳給機床部件，使各部位產生一定的溫升，同時零部件產生熱變形。由于零件結構、形狀等尺寸不一，將產生彎曲、扭曲、拉伸等變形即熱位移，這樣造成主軸、工作臺、刀具等部位發生一定的相對位移，從而影響加工精度，使加工精度大大降低。機床主軸等重要零部件熱變形誤差是引起機床變形誤差的重要因素，因此在機床設計時，減少主軸系統的發熱和加強散熱，降低溫升，通過對主軸系統進行熱特性分析是保證機床的加工精度的重要措施。
   
采用有限元數值分析軟件可以較好地分析這個問題。ANSYS Workbench 12是ANSYS求解實際問題的新一代產品。除了可以解決靜力學、動力學、非線性接觸分析以外還支持直接耦合場分析、如熱-結構、熱-電、熱-磁耦合等。熱-結構耦合分析是耦合物理場分析中比較常見和重要的一種，本文以柳州市數控機床研究所研發的CK6136數控車床主軸系統為研究對象，建立主軸系統熱特性分析的有限元模型，在有限元分析軟件ANSYS Workbench 12中對其進行熱穩態分析，到主軸的溫升和溫度場分布情況，以求得的溫度場為依據進行熱-結構耦合分析計算得到其熱變形，同時計算主軸不同轉速條件下對主軸系統產生的溫升及熱變形的影響，最后得出相關結論。
   
1 熱分析基本理論
   
對于穩定熱分析，如果系統的凈熱流率為0，即流入系統的熱量加上系統自身產生的熱量等于流出系統的熱量：q流入+ q生成－ q流出= 0，在穩態熱分析中任一節點的溫度不隨時間變化。穩態熱分析的能量平衡方程為：
    
式(1)中：

[K]為傳導矩陣，包含導熱系數，對流系數及輻射率和形狀系數；
{T}為節點溫度向量；
{Q}為節點熱流率向量，包含熱生成。
   
根據參考文獻，由邊界條件和載荷結合熱應力分析可以推導出整體剛度矩陣[K]和載荷向量{P}：
    
式(2)、式(3)中[K]e為單元剛度矩陣，{P}e為單元節點載荷，{P}eT為單元節點熱載荷。
   
根據位移模式，由平衡條件、變分原理及胡克定律[K]·{δ}={P}可以推導出節點位移{δ}，由單元節點位移{δ}eT及溫升{ΔT}e，計算獲得總應變{ε}及熱變形{ε}T，最終熱應力表達式為：
    
式(4)中[D]為彈性矩陣。
   
2 主軸系統有限元模型
   
本文研究對象CK6136 數控車床主軸系統主要由主軸、前后角接觸球軸承、床頭箱、封蓋前后端蓋等部件組成，對于復雜結構在進行CAD建模時需要對模型進行適當的簡化。例如模型中的圓角和倒角以及對分析結果不造成很大影響的小孔就應該盡量省略，還有將鍵槽、螺紋孔、擋油孔均按實體處理等，在三維軟件UG6.0建立好主軸系統結構簡圖(半剖視圖)如圖1所示，同時設置好UG6.0和ANSYS Workbench 12的無縫連接。

3 機床熱源及邊界條件計算
   
3.1 熱源分析
   
對于數控車床主軸系統熱源分析本文主要考慮前、后軸承的摩擦發熱，而忽略電動機發熱、切削熱發熱對主軸系統的影響，根據參考文獻主軸軸承熱源的計算公式如下：
   
式(5)中M為軸承摩擦力矩，指由軸承滾動摩擦、滑動摩擦和潤滑劑摩擦的總和產生的阻滯軸承運轉的阻力矩；n為軸承轉速。M的計算公式如下：
    
式(6)中M1為和軸承載荷大小、滾動體和滾道間接觸彈性變形量及滑動摩擦有關的摩擦力矩分量，M1的計算公式如下：
    
M2為與速度有關的選項，它反映潤滑劑的流體動力消耗，M2的計算公式如下：
   
(1)當vn＞3. 33×10-5m2·r /s2 時：
     
(1)當vn< 3. 33×10-5m2·r /s2 時：
   
上式(7)~式(9)中，f1為軸承類型和所受負荷有關的系數；p1為軸承摩擦力矩的計算負荷，f1和p1可查閱相關資料查得；f0為與軸承類型和潤滑方式有關的常數；v為潤滑劑的運動粘度；dm為軸承中徑；n為軸承轉速。
   
3.2 換熱系數的計算
   
在空轉且不考慮冷卻液的條件下主軸部件對流則是主軸部件表面與空氣之間的對流，根據參考文獻由努謝爾特準則，換熱系數計算公式為：

式中，Re為雷諾數，Pr為普朗特數。

(2)對于床頭箱
 

式中，C、n為常數由相關參考文獻查出；Gr為格拉曉夫準數；g為重力加速度；β為體膨脹系數：L為形狀特征尺寸；v為運動粘度；ΔT為流體與壁面溫差。

3.3 材料屬性設置
   
在ANSYS Workbench 12中重新添加床頭箱，軸承和主軸三種材料，其材料屬性見表1。

4 主軸系統熱特性有限元分析
   
4.1 熱穩態分析
   
將模型導人有限元分析軟件ANSYS Workbench 12中，對主軸系統模型進行四面體網格劃分，網格劃分后得到1807105個節點，1184866個單元，網格劃分模型如圖2所示。

本文研究以主軸n=3000r/min機床最高工作轉速為計算依據，空轉3小時達到熱平衡狀態，根據前述計算公式，可以得到前軸承(靠近刀具端)發熱量為56.145W，后軸承(遠離刀具端)發熱量為48.73W；主軸內表面傳熱系數為213.667W/m2℃ 主軸外表面的傳熱系數為107.189W/m2℃，床頭箱內其他表面傳熱系數為9W/m2℃。將以上邊界條件約束到有限元模型上，得到主軸系統的溫度場分布如圖3。

由上圖可知，主軸最高溫度發生在主軸前端軸承安裝處達52.09℃，后端軸承安裝處溫度達48.15℃，箱體溫度較低且溫度分布不均勻，前后軸承溫升不一致，最大溫升為32℃(空氣溫度假設為20℃情況下)，從溫度場分布可知由于在軸承安裝處溫升不同將導致主軸系統發生一定的熱膨脹，從而產生熱變形。

為了驗證模型計算的正確性，采用非接觸式紅外溫度測量儀對主軸系統作了溫度測試實驗分析，具體操作如下：在n=3000r/min，機床正常空轉3小時后，測溫點分布在主軸前后端軸承附近，取主軸圓周等分位置四點多次測量取平均值，同時還測量了箱體表面的溫度。實驗得到前端軸承附近溫度測量值為48.9℃，理論分析值誤差為6.5%；后端軸承附近溫度測量值為46.5℃，理論分析值誤差為3.5%；箱體側表面溫度測量值為43.1℃理論分析值誤差為5.6%。由以上分析得到的結果誤差基本上控制在7%以內，可以認為由模型得到的計算結果是可靠的，與實際情況相符。
   
4.2 熱-結構耦合分析
   
在有限元分析軟件ANSYS Workbench12中利用直接耦合順序耦合解方法，將得到的溫度場的溫度分布作為熱載荷加載到模型上，同時約束床頭箱體底部與床身接觸面的所有自由度，計算得到主軸系統的熱變形如圖4所示。

由上圖可知，機床主軸系統最大變形發生在床頭箱體底部一小范圍，達54.537μm，主軸系統主軸前端徑向位移為8.9μm、綜合位移為30.298μm；主軸后端徑向位移為12.2μm，綜合位移為42.417μm，有一定的變形，但還是符合加工精度要求。
   
4.3 主軸轉速變化對主軸系統溫升及熱變形的影響
   
由式(5)～式(9)，計算出不同主軸轉速條件下前后軸承的發熱量見表2。在有限元分析軟件中重新計算得到主軸轉速變化對主軸系統產生的溫升及熱變形的影響如圖5、圖6所示。