微小間隙空氣靜壓主軸角剛度分析與實驗**

周 亮 王寶瑞 陽 紅

(中國工程物理研究院機械制造工藝研究所，四川 綿陽 621900)

微小間隙空氣靜壓主軸角剛度分析與實驗**

周 亮 王寶瑞 陽 紅

(中國工程物理研究院機械制造工藝研究所，四川 綿陽 621900)

通過綜合氣體運動方程、連續性方程和狀態方程，建立小孔節流氣體潤滑主軸數學模型。采用有限單元法、超松弛迭代法和比例分割法對雷諾方程進行離散求解，解決了微小間隙空氣靜壓主軸跨尺度微偏轉的氣膜壓力分布數值計算難題。進行靜態角剛度測試實驗，驗證了所提方法的有效性。

空氣靜壓主軸；角剛度；有限元分析；小孔節流

靜壓主軸是超精密機床實現精密回轉運動的關鍵部件，它采用流體潤滑技術，借助于氣膜或液膜來分開運動件與支承件，以減小摩擦，提高運動穩定性。在超精密加工機床中，靜壓主軸是電動機與刀具或工件裝夾機構聯系的紐帶，對機床的加工精度起著決定性影響。

目前在靜壓主軸方面的研究主要集中在剛度(徑向剛度、軸向剛度、前端剛度)、回轉精度、動態特性以及熱穩定性等方面。如姚俊等介紹了各種機床主軸回轉精度測量方法及其原理以及測量數據的處理和誤差分析方法[1]。彭萬歡等采用兩點測量法深入分析了超精密空氣靜壓主軸回轉誤差測試過程中偏心的影響和作用原理[2]。陳東菊等研究了氣體稀薄效應對軸承承載力和剛度的影響，為微尺度下軸承性能的研究提供了一定的理論依據[3]。在主軸性能研究中，其軸線抵抗偏擺的能力——角剛度是一個很重要的指標。軸線偏擺不僅會影響主軸的回轉精度；在主軸運動過程中，如果角剛度較小，還可能產生擺振。然而，國內對角剛度的研究還不夠深入。侯國安等在研究靜壓主軸動態特性時討論了水平方向角剛度的理論計算方法，但未考慮實際工況下切削力以及重力作用存在的初始偏心對角剛度的影響，沒有給出詳細的驗證說明[4]。安晨輝等通過旋轉歐拉方程研究了主軸穩定性，具有借鑒意義，但也沒有驗證其方法的合理性[5]。楊福興[6]、夏歡[7]以及徐灝[8]等人結合軸承尺寸，分別給出了通過徑向剛度和軸向剛度計算角剛度的經驗公式，齊乃明等[9]推導了雙列徑向軸承在無偏心狀態下的角剛度計算；但此類公式不含角度參數，不能計算變角度情況下的角剛度。

目前，對靜壓主軸角剛度的研究，通常僅計算了無重力作用無偏心狀態下的角剛度，與實際工況差異較大，也沒有對角剛度進行實驗測試驗證，沒有對角剛度的形成機理進行深入系統的研究。本文在氣體運動方程、連續性方程和狀態方程的基礎上，引入重力和初始偏心距，建立起實際工況下小孔節流氣體靜壓主軸徑向軸承和止推軸承數學模型，采用有限單元法、超松弛迭代法和比例分割法進行離散、求解和計算氣膜內壓力分布，從而解決微小氣膜間隙下微偏轉跨尺度的數值計算問題。最后通過實驗測試靜態條件下靜壓主軸的角剛度，驗證數值計算的準確性。

1 空氣靜壓主軸角剛度建模

1.1 空氣靜壓主軸角剛度定義

空氣靜壓主軸依靠氣體靜壓效果實現支承和限制主軸轉子位置和姿態。徑向剛度和軸向剛度評價了轉子X、Y方向的位置穩定性，與此類似，角剛度評價了轉子姿態穩定性。因此，與剛度定義類似，角剛度為靜壓主軸力矩與偏擺角度的比值，計算式為

kθ=M/θ

(1)

式中：kθ為角剛度；M為主軸總力矩；θ為主軸偏轉角。

主軸通常由徑向軸承和止推軸承組成。如圖1所示，為某T型結構的空氣靜壓主軸模型。由于T型主軸結構非對稱，其角剛度采用經驗公式計算時誤差會顯著增大。因此，需采用有限元法來計算靜壓主軸的角剛度。

1.2 空氣靜壓主軸角剛度理論建模

如圖2所示，當T型軸繞過A點垂直紙面軸線發生微小偏轉時，主軸受力可以簡化為徑向合力F1和止推力偶M(F2，F2’)。徑向軸承力矩和止推軸承力矩可以分別計算，同時止推軸承產生的力矩可以投影為繞軸承圓盤中心的轉動。主軸總力矩為

M=M1+M2=F1L+2F2Rsinα=F1L+2F2r

(2)

1.2.1 空氣靜壓主軸徑向軸承建模

對于徑向軸承，其氣膜厚度遠小于其他兩個方向的尺寸，通常情況下，厚度方向的壓力和速度變化可以忽略；同時，當無轉動或者轉速很低時，雷諾數很小，與之相關的項可以忽略。經過簡化，徑向軸承的Reynolds方程無量綱化形式[10]為

(3)

(4)

式中：η為空氣粘滯系數；l、hm分別表示無量綱化對應參考長度、氣膜厚度；p0、pa、ρa分別表示供氣壓力、大氣壓力和大氣密度；mr為單位時間流過節流孔的質量流量。

將氣膜沿垂直方向切分展開，則壓力分布近似為平面問題。將氣膜研究域劃分為三角形單元體集，如圖3所示，則會產生對稱邊界、大氣邊界和節流孔邊界三種邊界條件。

每個單元體內任意點壓方都可由所在單元體的三節點壓方表示，設插值函數為

f=Nifi+Njfj+Ntft

(5)

按照變分原理對雷諾方程進行處理，同時將式(4)和(5)代入雷諾方程式(3)，對每個節點可以列寫出一個方程，則整個研究域可以寫出所有未知節點的方程，構成一個位置函數的矩陣方程：

KF=T

(6)

式中：F為待求節點壓力函數，為n×1維列陣，K為n×n維剛度矩陣，是Ni、Nj、Nt函數的矩陣形式；T為n×1維列陣，是與邊界條件和節流孔流量有關的非線性列陣。涉及到的流量項可按下式計算：

(7)

(8)

式中：Ar為節流面積；φ為流量系數，一般可取0.8； pd為節流孔出口壓力； k為氣體比熱比，空氣一般為1.4；其余同上。

當主軸發生偏轉后，氣膜厚度h也是位置的函數，如圖4所示。

h=h0+e·cosγ=h0+Lsinθ·cosγ+e0·cosγ

(9)

式中：h0為氣膜平均厚度；γ為單元體與豎直方向夾角；e為偏心距；e0為初始偏心距；L為單元體至偏轉中心距離；θ為主軸偏轉角。

聯立式(6)～(9)可以求解氣膜內壓力分布，但式(8)是非線性矩陣，故采用超松弛迭代法求解方程。若迭代后精度不滿足要求，則采用比例分割法計算節流孔處第i+1次迭代壓方初始值，如公式(10)。

(10)

求得氣膜壓力分布后，整個主軸的承載W1和力矩M1計算式為

(11)

(12)

式中：m為單元體總個數。

1.2.2 空氣靜壓主軸止推軸承理論建模

與徑向軸承類似，止推軸承氣膜厚度方向的壓力和速度變化可以忽略。靜態時，止推軸承的Reynolds方程無量綱化形式[10]為

(13)

(14)

式(14)與式(3)十分相似，對式(14)的處理與式(3)相同。止推軸承通常為圓盤形，止推盤偏轉時，氣膜厚度的變化如圖5所示，(r，β)單元體上氣膜厚度h為

h=h0+s·sinθ=h0+rcosβ·sinθ

(15)

式中：s為單元體與水平軸線的距離；β為單元體與豎直方向夾角：r為單元體所在圓半徑；h0為止推氣膜平均厚度；θ為主軸偏轉角。

同樣，用有限單元法求解止推軸承的承載W2與力矩M2的計算與徑向軸承類似：

(16)

(17)

式中：m為單元體總個數。

2 空氣靜壓主軸角剛度數值計算

2.1 空氣靜壓主軸角剛度計算流程

基于以上分析，我們在MATLAB 14.0基礎上開發了角剛度計算程序，程序流程如圖6所示。通過改變參數輸入的大小，可以實現不同尺寸主軸徑向軸承和止推軸承模型的快速建模和計算。

2.2 空氣靜壓主軸角剛度計算結果

以某T型靜壓主軸為例，該主軸由兩個徑向軸承和一個止推軸承組成。徑向軸承直徑D=100 mm，寬度B=100 mm，氣膜厚度h1=12 μm，兩軸承間距L’=50 mm，雙圈節流孔，每圈12孔，節流孔直徑d=0.15 mm；止推軸承內徑da=104 mm，外徑db=228 mm，氣膜厚度h2=10 μm，雙圈節流孔，第一圈圓半徑r1=68 mm，12孔，第二圈圓半徑r2=98 mm，24孔，節流孔直徑d=0.15 mm；整個轉子部分約重28.1 kg。

基于編寫的MATLAB程序，可以分別計算單個徑向和止推軸承的承載和力矩。計算中，供氣壓力為0.5 MPa，出口壓力為0.1 MPa，徑向軸承劃分為48×72個三角形單元，1 800個節點；止推軸承劃分為30×128個三角形單元，2 304個節點。首先，當偏角θ0=0，計算出主軸在重力作用下的初始偏心距為e0；然后計算一定偏角θ下的徑向和止推氣膜上的角剛度kθ1、kθ2。對于上述主軸，由該程序計算可以得到初始偏心距e0≈1.1 μm，在偏角θ=1×10-5rad時，靜壓主軸徑向軸承角剛度為kθ1=1.2×105N·m/rad，止推軸承的角剛度為kθ2=16.4×105N·m/rad，故整個靜壓主軸的角剛度kθ=kθ1+kθ2=17.6×105N·m/rad。

2.3 經驗公式計算

上述主軸設計的徑向剛度為405 N/μm，軸向剛度為507 N/μm。文獻[8-9]給出了通過徑向剛度和止推剛度計算角剛度公式：

徑向軸承的角剛度kθ1為

kθ1=k1B2/16

(18)

式中：k1為徑向剛度，B為徑向軸承寬度。

止推軸承的角剛度kθ2為

kθ2=k2((da+db)/2)2/8

(19)

式中：k2為軸向剛度；da、db為止推盤內外徑。

由式(18)和式(19)可以求得該主軸的徑向軸承角剛度kθ1=2.5×105N·m/rad，止推軸承角剛度kθ2=17.5×105N·m/rad，整個主軸角剛度kθ=kθ1+kθ2=20.0×105N·m/rad。

3 空氣靜壓主軸角剛度實驗測量

3.1 力矩—位移法測量角剛度

針對上述靜壓主軸，為了驗證角剛度理論，設計了角剛度力矩—位移測量方法。實驗原理如圖7所示。力矩M=FL1，角度θ≈(s1+s2)/L2，角剛度kθ=M/θ。被測主軸和測量裝置如圖8所示。由于主軸轉子部分在主軸座內部，無法直接加載和測量，故力矩加載通過過渡接頭連接，位移測點分布在連接轉子的測量片上。測量中保證兩力大小相等，方向相反。通過實驗，我們測得了該主軸在微小偏轉時力矩與角度的關系(見表1)，繪制了力矩—角度曲線，如圖9所示。

采用最小二乘法通過MATLAB對實驗數據進行一次項擬合，可以得到力矩—角度關系式為：

y=14.4x-1.091

(20)

對上式求導，可得該主軸在所測范圍內的角剛度約等于14.4 N·m/rad。

表1 靜壓主軸力矩與角度測量

力矩/(N·m)01．0111．5022．1162．813角度/(×10-5rad)00．0740．1220．2720．278力矩/(N·m)3．4944．2574．9145．6256．144角度/(×10-5rad)0．3890．4260．4870．5150．544力矩/(N·m)7．4808．2999．20710．136角度/(×10-5rad)0．5610．6110．6300．728

3.2 理論計算與實驗對比分析

對比實驗結果與計算分析的數值，實驗與數值計算結果相近。并且，有限單元法的計算結果比經驗公式更接近實驗測量值。同時，通過有限單元法可以計算任意偏角下的角剛度，彌補了經驗公式的不足。

數值計算比實驗結果偏大，主要有以下幾方面的原因。一是數值計算將氣膜內氣體的流動均視為層流，而實際氣流中存在紊流，此會導致數值計算結果偏大；二是主軸受力會產生形變，位移測量受形變誤差影響使計算的偏角增大，此會導致測量的角剛度偏?。蝗侵鬏S結構尺寸在供氣壓力作用下也會改變，如小孔直徑在較大供氣壓力作用下可能增大，四是由于傳感器和儀器的粘滯特性，測量過程中存在間隙誤差。角剛度測量的誤差來源較復雜，目前測量的角剛度是各種因素作用下的綜合角剛度。因此，在下一步的研究中，角剛度實驗測量的誤差分離將會是研究的重點。

4 結語

(1)采用有限單元法編寫的角剛度數值計算方法與實驗結果比較吻合，運用該算法能夠實現各種尺寸結構的靜壓主軸角剛度計算。

(2)有限單元法比經驗公式計算結果更接近實驗測量值；并且有限單元法克服了經驗公式無法確定角剛度與偏轉角度的對應關系，具有更高的靈活性和更廣的適應性。

(3)角剛度的測量存在復雜的誤差來源，高精度的測量需要進行嚴格的誤差分離。

[1]姚俊，王平. 主軸回轉精度測量方法[J]. 制造技術與機床，2011(12)：176-180.

[2]彭萬歡，趙午云，陶繼忠，等. 超精密空氣靜壓主軸徑向回轉誤差的測試研究[J]. 制造技術與機床，2008(10)：55-58.

[3]陳東菊，周帥，楊智，等. 稀薄效應對空氣靜壓止推軸承性能影響[J]. 四川大學學報，2016，48(1)：194-199.

[4]侯國安. 流體靜壓支承對超精密金剛石車床動態特性影響的研究[D]. 哈爾濱：哈爾濱工業大學，2013.

[5]安晨輝，許喬，張飛虎. 空氣靜壓主軸的壓強分布與動態特性分析[J]. 納米技術與精密工程，2009，7(5)：459-468.

[6]楊福興，董申. 超精密空氣靜壓主軸靜態性能的數值分析[J]. 機械設計與研究，2004，20(2)：54-56.

[7]夏歡，吉方，陶繼忠. 空氣彈簧調平系統研究[J]. 現代制造工程，2012(11)： 104-107.

[8]徐灝，邱宣懷，蔡春源，等. 機械設計手冊[M]. 北京：機械工業出版社，1995.

[9]Qi Naiming，Li Zhongying，Yang Guojun，et al. Study on high stiffness gas bearing for precision centrifuger[J]. Journal of HarBin Institute of Technology，1998，5(4)：47-50.

[10]劉暾，劉育華，陳世杰. 靜壓氣體潤滑[M]. 哈爾濱：哈爾濱工業大學出版社，1990.

如果您想發表對本文的看法，請將文章編號填入讀者意見調查表中的相應位置。

Analysis and test of angular stiffness on aerostatic bearing

ZHOU Liang, WANG Baorui, YANG Hong

(Institute of Mechanical Manufacturing Technology, China Academy of Engineering Physics, Mianyang 621900, CHN)

The paper synthesizes gas motion equation, continuity equation and state equation, to establish mathematical model of orifice gas lubricated journal bearing. The Reynolds equation is discretized by the finite element method, then solved by the super relaxation iteration method and the proportional division method. The paper solves numerical problems of cross scale and micro deflection on small gas film clearance of aerostatic bearing spindle. Static angular stiffness experiments were carried out, which proves the validity of the proposed method.

aerostatic bearing spindle; angular stiffness; finite element analysis; orifice restriction

*國家自然科學基金(51305413)；中物院科技專項資助(9120602)

TH133.35

A

10.19287/j.cnki.1005-2402.2016.12.006

周亮，男，1990年生，碩士研究生，主要研究方向為精密加工及測量。

(編輯 汪 藝)

2016-09-26)

161214