王學(xué)敏，白紅宇，莊明

(中國科學(xué)院等離子體物理研究所，安徽合肥230031)

由于氣體潤滑介質(zhì)具有摩擦功率小、無污染、運動精度高等優(yōu)點，特別是近年來隨著核能和航空航天技術(shù)的發(fā)展，其應(yīng)用前景廣闊，通過求解Reynolds方程準(zhǔn)確預(yù)測軸承性能是氣體軸承應(yīng)用的關(guān)鍵.目前，已經(jīng)有多種數(shù)值方法應(yīng)用于對Reynolds方程的求解［1］，其中由于有限差分法編程序靈活、成熟，其應(yīng)用最為廣泛.國內(nèi)外雖然有多篇關(guān)于利用有限差分法求氣體軸承性能的文章［2-4］，但給出的大都是未考慮偏位角影響的靜態(tài)算例，在求解過程中對一些影響計算精度的因素也并未加以討論，本文對用有限差分求解軸頸軸承性能過程中若干影響計算精度的問題進(jìn)行了細(xì)致討論，并得出了若干可以提高計算結(jié)果精度的結(jié)論.

1 Reynolds方程無量綱化及邊界條件

很多文獻(xiàn)中都有詳細(xì)的Reynolds方程推導(dǎo)過程［5-6］，完整的非定?？蓧嚎s雷諾方程如下:

以上代入式(1)，得到穩(wěn)態(tài)無量綱化Reynolds方程式如下:

式中:H=1+εcos θ，ε 為偏心率，式(2)的無量綱邊界條件為:

2)氣膜被剪開的兩邊界上:P(θ=0)=P(θ =2π)，即 P1，j=PΠ，j，其中 Π 為最大網(wǎng)格節(jié)點數(shù)，dP/dx(θ=0)=dP/dx(θ=2π)，即(P2，j-P1，j)/ΔX=(PΠ，j-PΠ-1，j)/ΔX.

2 Reynolds方程離散化

由于氣膜厚度相對于其他幾何尺寸小103～104倍，因此可以略去圓柱表面率的影響，為了便于程序編寫，一般都是把圓柱軸承的氣膜展開成如圖1的一個平面［7-8］，將軸承面按兩供氣孔的中心線劃分為N個區(qū)域，其中N為供氣孔數(shù)(本文以雙排8個供氣孔為例).

圖1 氣膜劃分區(qū)域示意圖Fig.1 Computational domain of journal gas bearings

采用具有二階精度的中心差分格式的超松弛迭代法(SOR)［9］對式(2)進(jìn)行離散，得到離散化Reynolds方程式為

單個區(qū)域流出的質(zhì)量流量包括:

軸承承載能力計算為

式中，Wn和Wt分別表示氣膜內(nèi)壓力在偏心方向上對面積積分和與其成直角方向上對面積積分.

3 有限差分法求解過程若干問題

為了比較算法的有效性，本文利用具有解析解的ε=0時軸頸軸承模型，其模型如圖2，具體參數(shù)如下:軸承直徑D=14 mm，長度L=33 mm，兩排環(huán)面節(jié)流供氣小孔，每排8個，小孔直徑d=0.4 mm，小孔至端面距離L1=9.5 mm，軸承平均半徑間隙Cr=0.022 5 mm，偏心率ε=0，氣體工質(zhì)為氦氣.

圖2 ε=0時理想模型Fig.2 Ideal gas bearing model with ε =0

理想情況下，當(dāng)滿足收斂條件時，氣膜內(nèi)壓力分布應(yīng)該左右上下完全對稱，承載能力=0.但實際情況會有以下原因造成載荷不為0，即存在一個小的載荷誤差Wσ，正常情況下，載荷誤差Wσ相比有偏心情況下的載荷為小值，不會對結(jié)果產(chǎn)生本質(zhì)影響.

1)邊界條件2)假設(shè)被剪開兩條邊上的壓力其與相鄰的兩條邊上壓力相等，當(dāng)網(wǎng)格數(shù)不足夠多時，相鄰兩條邊總會有一定壓力差，會導(dǎo)致壓力分布并不完全對稱.

2)本程序利用辛樸森公式求壓力對面積的積分得到承載能力也會一定誤差.

3.1 前端角γ對結(jié)果影響

在求軸頸軸承靜態(tài)特性時，將圓柱形氣膜展開成如圖1長方形計算區(qū)域過程中，從什么位置剪開，即氣膜圓柱面剪開線與其垂直軸線組成的平面和按豎直軸線剪開平面的夾角γ(前端角)，對計算結(jié)果是否會產(chǎn)生影響進(jìn)行探討.本文按照 10°、22.5°和30°方向分別展開，并且都以兩節(jié)流小孔的中心線將計算域劃分為8個區(qū)域.圖3為3種情況下的計算結(jié)果，從圖中可以看出，載荷誤差Wσ都在一個相對比較小的范圍，并且沒有遞增或遞減趨勢，3種情況下得到的氣體消耗量Q基本上都相等，因此說明前端角γ對求軸頸軸承靜態(tài)性能并沒有影響，即可以從任意位置將圓柱軸承展開成平面計算區(qū)域.

此規(guī)律為劃分網(wǎng)格區(qū)域帶來了便利，以8個供氣孔為例，選擇22.5°位置將氣膜展開成長方形計算區(qū)域，可以得到大小相等的8個區(qū)域，使得在節(jié)流小孔位置布網(wǎng)格點更方便，另外還可以避免被剪開的兩條邊上碰到有節(jié)流小孔，使程序處理變得復(fù)雜.

圖3 前端角γ對計算結(jié)果影響Fig.3 Load error and gas consumption versus pressure plot

3.2 網(wǎng)格數(shù)與收斂精度δ對結(jié)果影響

圖4 網(wǎng)格數(shù)和收斂精度δ對計算結(jié)果影響Fig.4 Load error and gas consumption versus pressure plot

一般數(shù)值計算問題中，網(wǎng)格數(shù)與計算精度成正比，即網(wǎng)格數(shù)越多計算結(jié)果精度越高.但本文發(fā)現(xiàn)此規(guī)律在超松弛迭代法(SOR)中并不適用，網(wǎng)格數(shù)與計算結(jié)果精度并不是成正比，而是與收斂精度δ有關(guān)，當(dāng)收斂精度8選擇不當(dāng)，網(wǎng)格數(shù)越多計算結(jié)果越不精確，有的甚至得到錯誤的結(jié)果，本文對網(wǎng)格數(shù)與收斂精度δ對計算結(jié)果的影響進(jìn)行了研究.

圖4表明，收斂精度取δ=10-4時，當(dāng)長寬的網(wǎng)格數(shù)m=n=128和256時，雖然程序也滿足了收斂條件，但載荷誤差Wσ很大，特別是當(dāng)m=n=256時，載荷誤差δ能達(dá)到約20 N，很顯然計算結(jié)果是不可信的.相比而言，網(wǎng)格數(shù)m=n=32和64時，載荷誤差Wσ在一個比較小的范圍，得到的結(jié)果更接近實際.當(dāng)網(wǎng)格數(shù)m=n=128不變，提高收斂精度到δ=10-5，載荷誤差Wδ在一個比較小的范圍內(nèi)變化，氣體消耗量與網(wǎng)格數(shù)m=n=32和64時很接近，繼續(xù)提高收斂精度到δ=5×10-6，載荷誤差Wσ并沒有繼續(xù)減小，而也是在相同小的范圍內(nèi)變化.

因此，在用有限差分解Reynolds方程過程中，網(wǎng)格數(shù)需要與收斂精度δ相匹配，以本文為例，本程序采用具有二階精度的中心差分格式離散Reynolds方程，數(shù)值解精度可以達(dá)到O(Δx2)，步長Δx由網(wǎng)格點數(shù)(m，n)和軸承特征長度(或?qū)挾?決定，收斂精度δ值至少要小于Δx2數(shù)量級才能保證超松弛迭代法(SOR)也具有二階精度，本文中如果是m=n=64，δ值取1 ×10-4合適，如果 m=n=128，δ值至少需取5×10-5合適，依此類推.網(wǎng)格數(shù)增多需要同時提高收斂精度δ，網(wǎng)格數(shù)一定，當(dāng)收斂精度δ值小于某臨界值以后，再小對提高數(shù)值解精度并沒有太大幫助，并且還增加了計算成本，因此一般情況下δ值取到1×10-6即可.

圖5 程序流程圖Fig.5 Flow chart for calculation program

3.3 偏位角α影響

偏位角α由負(fù)載和轉(zhuǎn)速唯一確定，當(dāng)轉(zhuǎn)子的轉(zhuǎn)速不為零時，軸頸的偏位角也不為零［10］，此時的氣膜厚度表達(dá)式應(yīng)該按實際偏位角方向來計算，而大多數(shù)求解過程都是忽略了偏位角α對氣膜厚度的影響.

本程序先假設(shè)一個小的偏位角度，在偏位角方向上將氣膜展開，此時的氣膜厚度表達(dá)式即為實際的氣膜厚度，當(dāng)滿足收斂條件以后，求得實際的偏位角α'=arctan(Wt/Wn)，判斷α'和α是否相等，不相等則改變初設(shè)的α重新返回計算，直至α'和α的差值與α的比小于ζ，此時求得的α'即為實際情況下的偏位角.

為了節(jié)約篇幅，本文以2個節(jié)流孔來說明其迭代的流程圖如圖5，節(jié)流孔數(shù)增多就是再添加如區(qū)域I，Π的質(zhì)量守恒循環(huán)，流程圖中虛線部分表示有轉(zhuǎn)速時需要加入的判斷偏位角是否合適的嵌套循環(huán)(如圖中A位置)，在求軸頸軸承靜態(tài)性能時按實線流程直接計算軸承流量.

4 算例

文獻(xiàn)［11］利用Marc軟件中的軸承分析模塊求解了如下參數(shù)的氣體軸承靜態(tài)性能，并且和試驗進(jìn)行了對比，具體參數(shù)如下:軸承直徑D=45 mm，長度L=35 mm，共有兩排環(huán)面供氣孔，每排8個，節(jié)流孔直徑d=0.3 mm，節(jié)流孔至斷面距離l=10 mm，軸承平均間隙Cr=0.02 mm，供氣壓力ps=6.5×105Pa.

氣體消耗量和承載能力隨偏心率變化的對比關(guān)系分別如圖6、7.

與文獻(xiàn)中計算結(jié)果相比，誤差在5%以下，產(chǎn)生誤差的原因除了程序本身以外，還與各種方法對文獻(xiàn)中未提供的其它參數(shù)取值差異有關(guān)，另外收斂精度和網(wǎng)格數(shù)的選擇對計算結(jié)果也有一定影響.

圖6 氣體消耗量計算結(jié)果比較Fig.6 Comparison of the results for gas consumption

圖7 承載能力計算結(jié)果比較Fig.7 Comparison of the results for load capacity

5 結(jié)論

1)本文利用有限差分法編寫了一套可以求靜壓軸頸軸承動靜態(tài)性能的程序，引入了零偏心情況下的理想模型，討論了前端角γ對求軸頸軸承靜態(tài)性能計算結(jié)果的影響，結(jié)果證明前端角γ對求軸頸軸承靜態(tài)性能計算結(jié)果沒有影響;

2)對網(wǎng)格數(shù)與收斂精度δ進(jìn)行研究發(fā)現(xiàn)，網(wǎng)格數(shù)需要與收斂精度δ相匹配，否則對計算結(jié)果會造成很大影響;

3)引入了偏位角α的嵌套循環(huán)，考慮了實際氣膜厚度表達(dá)式對軸承性能的影響，計算結(jié)果與有實驗驗證的算例比較，兩者誤差在5%以下，該方法已成功用于對俄羅斯氦透平膨脹機(jī)靜壓氣體軸承端的改造上，效果良好.

［1］吳華鵬，孔憲梅.GMRES算法在雷諾方程數(shù)值解法中的應(yīng)用［J］.潤滑與密封，2000(1):2-4.

WU Huapeng，KONG Xianmei.The application of GMRES algorithm in calculation Reynolds equation［J］.Lubrication Engineering，2000(1):2-4.

［2］彭萬歡，蒲如平.氣體靜壓軸頸軸承的靜態(tài)特性分析及其試驗研究［J］.潤滑與密封，2006(12):204-206.

PENG Wanhuan，PU Ruping.A static performance analysis and experimental study of aerostatic journal bearing［J］.Lubrication Engineering，2006(12):204-206.

［3］LO Chengying，WANG Chengchi，LEE Yuhan.Performance analysis of high-speed spindle aerostatic bearings［J］.Tribology International，2005，38:5-14.

［4］張鳴，朱煜，段廣洪.基于FEMLAB的氣浮軸承靜態(tài)性能求解方法［J］.潤滑與密封，2006(2):63-65.

ZHANG Ming ，ZHU Yu ，DUAN Guanghong.A new method based FEMLAB to predictstatic performance ofair bearing［J］.Lubrication Engineering，2006(2):63-65.

［5］十合晉一.氣體軸承設(shè)計、制作與應(yīng)用［M］.韓渙臣，譯.北京:宇航出版社，1988:12-20.

［6］POWELL J W.Design of aerostatic bearings［M］.London:Machinery Publishing Co.Ltd，1970:43-47.

［7］王云飛.氣體潤滑理論與氣體軸承設(shè)計［M］.北京:機(jī)械工業(yè)出版社，1999:19-30.

［8］劉暾，劉育華，陳世杰.靜壓氣體潤滑［M］.哈爾濱:哈爾濱工業(yè)大學(xué)出版社，1990:40-60.

［9］黃云清.數(shù)值計算方法［M］.北京:科學(xué)出版社，2009:70-85.

［10］池長青.氣體動靜壓軸承的動力學(xué)和熱力學(xué)［M］.北京:北京航空航天大學(xué)，2008:60-74.

［11］張恩龍.高速全支撐氣體靜壓電主軸的承載特性研究［D］.廣州:廣州工業(yè)大學(xué)，2005:65-86.