基于MATLAB與CFX的動壓氣體軸承計算對比

吳炎，謝振宇，陳李成，郝建勝，李超

(南京航空航天大學 直升機傳動技術重點實驗室，江蘇 南京 210016)

0 引言

人字槽徑向動壓氣體軸承潤滑介質是氣體，具有阻尼小、功耗低、污染少、壽命長等優點，適用于相對極端的環境，應用領域較廣[1-2]。但該軸承的承載力偏低，剛度偏弱[3-4]，理論分析難度較大，因此相關理論計算數據和試驗數據比較少。

1886年REYNOLDS O推導出了不可壓縮流體的Reynolds方程，描述了相對運動表面的運動速度、介質黏度、間隙形狀與介質膜壓力分布的關系，為流體潤滑技術奠定了理論基礎[5]。1962年-1965年間，CASTELLI V、REDDY M M等在利用有限差分法和有限元法等數值方法在求解靜態雷諾方程上作出了杰出貢獻[6-7]。有限差分法和有限元法廣泛應用于求解動壓軸承的靜態性能。

計算機技術的高速發展為Reynolds方程的求解提供了極大的助力。在MATLAB上對Reynolds方程進行編程，利用有限差分法求解動壓氣體軸承的壓力分布逐漸成為一種主流形式，也是CFD計算最早采用的方法。而一些專業的CFD軟件逐漸崛起，其中ANSYS軟件的CFX和Fluent模塊使用最廣泛。

二者計算方式存在較大差異，這意味著計算結果也可能存在較大差異，因此分析不同計算方法的差異是非常必要的。

本文以人字槽徑向動壓氣體軸承為研究對象，對比分析ANSYS軟件中CFX模塊與基于MATLAB的有限差分法計算結果的差別，為使用兩種方法分析動壓氣體軸承提供參考。

為方便敘述，在后文中，將人字槽徑向動壓氣體軸承簡稱為氣體軸承。

1 氣體軸承模型

取氣體軸承的直徑D為55×10-3m，氣體軸承寬度L為24×10-3m，平均氣膜厚度cr為3×10-5m，槽數Nc為16，槽寬bc為4.32×10-3m，槽深hc為5×10-5m，單邊槽長lc為6×10-3m，槽角為19°，額定轉速為30000 r/min，偏心率為0.6。氣體軸承模型如圖1所示。

圖1 氣體軸承外壁模型

2 靜態特性的MATLAB求解

2.1 雷諾方程建立與離散化

從N-S方程出發進行推導，得到柱坐標系下可壓縮氣體的廣義曲面支承靜態二元流動雷諾方程：

(1)

式中：R為氣體軸承半徑，m；θ為圓周方向弧度數，rad；p為氣膜壓力，Pa；h為氣膜厚度，m；μ為氣體黏度，N·s/m2；z為氣體軸承軸向坐標，m；ω為工作角速度，rad/s。將式(1)無量綱化可得：

(2)

式中：Z=z/L；H=h/cr；P=p/p0；λ=L/D；Λ為氣體的可壓縮系數；L為氣體軸承軸向寬度，m；cr為氣體軸承平均氣膜厚度，m；p0為環境壓力，Pa；μ0為25 ℃時空氣的動力黏度；ν為擾動角速度，rad/s；D為氣體軸承直徑，m。利用中心差分法對式(2)進行離散，化簡后提取系數矩陣得

ai,jδi-1,j+bi,jδi+1,j+ci,jδi,j+di,jδi,j-1+ei,jδi,j+1=-Si,j

(3)

式中：ai,j、bi,j、ci,j、di,j、ei,j、Si,j為與氣膜厚度、氣膜壓力有關的系數矩陣；δi,j=(Pn+1-Pn)i，j表示在第(i，j)節點處，第n+1次迭代與第n次迭代壓力差。

迭代計算的收斂條件有很多種，其中最常用的是相對精度判斷(相對收斂準則)，相應的收斂公式為

(4)

2.2 氣膜厚度方程建立

如圖2所示，設兩中心連線O1O2與載荷W方向為偏位角φ，偏心距e與平均氣膜厚度cr的比值為偏心率ε。忽略無窮小量，整理后可得無量綱氣膜厚度為

(5)

其中無槽處槽深hc=0。

圖2 轉子相對位置示意圖

將氣膜沿圓周方向展開，其分布模型如圖3所示。

圖3 氣膜軸向展開示意圖

對氣膜模型進行網格劃分，并對網格節點進行分類，用于判斷網格節點是否處于槽區。氣體軸承網格劃分如圖4所示。

圖4 氣體軸承網格劃分示意圖

圖4中，lc為單邊槽長，bc為槽寬，α為槽角。一組槽可分為兩個槽區，分別由四條直線圍成。為方便計算，使某一組槽的一個邊界點置于原點處，并使整個槽區處于第一象限，設k=tanα，則有：

節點處于左半部分槽區條件

(6)

節點處于右半部分槽區條件

(7)

其中：Lc=lc/L為無量綱槽長；B=bc/R為無量綱槽寬。當節點滿足式(6)或式(7)的條件時，可認為該節點落在槽區中，由此可以判斷節點所在位置。

采用MATLAB進行數值計算的編程，動壓氣體軸承的性能求解程序框圖如圖5所示。

圖5 MATLAB求解流程圖

3 靜態特性的CFX求解

3.1 氣膜模型建立與網格劃分

Reynolds方程是N-S方程的簡化版，推導基礎為流體的動量守恒。而CFX的計算則同時兼顧了動量守恒和質量守恒，融合了有限元法和有限體積法，內置了非常豐富的物理模型，很適合用于氣體軸承這種氣膜幾何尺寸跨度極大的模型計算。

利用UG建立氣體軸承外壁與轉子部分的三維模型，導入Workbench的DM模塊中，對氣體軸承的內流域進行抽取，得到氣膜模型。氣體軸承的氣膜模型如圖6所示。

圖6 氣膜三維幾何模型

由于氣膜厚度與氣體軸承直徑相差兩個數量級，而氣體軸承存在偏心和人字槽，模型不能應用周期陣列，且在模型傳輸過程中易出現幾何畸變，生成網格時易產生負網格。ANSYS內置的ICEM模塊內置強大的幾何修補功能，且能對模型進行結構化網格的劃分，保證網格質量，提高計算精度。

針對氣膜厚度與氣體軸承直徑相差過大的問題，可利用ICEM的塊映射，通過雕塑法在拓撲空間進行模型的塊劃分，利用擠出得到與槽區相合的塊。由于模型存在偏心，需要使每一個塊與實體模型間建立點、線、面的關聯，以保證每一個塊與物理模型形成完美映射，生成光滑貼體的高質量結構網格。而結構網格在生成速度、計算精度、抗畸變度等方面，往往更加優越。該模型的塊劃分結果如圖7所示，生成的結構化網格(局部)如圖8所示。

圖7 氣膜模型的塊劃分

圖8 氣膜模型的局部結構化網格

3.2 邊界條件及求解設置

將網格導入CFX求解器，將氣體軸承的潤滑介質定為25℃下滿足理想氣體條件的空氣。由于氣體軸承的雷諾系數很小，模型的domain選為層流模型。氣體軸承氣膜兩側界面同時充當入口和出口，因此該界面為開放式邊界，采用Entrainment作為邊界氣體交換條件。轉子壁面設為旋轉面，氣體軸承外壁設為固定面。為保證計算精度，采用高階求解對氣體軸承性能進行計算。

CFX中對殘差的定義為：

rn=c-Aθn

(8)

其中：rn為第n次迭代計算后所得的殘差；c為CFX求解方程中質量守恒微分方程與動量守恒微分方程離散、線性化后的常數項；A為系數矩陣；θn為待求量(壓力、速度)第n次迭代計算值。

為保證計算精度，取殘差<10-6時停止計算。

4 計算結果分析

4.1 壓力分布

氣體軸承轉速為30000 r/min、偏心率為0.6時，MATLAB計算所得壓力分布如圖9所示，CFX計算所得壓力分布如圖10所示。

圖9 MATLAB氣膜壓力分布

圖10 CFX氣膜壓力分布圖

可以看出，通過MATLAB進行數值計算所得結果與CFX計算結果在分布情況上基本一致。在氣體軸承兩邊以及氣膜較厚處，壓力大小與外界壓強基本一致，在氣膜較薄處壓力相對較大，且槽根處壓力相對較大，說明氣體由于黏性在氣膜內被不斷壓縮，聚積在氣膜較薄處和人字槽槽根處，形成動壓效應和泵吸效應，產生支承力支承轉子。

4.2 靜態承載力結果對比

轉速為30000 r/min時承載力隨偏心率變化情況如圖11所示，偏心率為0.6時承載力隨轉速變化情況如圖12所示。

圖11 承載力隨偏心率變化關系

圖12 承載力隨轉速變化關系

可以看出，偏心率較小時，二者計算結果基本一致，CFX計算結果略高于MATLAB；當偏心率>0.6時，MATLAB計算結果更大，在極限偏心率下二者差別進一步拉大。而隨著轉速的升高，二者計算結果幾乎沒有明顯區別。

5 結語

通過不同偏心率和轉速下MATLAB數值計算與CFX仿真求解所得結果的對比，可以得到以下結論：

1)動壓氣體軸承的承載力分別隨偏心率和轉速的增大而增大，MATLAB數值計算與CFX仿真分析所得結果趨勢變化基本一致，非極限偏心率下二者所得結果基本一致。

2)在極限偏心率下MATLAB數值計算結果與CFX仿真結果存在較大偏差。

本文研究的算例雖然有限，但所得結論仍可為動壓氣體軸承的計算提供有益的參考。