基于流體域仿真的高速風機氣浮箔片軸承動靜特性研究

張晨，許亮，孫丹，李雨薇，高飛，荊建平

(1.上海交通大學機械系統與振動國家重點實驗室，上海 200240；2.中國船舶集團有限公司第七〇三研究所， 黑龍江哈爾濱 150078；3.新鄉航空工業(集團)有限公司，河南新鄉 453000)

高速風機作為航空、艦船、化工等工業領域常見的流體機械，其制造產業的提升對工業流體機械的發展尤為重要。對于風機而言，隨著工作時間的增加，其動平衡會被逐漸破壞，軸承部位的振動會不斷加大。一旦軸承的振動大于某一限值，就會對風機造成巨大的影響，必須停機修理[1]。隨著科學技術水平的不斷提升，風機應用的環境多變復雜，風機向著體積小、質量輕的方向不斷發展，這對風機轉子系統所使用的軸承的精度提出了更高的要求。氣浮軸承憑借著精度高、發熱小、噪聲低等各項優點，正在逐步替代傳統的機械軸承，成為風機等高精度設備中的關鍵部件[2]。如今，氣浮軸承應用于風機的轉子系統中，針對風機類產品轉子系統特性的分析、臨界轉速的預測等研究便是所需要解決的問題。目前市面上常見的商業計算軟件均缺失較為系統的氣浮軸承特性計算方法，氣浮軸承的研究也在不斷地發展與完善。

工業中的高速風機多數選擇使用高速氣浮軸承，結構如圖1所示，該類軸承的特性計算難度更大[3]。并且，為了提高氣浮軸承的抗沖擊能力與自適應能力，在軸與軸承之間加入了具有柔性的金屬箔片(見圖2)，其能夠根據轉速與載荷的變化建立起合適的工作氣膜厚度，同時箔片軸承的結構阻尼效應能夠更好地抑制轉子的振動，從而有更好的穩定性。因此對多葉箔片軸承的特性進行理論分析有更加重要的現實意義[4]。

圖2 氣浮軸承箔片Fig.2 Foil for air bearing

WALOWIT等[5-6]是較早研究波箔型箔片軸承的學者，他們將起支撐作用的箔片從軸承中分離，分析單個箔片的受力情況，建立對應的軸承模型，求解得到箔片軸承的靜態特性。以HESHMAT為代表的一批學者對波箔型箔片軸承進行了比較深入的理論和實驗研究，采用有限差分法來計算彈性流體動力潤滑(EHD)問題，形成研究波箔型箔片軸承剛度和阻尼特性的經典理論方法[7-9]。

在國內，龔煥孫[10]最早開展波箔型箔片軸承剛度和阻尼特性的研究，他基于Krichhoff假設，施加載荷與邊界條件進行計算，得到軸承靜態特性，同時搭建實驗臺，對其靜態特性進行試驗研究，與理論分析的結果相符。宋來運[11]從理論和實驗驗證2個層面對動靜壓氣浮軸承轉子系統的靜態性能、動態性能和動態穩定性等綜合性能進行了系統性的研究，為設計研發高速高精度動靜壓氣浮主軸系統提供理論基礎和技術支持。劉佳琦等[12]從彈性箔片軸承的結構特性、涂層特性和建模方法3個方面對氣體彈性箔片軸承的研究進展進行綜述，展望彈性箔片軸承在高精尖等領域的前景，并對彈性箔片軸承的研究方向提出了建議。劉恒等人[13]對箔片含有預變形的多葉箔片氣體動壓軸承進行研究，將箔片簡化為彎曲梁，推導出多葉箔片的變形方程和考慮預變形后的氣膜厚度方程，計算研究轉速等參數對軸承承載力的影響。

然而，當前的氣浮軸承的特性計算多是基于一定的假設條件，通過簡化軸承特性方程來得到軸承特性。郭軍剛等[14]建立多葉油潤滑箔片軸承Reynolds方程、對應的邊界條件以及收斂條件，通過傳遞矩陣計算箔片軸承軸系臨界轉速，同時，建立對應的試驗系統進行評估與驗證。

由于眾多假設條件的存在，計算結果往往與實際結果存在偏差。因此，為了提高高速風機轉子動力學分析精度，減少理論計算與實際情況之間的偏差，本文作者針對現有的一類多箔片氣浮軸承，建立流體域模型，通過使用商業仿真軟件分析風機產品內高速氣浮箔片軸承的靜態特性與動態特性隨著轉速提升的變化規律，從而更加精準地預估風機模態、臨界轉速以及失穩轉速。在設計使用氣浮軸承的轉子系統時，依據氣浮軸承自身的動態特性變化趨勢來設計對應的轉子系統結構，從而能夠更大程度地提高轉子系統的穩定性與安全性，這是開展軸承作為支承部件的轉子系統動力學特性研究的基礎。

1 計算方法與流體域模型

針對使用氣浮軸承的風機轉子系統，氣浮軸承中的多葉箔片在轉子系統啟動時，起到支撐作用；隨著軸頸轉速的不斷提升，帶動氣浮軸承的流體域內的流體高速流動，由收斂楔進入發散楔，產生氣膜；氣膜可視作具有較大剛度的支撐體，具有足夠的承載能力，隨著氣膜的動壓力不斷提升，將軸頸慢慢托起，與箔片中間形成空隙，不產生接觸。當軸頸的載荷與氣膜所能夠提供的承載能力實現相互平衡時，即為氣浮軸承的工作狀態。

文中提出的軸承特性求解的方法原理如圖3所示，主要過程為：首先建立氣浮軸承的流體域模型，針對流體域模型進行網格劃分，通過商用仿真軟件模擬一定轉速下氣浮軸承流體域的分布情況；移動軸頸使得軸頸受到流體域的作用力與所受到的載荷相互平衡，從而得到氣浮軸承在一定載荷下的軸心軌跡圖，即為氣浮軸承的靜特性計算；在靜特性的計算結果的基礎上，賦予軸頸與氣浮軸承轉速相同的同頻擾動，讀取在擾動的作用下流體域對軸頸的作用力，并沿軸頸表面進行積分，得到軸頸沿X、Y軸方向的受力變化，從而計算氣浮軸承的剛度與阻尼系數，即為氣浮軸承的動特性計算。

圖3 軸承特性求解流程Fig.3 Flow for solving bearing characteristics

建立氣浮軸承的流體域模型時，為了更加便于仿真計算與后續的模型修正，基于流體域沿軸向分布變化較小的情況，建立流體域二維平面模型。使用該模型進行模擬仿真，得到流體域結果后，通過使用修正系數，將二維模型計算結果進行修正，得到三維流體域模型的結果。

文中研究的氣浮軸承流體域內的流體為可壓縮氣體，使用商業仿真軟件根據氣體的連續方程與動量方程等進行求解，通過仿真模擬計算得到氣浮軸承流體域內的氣體分布情況以及軸頸表面的壓力分布，列出軸頸的受力方程，從而可以求解得到8個動特性系數[11]。

在轉子軸心附近產生一個小擾動Δx、Δy，流體域對轉子的作用力ΔFx、ΔFy為

(1)

將位移與受力取諧波分量可得：

(2)

(3)

將諧波分量代入式(1)，可得：

(4)

求解該方程組，即可得到氣浮軸承的剛度系數與阻尼系數。

2 模型網格劃分與邊界條件設定

2.1 模型網格劃分

由于氣浮軸承結構特殊性，為滿足CFD仿真模擬計算得到的流體域結果的精確度同時提高計算過程效率，選擇使用流體域二維模型進行仿真計算[11]，同時需要對流體域進行高質量的網格劃分，得到更合適的網格結構。網格分為結構化網格與非結構化網格，因為該類氣浮軸承的流體域模型比較簡單，且比較規則，沒有太多復雜的流體區域，因此，選擇劃分結構化網格。該類網格生成速度比較快，并且數據結構簡單，可以滿足氣浮軸承流體域的仿真計算。網格單元越小的情況下，得到的計算結果精度更高，但同時會增加仿真模擬的計算時間；且在仿真模擬計算的過程中，需要控制軸頸移動，導致流體域隨之改變，涉及動網格設置問題，對網格的質量要求也隨之提高。因此，需要選取滿足計算精度的網格尺寸與網格數量。文中采用的結構化網格劃分流體域模型(見圖4)，經過網格無關性驗證，采用網格數量分別為51 090、239 820、420 000與802 670的流體域模型。在該網格數下分別移動軸頸至工作轉速下的位置，得到的誤差如圖5所示?？梢钥闯?，節點數為399 800，網格數量為420 000的網格劃分得到的模型即可滿足仿真計算求解的精度要求。質量較高的流體域模型網格為后續的氣浮軸承靜特性與動特性的計算求解奠定了基礎。

圖4 氣浮軸承流體域模型及網格劃分示意Fig.4 Fluid domain model and grid generation of air bearing

圖5 網格無關性驗證Fig.5 Grid independence verification

2.2 模型邊界條件設定

使用該類氣浮軸承的風機轉子系統的工作轉速為100 000 r/min，軸頸所受到的載荷為轉子自身的重量3 N，流體域內的材料為可壓縮氣體，氣體的密度為8.13 kg/m3，氣體的壓縮系數為1.4，在仿真模擬的計算過程中不考慮氣體黏性的影響，氣體的壓力為800 kPa。模擬仿真的求解選擇Pressure-Based，Velocity Formulation選擇Absolute，Time選擇Transient，2D Space選擇Planar。

模型的設置中，Energy選擇On，Viscous選擇SSTk-ε，該模型可以較好地處理求解邊界層網格數據，在近壁面的位置能夠更好地進行處理。因為氣浮軸承的流體域模型的特點，可以看出，k-ε模型能夠更高精度地處理求解流體域內的流動情況；同時，選用SST，相較于BSL，依然保留了在近壁面采用k-ε模型，但SST簡化了模型，減小計算量，對壁面的距離更加敏感，能夠進行更加精確地仿真求解；同時，壁面的設置中加入軸頸的轉速與流體域的壓力邊界條件，以及考慮移動軸頸的動網格，根據軸頸位置的改變來調整網格結構。最后，根據網格的尺寸與軸頸的移動設定合適的時間步長與時間步數來得到更加精確的結果。得到靜特性的計算結果之后，賦予軸頸同頻擾動，求解得到軸頸表面的壓力分布，對軸頸表面進行積分得到軸頸沿X、Y軸方向受到的作用力，從而求解得到氣浮軸承的剛度系數與阻尼系數。完成100 000 r/min工作轉速下的氣浮軸承的特性計算之后，繼續求解60 000、80 000 r/min轉速下的氣浮軸承的特性，從而得到氣浮軸承從開始啟動到轉速提升至工作轉速這一過程中軸心位置的變化以及氣浮軸承的靜特性與動特性的變化趨勢。

3 仿真結果及分析

對氣浮軸承流體域模型進行計算流體力學求解并進行仿真計算，可以得到在不同的轉速下，氣浮軸承流體域內的氣體速度分布、湍動能分布以及氣體壓力分布情況。圖6、7、8所示為100 000 r/min轉速下氣浮軸承流體域內氣體的速度、壓力分布情況以及流體域內的湍動能分布。

從圖6中可以看出，在流體域內靠近軸頸邊界的氣體速度較大，隨之速度沿徑向方向不斷下降，在箔片邊界附近的氣體速度最低，與實際情況相比較為符合。從圖7中可以看出，流體域內氣體的壓力最大值發生在流體域厚度最小附近，該現象產生的原因可由流體動壓潤滑理論解釋。流體動壓潤滑是指軸頸在旋轉時，在軸頸與軸承的摩擦表面之間會充滿流體，由于流體自身具有黏性，當軸頸達到足夠高的旋轉速度時，流體會被軸頸帶入軸頸與軸承表面之間的空間。對于文中所研究的氣浮軸承，即為將流體帶入軸頸與箔片之間的流體域內，在軸頸與箔片之間形成楔形間隙，從而形成流體動壓效應；被帶入楔形間隙的流體會產生壓力，該壓力與軸頸賦予的載荷(即為外載荷)平衡時，軸頸與箔片之間形成穩定的氣膜，從而實現流體的動壓潤滑。因此，流體域內氣體的壓力分布情況與理論分析結果相符合。如圖8所示，軸承的流體域內變化較大的湍動能主要發生在箔片搭接的位置，即為流體域內氣體厚度較大的位置，因為該處為箔片搭接，流體域邊界發生突變，邊界的變化會導致流體域厚度隨之改變，從而在此處更易發生湍流。

圖6 100 000 r/min轉速下流體域速度分布Fig.6 Velocity distribution in fluid domain at speed of 100 000 r/min

圖7 100 000 r/min轉速下流體域壓力分布Fig.7 Fluid domain pressure distribution at speed of 100 000 r/min

圖8 100 000 r/min轉速下流體域湍動能分布Fig.8 Turbulent kinetic energy distribution in fluid domain at speed of 100 000 r/min

圖9—11所示為氣浮軸承在80 000和60 000 r/min轉速下流體域內的氣體速度分布、湍流分布以及氣體壓力分布情況。

圖9 不同轉速下氣體速度分布Fig.9 Speed distribution of gas at different speeds：(a) 80 000 r/min；(b)60 000 r/min

圖10 不同轉速下氣體壓力分布Fig.10 Pressure distribution of gas at different speeds： (a) 80 000 r/min；(b)60 000 r/min

圖11 不同轉速下氣體湍動能分布Fig.11 Turbulent kinetic energy distribution at different speeds： (a) 80 000 r/min；(b)60 000 r/min

從圖9—11中可以看出不同轉速下，速度分布、壓力分布以及湍動能分布情況以及變化趨勢都大致相同，不同之處在于轉速的變化，軸心位置會隨之改變，承載相同壓力的氣膜形狀發生變化，但氣膜對軸頸的壓力合力仍然保持與軸承的載荷相平衡，從而保證氣膜能夠穩定存在，軸頸能夠正常運轉。

根據賦予軸頸同頻擾動后求解得到軸頸表面的壓力分布情況，對軸頸表面進行積分得到軸頸沿X、Y軸方向受到的作用力，圖12與圖13所示為賦予軸頸同頻擾動后，軸頸X、Y軸方向的位置變化與軸頸所受到的2個方向的壓力變化情況。

圖13 軸頸表面壓力變化Fig.13 Change of journal surface pressure

從圖12、13中可以看出，軸頸表面所受到的壓力隨著軸心位置的變化呈現周期性變化的規律。因此，需要通過數據處理，從壓力變化曲線中找到與軸心位置變化頻率相同的曲線；同時求解軸心位置變化曲線與壓力變化曲線之間的相位差，根據頻率與相位差值即可得到求解剛度系數與阻尼系數所需要的壓力變化表達式。使用相同的方法，賦予同一轉速2個不同的同頻擾動，處理數據得到對應的壓力變化表達式，將表達式代入方程組內，從而求解得到氣浮軸承的剛度系數與阻尼系數。圖14與圖15所示為剛度系數與阻尼系數隨著轉速提升的變化情況。

圖15 阻尼系數隨轉速的變化Fig.15 Variation of damping coefficient with speed

從圖14、15中可以看出，X軸與Y軸方向的直接剛度系數隨著轉速的提升不斷下降，直接阻尼系數的變化趨勢則是先下降而后上升；交叉剛度系數隨著轉速的提升先下降而后上升，X軸與Y軸的交叉阻尼系數的變化趨勢則相反。這是因為根據動壓潤滑理論與方程可知，任意一點的氣膜壓力與該點的速度梯度的導數有關。因此，對于某一固定位置，當軸頸的轉速較低時，軸頸偏心率較大；隨著轉速的提升，軸頸偏心率會逐漸減小，在此過程中，該點的速度梯度對應發生變化，從而產生的氣膜壓力也隨之發生改變，代入方程求解得到的剛度系數與阻尼系數呈現不斷下降的趨勢。隨著軸頸轉速的不斷提升，軸頸的偏心率變化會不斷變小，該位置的速度梯度會隨之增大，產生的氣膜壓力也不斷提高，從而會使得剛度系數與阻尼系數隨著轉速的不斷提升呈現前后不同的變化趨勢。

圖14 剛度系數隨轉速的變化Fig.14 Variation of stiffness coefficient with speed

4 結論

建立了高速氣浮軸承流體域模型，通過仿真計算方法求解了軸承特性，分析了氣浮軸承工作狀態下穩定時的軸心位置隨著轉速提升對應的變化情況，以及氣浮軸承的動態特性系數與轉速的關系，發現隨著轉速的不斷提升，穩定時的軸心位置會逐漸向氣浮箔片軸承的圓心靠近；同時，軸承的直接剛度系數在一定的工作轉速范圍內隨著轉速的提升而下降，直接阻尼系數則在工作轉速范圍內下降到達低值后上升。研究結果可為預估風機的失穩轉速與臨界轉速提供理論支撐。