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

(南京航空航天大學(xué) 直升機(jī)傳動(dòng)技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210016)

0 引言

人字槽徑向動(dòng)壓氣體軸承潤(rùn)滑介質(zhì)是氣體，具有阻尼小、功耗低、污染少、壽命長(zhǎng)等優(yōu)點(diǎn)，適用于相對(duì)極端的環(huán)境，應(yīng)用領(lǐng)域較廣[1-2]。但該軸承的承載力偏低，剛度偏弱[3-4]，理論分析難度較大，因此相關(guān)理論計(jì)算數(shù)據(jù)和試驗(yàn)數(shù)據(jù)比較少。

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

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

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

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

為方便敘述，在后文中，將人字槽徑向動(dòng)壓氣體軸承簡(jiǎn)稱為氣體軸承。

1 氣體軸承模型

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

圖1 氣體軸承外壁模型

2 靜態(tài)特性的MATLAB求解

2.1 雷諾方程建立與離散化

從N-S方程出發(fā)進(jìn)行推導(dǎo)，得到柱坐標(biāo)系下可壓縮氣體的廣義曲面支承靜態(tài)二元流動(dòng)雷諾方程：

(1)

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

(2)

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

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為與氣膜厚度、氣膜壓力有關(guān)的系數(shù)矩陣；δi,j=(Pn+1-Pn)i，j表示在第(i，j)節(jié)點(diǎn)處，第n+1次迭代與第n次迭代壓力差。

迭代計(jì)算的收斂條件有很多種，其中最常用的是相對(duì)精度判斷(相對(duì)收斂準(zhǔn)則)，相應(yīng)的收斂公式為

(4)

2.2 氣膜厚度方程建立

如圖2所示，設(shè)兩中心連線O1O2與載荷W方向?yàn)槠唤铅?，偏心距e與平均氣膜厚度cr的比值為偏心率ε。忽略無窮小量，整理后可得無量綱氣膜厚度為

(5)

其中無槽處槽深hc=0。

圖2 轉(zhuǎn)子相對(duì)位置示意圖

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

圖3 氣膜軸向展開示意圖

對(duì)氣膜模型進(jìn)行網(wǎng)格劃分，并對(duì)網(wǎng)格節(jié)點(diǎn)進(jìn)行分類，用于判斷網(wǎng)格節(jié)點(diǎn)是否處于槽區(qū)。氣體軸承網(wǎng)格劃分如圖4所示。

圖4 氣體軸承網(wǎng)格劃分示意圖

圖4中，lc為單邊槽長(zhǎng)，bc為槽寬，α為槽角。一組槽可分為兩個(gè)槽區(qū)，分別由四條直線圍成。為方便計(jì)算，使某一組槽的一個(gè)邊界點(diǎn)置于原點(diǎn)處，并使整個(gè)槽區(qū)處于第一象限，設(shè)k=tanα，則有：

節(jié)點(diǎn)處于左半部分槽區(qū)條件

(6)

節(jié)點(diǎn)處于右半部分槽區(qū)條件

(7)

其中：Lc=lc/L為無量綱槽長(zhǎng)；B=bc/R為無量綱槽寬。當(dāng)節(jié)點(diǎn)滿足式(6)或式(7)的條件時(shí)，可認(rèn)為該節(jié)點(diǎn)落在槽區(qū)中，由此可以判斷節(jié)點(diǎn)所在位置。

采用MATLAB進(jìn)行數(shù)值計(jì)算的編程，動(dòng)壓氣體軸承的性能求解程序框圖如圖5所示。

圖5 MATLAB求解流程圖

3 靜態(tài)特性的CFX求解

3.1 氣膜模型建立與網(wǎng)格劃分

Reynolds方程是N-S方程的簡(jiǎn)化版，推導(dǎo)基礎(chǔ)為流體的動(dòng)量守恒。而CFX的計(jì)算則同時(shí)兼顧了動(dòng)量守恒和質(zhì)量守恒，融合了有限元法和有限體積法，內(nèi)置了非常豐富的物理模型，很適合用于氣體軸承這種氣膜幾何尺寸跨度極大的模型計(jì)算。

利用UG建立氣體軸承外壁與轉(zhuǎn)子部分的三維模型，導(dǎo)入Workbench的DM模塊中，對(duì)氣體軸承的內(nèi)流域進(jìn)行抽取，得到氣膜模型。氣體軸承的氣膜模型如圖6所示。

圖6 氣膜三維幾何模型

由于氣膜厚度與氣體軸承直徑相差兩個(gè)數(shù)量級(jí)，而氣體軸承存在偏心和人字槽，模型不能應(yīng)用周期陣列，且在模型傳輸過程中易出現(xiàn)幾何畸變，生成網(wǎng)格時(shí)易產(chǎn)生負(fù)網(wǎng)格。ANSYS內(nèi)置的ICEM模塊內(nèi)置強(qiáng)大的幾何修補(bǔ)功能，且能對(duì)模型進(jìn)行結(jié)構(gòu)化網(wǎng)格的劃分，保證網(wǎng)格質(zhì)量，提高計(jì)算精度。

針對(duì)氣膜厚度與氣體軸承直徑相差過大的問題，可利用ICEM的塊映射，通過雕塑法在拓?fù)淇臻g進(jìn)行模型的塊劃分，利用擠出得到與槽區(qū)相合的塊。由于模型存在偏心，需要使每一個(gè)塊與實(shí)體模型間建立點(diǎn)、線、面的關(guān)聯(lián)，以保證每一個(gè)塊與物理模型形成完美映射，生成光滑貼體的高質(zhì)量結(jié)構(gòu)網(wǎng)格。而結(jié)構(gòu)網(wǎng)格在生成速度、計(jì)算精度、抗畸變度等方面，往往更加優(yōu)越。該模型的塊劃分結(jié)果如圖7所示，生成的結(jié)構(gòu)化網(wǎng)格(局部)如圖8所示。

圖7 氣膜模型的塊劃分

圖8 氣膜模型的局部結(jié)構(gòu)化網(wǎng)格

3.2 邊界條件及求解設(shè)置

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

CFX中對(duì)殘差的定義為：

rn=c-Aθn

(8)

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

為保證計(jì)算精度，取殘差<10-6時(shí)停止計(jì)算。

4 計(jì)算結(jié)果分析

4.1 壓力分布

氣體軸承轉(zhuǎn)速為30000 r/min、偏心率為0.6時(shí)，MATLAB計(jì)算所得壓力分布如圖9所示，CFX計(jì)算所得壓力分布如圖10所示。

圖9 MATLAB氣膜壓力分布

圖10 CFX氣膜壓力分布圖

可以看出，通過MATLAB進(jìn)行數(shù)值計(jì)算所得結(jié)果與CFX計(jì)算結(jié)果在分布情況上基本一致。在氣體軸承兩邊以及氣膜較厚處，壓力大小與外界壓強(qiáng)基本一致，在氣膜較薄處壓力相對(duì)較大，且槽根處壓力相對(duì)較大，說明氣體由于黏性在氣膜內(nèi)被不斷壓縮，聚積在氣膜較薄處和人字槽槽根處，形成動(dòng)壓效應(yīng)和泵吸效應(yīng)，產(chǎn)生支承力支承轉(zhuǎn)子。

4.2 靜態(tài)承載力結(jié)果對(duì)比

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

圖11 承載力隨偏心率變化關(guān)系

圖12 承載力隨轉(zhuǎn)速變化關(guān)系

可以看出，偏心率較小時(shí)，二者計(jì)算結(jié)果基本一致，CFX計(jì)算結(jié)果略高于MATLAB；當(dāng)偏心率>0.6時(shí)，MATLAB計(jì)算結(jié)果更大，在極限偏心率下二者差別進(jìn)一步拉大。而隨著轉(zhuǎn)速的升高，二者計(jì)算結(jié)果幾乎沒有明顯區(qū)別。

5 結(jié)語

通過不同偏心率和轉(zhuǎn)速下MATLAB數(shù)值計(jì)算與CFX仿真求解所得結(jié)果的對(duì)比，可以得到以下結(jié)論：

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

2)在極限偏心率下MATLAB數(shù)值計(jì)算結(jié)果與CFX仿真結(jié)果存在較大偏差。

本文研究的算例雖然有限，但所得結(jié)論仍可為動(dòng)壓氣體軸承的計(jì)算提供有益的參考。