基于Workbench的高壓圓盤氣體軸承共軛傳熱研究

郭良斌，吳永良

(1.武漢科技大學冶金裝備及其控制教育部重點實驗室，湖北武漢 430081；2.武漢科技大學機械傳動與制造工程湖北省重點實驗室，湖北武漢 430081)

隨著氣體軸承向高轉速、高供氣壓力方向發展，氣膜中熱量的產生和遷移愈加顯著，不僅導致軸承及相關部件的熱變形，而且還改變了軸承在實際工況下的工作間隙，進而使軸承的實際特性偏離理想設計狀態[1]。因此，氣體軸承的熱穩定性逐漸成為制約其工作性能正常發揮的一個關鍵因素[2-6]。郭良斌和徐凡[7]從二維氣體動力學方程組出發，提出了一種雙對稱收縮段設計理論，采用該收縮段向靜壓圓盤氣體軸承供氣，可以使圓盤間隙出口處的氣體流速達到超音速，軸承承載力可隨供氣壓力的提高線性增加，供氣壓力可以達到2.03 MPa(20 atm)以上。供氣壓力的提高使得圓盤間隙中的氣體壓力和流速迅速增加，軸承流道內氣體從低速發展到高亞音速、甚至超音速狀態，氣膜間隙內部各部分之間、以及氣膜與圓盤固壁之間均存在明顯的溫度差，軸承圓盤內部熱傳導與流道間隙內高亞音速或超音速氣流的對流換熱過程通過圓盤壁面緊密耦合在一起，是一個典型的共軛傳熱問題。

共軛傳熱是指包含固體內部熱傳導與流經該固體表面的流體的自由、強制或混合對流之間的相互作用的一種復雜傳熱過程[8]。共軛傳熱現象在科學和工程領域中大量存在，如核反應系統[9-12]，復合式氣膜沖擊冷卻[13-16]，高速飛行器[17-20]，以及油膜軸承等[21-22]。

目前共軛傳熱方法在工程領域有一定程度的應用[10-11，23-24]，受到越來越多的關注[25]。高壓氣體軸承流道間隙內氣流溫度變化對軸承圓盤的變形有直接的影響，因此研究軸承流道間隙內部和軸承圓盤的溫度場，以及流固耦合界面上的熱流密度分布規律是十分必要的。

本文作者使用ANSYS Workbench 2020R2工作平臺的Fluid Flow(Fluent)模塊，對高壓圓盤氣體軸承進行流體域與固體域的網格劃分，設置邊界條件并使用該模塊的Fluent程序進行共軛傳熱數值模擬，得到了該軸承共軛傳熱條件下流場參數分布、流體域與固體域的溫度分布及流固耦合壁面上的熱流密度分布，并將其與非共軛傳熱恒溫壁面條件下流體域的流場參數計算結果進行了對比，得到了高壓圓盤氣體軸承共軛傳熱的一些基本特性，為該類軸承的設計和制造提供了有益的指導。

1 軸承結構與邊界命名

1.1 軸承結構

高壓圓盤氣體軸承的結構如圖1所示。高壓圓盤氣體軸承由上、下工作圓盤組件組成，兩圓盤設計工作間隙h1為0.2 mm。上、下工作圓盤組件關于氣膜對稱面14呈平面對稱，同時關于對稱軸線1呈軸對稱。以上工作圓盤組件為例，其由上穩流管3、上密封墊圈5和上圓盤6通過螺紋連接緊固在一起，上穩流管3外端與供氣管道的管接頭相連。上穩流管3內部加工有軸對稱收縮段2，目的是對供氣管道輸送進來的氣體進行整流，使氣體以均勻的速度進入上供氣孔4。氣體從穩流管3、10進入，經供氣孔4、11和上、下圓盤收縮段7、13及平行氣膜區15，最終從氣膜間隙出口流出至環境大氣。上、下圓盤收縮段7、13的作用是使軸向來流無分離地轉變為平行于圓盤壁面的徑向均勻氣流，逐漸加速后再進入平行氣膜間隙。穩流管3、10的入口直徑d2為20 mm，高度h3為70 mm；供氣孔4、11的直徑d1為10 mm，圓盤直徑D為120 mm，圓盤厚度h2為30 mm。

圖1 高壓圓盤氣體軸承結構示意

文中采用Solidworks專業建模軟件對高壓圓盤氣體軸承進行三維建模，建模時忽略穩流管與工作圓盤之間的密封墊及穩流管端面倒角，以及工作圓盤上、下端面的倒角，并在氣膜間隙出口建立封閉面，方便抽取流體域。簡化后的高壓圓盤氣體軸承3D幾何模型如圖2所示。

圖2 高壓圓盤氣體軸承的簡化三維結構

1.2 軸承邊界命名

為了在保證計算精度的同時，盡量減少網格數量，以減少計算時間，文中采用氣體軸承三維模型的1/2作為計算域。在DesignModeler建模工具中導入Solidworks生成后綴名為x_t的模型文件，使用Fill功能對軸承內流域進行抽取，并利用布爾運算功能將穩流管和工作圓盤合并成一體。其流體域與固體域的邊界命名如圖3所示。

圖3 氣體軸承流體與固體計算域的邊界命名

2 網格劃分與邊界條件設置

2.1 網格劃分

在Fluid Flow(Fluent)模塊中，使用Mesh網格劃分工具，對高壓圓盤氣體軸承的流體域與固體域進行結構化網格劃分，其最終劃分網格結果如圖4所示。

圖4 三維軸承模型的固體域網格與圓盤收縮段附近的流體域局部網格

流體域網格劃分方法選擇Multizone，采用尺寸函數類型為Curvature，網格平滑度為中等，曲率法向角度為3°。對流體域近壁區添加5層膨脹層，膨脹層設置方法選擇為Smooth Transition，增長率設置為1.05，固體域網格劃分方法選擇Body sizing，且采用與流體域相同的單元尺寸。

2.2 單元基本尺寸的選擇

單元基本尺寸是控制模型網格數量的主要參數，基本尺寸越小，網格數量越大；單元基本尺寸的大小也會影響到網格質量，進而對計算結果產生影響。詳細分析不同單元基本尺寸算例的網格質量，對選擇合適的單元尺寸作為最終分析算例的設置參數具有指導作用。單元基本尺寸從1 mm變化到0.42 mm時，5種算例對應的網格數量和質量參數如表1所示。

表1 不同單元基本尺寸對應的網格數量和質量

從表1可知，單元基本尺寸從1 mm變化到0.42 mm時，網格數量從100萬增大到340萬；網格的平均質量、長寬比、雅可比值均變化較小，但歪斜率和扭曲因子變化較大；歪斜度從0.118 49減小到0.061 04，扭曲因子從2.378 1×10-4減小到1.121 6×10-4，兩者都幾乎減小了1/2。因此，表1中算例4和算例5的劃分參數是較為合適的參數。圖5進一步比較了5個算例在上圓盤流固耦合面(流體側)上的熱流密度。

圖5 不同網格數量時耦合面的熱流密度

由圖5可知，在算例4和算例5中，算例4的上圓盤流固耦合面(流體側)上的熱流密度計算值最大。對于具有相同導熱系數的同種材料，熱流密度越大，對應的溫度梯度也越大，圓盤內部溫度分布不均勻的程度也越大，進而其熱變形也越嚴重。因此，采用算例4的參數進行共軛傳熱求解，對后續的軸承圓盤變形分析、材料選擇和幾何參數設計而言，是保守的方案，得到的結果具有更大的安全裕量。故文中采用290萬網格的模型進行共軛傳熱求解。

2.3 邊界條件設置

如圖3所示，將上進氣口、下進氣口均設置為壓力入口，設置參數保持一致，即：入口總壓為2.03 MPa，入口初始壓力為2.02 MPa，湍流強度為5%，水力直徑為20 mm，入口總溫為300.15 K；將氣膜出口設置為壓力出口，出口壓力為0.10 MPa，湍流強度為5%，水力直徑為120 mm，出口總溫為300.15 K；將流體域中的上、下圓盤流固耦合面(流體側)與固體域中的上、下圓盤流固耦合面(固體側)設置為固定無滑移耦合面；將流體計算域對稱面、上圓盤對稱面、下圓盤對稱面均設置為對稱面；將上穩流管外側面、下穩流管外側面、上圓盤外表面、下圓盤外表面設置為與空氣自然對流的壁面，對流換熱系數取10 W/(m2·K)，自由流溫度取300.15 K；考慮到穩流管將與管接頭連接，故將上、下穩流管端面設置為絕熱壁面。

設置固體材料為Steel，其屬性保持默認；設置流體材料為Air，其密度采用理想氣體狀態方程，黏度采用三參數的薩瑟蘭公式，比熱容和導熱系數均采用溫度校正；設置操作壓力為0；采用SIMPLEC算法，其各項差分格式保持默認。能量松弛因子取0.3，其他松弛因子保持默認。

3 導熱方程與湍流模型

3.1 導熱方程

高壓圓盤氣體軸承工作時，一方面，流道內的氣流與圓盤壁面之間存在溫差，會產生對流換熱；另一方面，由于軸承圓盤內部存在著非均勻溫度場，固體圓盤內部也會產生熱傳導。流體計算域的能量方程[26]為

(1)

高壓圓盤氣體軸承固體計算域的導熱方程如下：

(2)

3.2 湍流模型

文中使用Realizablek-ε湍流模型，它是帶旋流修正的k-ε模型，對比標準的k-ε模型，該模型增加了一個與湍流黏性有關的公式，為湍流耗散率增加了新的傳輸方程。其優點之一是對平板和圓柱射流的發散比率有著更精確的預測。Realizablek-ε湍流模型的方程[26]如下：

Gb-ρε-YM

(3)

(4)

3.3 共軛傳熱問題的計算流程

Fluid Flow(Fluent)模塊中共軛傳熱問題的計算流程與單一流動問題類似，區別在于計算域同時包括流體域和固體域，設置邊界條件時，既要設置流體域的邊界條件，也要設置固體域的邊界條件。計算時對流體域和固體域進行整體迭代求解，得到耦合面的參數。計算流程如圖6所示。

圖6 共軛傳熱問題的計算流程

4 數值結果及分析

4.1 總體速度場和壓力場

圖7所示為高壓圓盤氣體軸承間隙流場的速度流線圖及圓盤收縮段附近的流場速度矢量圖。由圖7(a)可知，氣體從上、下穩流管入口流入，流經穩流管收縮段、供氣孔、圓盤收縮段及平行氣膜區，最終在氣膜間隙出口處流出至環境大氣。穩流管內的氣流速度普遍較低，而平行氣膜區氣流速度普遍較高，氣流從進入穩流管開始逐漸加速，并在工作圓盤收縮段出口處達到極大值，而后在平行氣膜區速度有所降低，在氣膜間隙出口處，氣流速度再次加速達到最大值。盡管進入上、下進氣口的氣流都存在一定的湍流度，但在上、下穩流管軸對稱收縮段的整流作用下，氣流形成了均勻的流線從上、下供氣孔進入流道間隙。由圖7(b)可以直觀地觀察到，流道中心區域的氣流速度矢量分布及圓盤收縮段壁面邊界層的速度梯度變化，此處氣流經供氣孔流入，在上、下圓盤收縮段的作用下，將軸向來流無分離地轉變為徑向流動，并沿半徑R方向逐漸加速，最后從氣膜間隙出口流出。

圖7 氣體軸承間隙流場的速度流線圖及圓盤收縮段附近的速度矢量分布

使用Fluent軟件自帶的Iso-surface截線工具可以截取高壓圓盤氣體軸承的氣膜對稱線，并提取對稱線上的數據。如圖8所示為高壓圓盤氣體軸承氣膜對稱線上的馬赫數和靜壓分布。

圖8 氣膜對稱線上的馬赫數和靜壓分布

由圖8中的馬赫數曲線可知：氣流馬赫數在供氣孔區域a從接近于0逐漸增大至0.11，并在收縮段區域b快速增大至超音速1.13；而后又在平行氣膜區c先緩慢下降，平行氣膜區c的馬赫數整體維持在0.86左右的高亞音速狀態；最后在接近氣膜間隙出口時由于壓力遠大于環境背壓，從而使馬赫數快速增加到超音速1.3。

由圖8中的靜壓曲線可知：氣流壓力在供氣孔區域a緩慢下降至1.95 MPa，而后在圓盤收縮段區域b快速下降至0.827 MPa，在到達平行氣膜區c后，氣流壓力經過小幅上升后再次連續下降至氣膜出口處的0.349 MPa。因為超音速流動的基本特性是下游的低壓擾動不會影響上游的壓力分布，故外界環境的低壓擾動不會影響軸承流道內的壓力分布，所以在高壓圓盤氣體軸承流道內的絕大部分區域，壓力能維持在0.592 MPa以上，較好地實現了高壓氣體潤滑。

4.2 傳熱分析

4.2.1 溫度場分析

如圖9所示為軸承對稱面上的溫度云圖，軸承對稱面由流體計算域對稱面、上圓盤對稱面、下圓盤對稱面組合而成。

圖9 軸承對稱面上的溫度云圖

在來流總溫不變的條件下，馬赫數越大，溫度越低。結合圖8、9可知，供氣孔區域氣流速度較低，該區域的溫度與環境溫度相比變化很小，故上、下穩流管的溫度變化不明顯；圓盤收縮段出口處及氣膜間隙出口處的馬赫數較高，相應的溫度也較低；而在平行氣膜區絕大部分區域氣流仍然保持在0.86馬赫以上的高亞音速狀態，故整個軸承圓盤在鄰近平行氣膜區的部分溫度較低，而遠離平行氣膜區的部分溫度較高。

圖10所示為與軸承對稱面相垂直的截面上的溫度等值線圖。

圖10 與軸承對稱面相垂直的截面上溫度等值線

由圖10可知，上、下圓盤的溫度以氣膜對稱面為對稱平面呈對稱式分布，圓盤收縮段出口處及氣膜間隙出口處的溫度最低，這兩處周圍的等溫線較為密集，溫度梯度較大，導熱熱流密度也較大。軸承圓盤內的溫度梯度矢量以收縮段出口處及氣膜間隙出口處為中心，向圓盤外端面及穩流管方向逐漸減小。熱流密度矢量也相應減小，其方向與溫度梯度矢量的方向相反。

氣膜對稱線及流體域耦合壁面上的溫度分布如圖11所示。

圖11 氣膜對稱線及流體域耦合壁面上的溫度分布

由圖11可知，在氣膜對稱線上，氣流在供氣孔區域a的溫度變化并不明顯，而在圓盤收縮段區域b的前半部分變化較小，后半部分變化較大，收縮段出口處的溫度已陡降至241 K。而后在平行氣膜區c由于氣流的速度有所下降，其溫度有明顯上升，在平行氣膜區c的絕大部分區域，溫度都維持在260 K以上。當氣流到達氣膜間隙出口處時，由于速度再次增大至超音速，其溫度驟降至235 K。圖11中流體域耦合壁面上的溫度變化趨勢與氣膜對稱線上氣流溫度的分布趨勢相類似，但耦合壁面在工作圓盤收縮段出口處及氣膜出口處的溫度明顯高于氣膜對稱線上相應位置的氣流溫度，說明高亞音速特別是超音速狀態下圓盤壁面溫度邊界層中可以存在較高的溫度差。總的說來，共軛傳熱條件下，流固耦合壁面上的溫度分布仍然取決于間隙內的氣流速度分布，即軸承初始的速度型設計。

4.2.2 熱流密度分析

由于耦合的流體域壁面與固體域壁面上的熱流密度大小相等方向相反，故僅分析流體域耦合壁面的熱流密度即可。如圖12所示為高壓圓盤氣體軸承流體域耦合壁面上的熱流密度分布。

圖12 流體域耦合壁面上的熱流密度分布

由圖12可知，圓盤收縮段區域b的流體域耦合壁面上熱流密度先小幅下降后再上升，然后再經過平緩的下降后再明顯上升，并在收縮段區域b的出口處達到峰值，其值為73 467.1 W/m2，且收縮段區域 b的部分區域流體壁面熱流密度為負，說明在收縮段區域b內的不同位置，同時存在流體吸熱和放熱的現象，這是由于軸承圓盤的熱傳導與間隙內氣體對流換熱具有的復雜共軛作用機制造成的。平行氣膜區c的換熱達到熱充分發展狀態，氣流吸熱，熱流密度為正。氣膜間隙出口處的流體域壁面熱流密度達到最大值135 966 W/m2。

4.3 共軛傳熱與非共軛傳熱的對比分析

為了研究高壓圓盤氣體軸承共軛傳熱解與非共軛傳熱解的差異性，以圖3(a)為研究模型，獲取高壓圓盤氣體軸承在恒溫壁面條件下的計算結果，此時整個計算域只包含單一的流體域，可看作是一類非共軛傳熱解。即先在Fluid Flow(Fluent)模塊的Mesh界面中，將單一流體域模型輸出成后綴名為msh的網格文件，然后單獨啟動Fluent 2020R2軟件，導入前述的msh網格文件；將流體域上、下圓盤壁面設置為固定無滑移恒溫壁面，壁溫300.15 K，其余設置同2.3節。

將單一流體域的計算結果與共軛傳熱條件下的計算結果進行對比，圖13所示為氣膜對稱線上的馬赫數、靜壓和靜溫對比。可知，共軛傳熱條件下高壓圓盤氣體軸承氣膜對稱線上的馬赫數、靜壓、靜溫和非共軛傳熱恒溫壁面條件下的變化趨勢是一致的，但在局部區域數值結果存在明顯差異。由圖13(a)可知，在供氣孔區域a和圓盤收縮段區域b的絕大部分區域，兩者僅存在微小差異，但從收縮段出口至平行氣膜區c的絕大部分區域，兩者存在明顯的差異；其中最大差異發生在圓盤收縮段出口處，兩者相差達到10.76%。

圖13 氣膜對稱線上的流場參數對比

由圖13(b)可知，共軛傳熱和非共軛傳熱恒溫壁面條件下軸承氣膜對稱線上的靜壓在供氣孔區域a、圓盤收縮段區域b和平行氣膜區c均存在較為明顯的差異，且最大差異發生在收縮段出口處，兩者相差15.78%。由圖13(c)可知，共軛傳熱和非共軛傳熱恒溫壁面條件下氣膜對稱線上的靜溫在供氣孔區域a和圓盤收縮段區域b的絕大部分區域存在微小差異，但在收縮段出口處及平行氣膜區c均存在較大差異，其中收縮段出口處溫度差為13.05 K，相對誤差為5.14%；最大差異發生在氣膜間隙出口處，溫度差為18.14 K，相對誤差為7.16%。綜上，造成2種條件下流場參數計算結果存在明顯差異的原因，在于2種情況下流固耦合壁面上的溫度分布不同，進而使軸承間隙流道內，尤其是平行氣膜區的馬赫數發生變化，進而改變了間隙內的壓力和溫度分布。

如圖14所示為高壓圓盤氣體軸承共軛傳熱與非共軛傳熱恒溫壁面條件下流體域耦合壁面上的熱流密度的對比。可知，2種情況下流體域耦合壁面上熱流密度的變化趨勢相類似，但圓盤收縮段b及平行氣膜區c的數值結果均存在較大差異，其中氣膜間隙出口處的熱流密度差異為60.06%，最大差異發生在收縮段出口處，兩者相差79.28%。非共軛傳熱恒溫壁面條件下，間隙內的氣體只吸熱，流體域耦合壁面上的熱流密度均為正值。

圖14 流體域耦合壁面的熱流密度對比

共軛傳熱和非共軛傳熱恒溫壁面條件下，流固耦合壁面上熱流密度計算結果產生顯著變化的主要原因，在于恒溫壁面假設在很大程度上高估了圓盤壁面氣膜側的溫度。由于在流固耦合面附近，軸承圓盤的熱傳導與氣膜內氣體的對流換熱發生了復雜的共軛作用機制，這與恒溫壁面簡化下的傳熱機制有本質的不同。軸承在實際工作中，流固耦合壁面上的溫度不可能保持恒溫，因此恒溫壁面假設使計算誤差偏大，共軛傳熱計算可以得到更為符合實際的結果[23]。

5 結論

(1)由于存在著復雜的共軛作用機制，高壓圓盤氣體軸承流固耦合壁面上的熱流密度及溫度分布曲線，均是隨R的增大發生劇烈變化的曲線；流固耦合壁面上既非恒熱流密度，也非恒壁溫，這是非共軛傳熱視角下難以獲得的結果；高壓圓盤氣體軸承的傳熱問題，是一個典型的共軛傳熱問題。

(2)共軛傳熱時，高壓圓盤氣體軸承流固耦合壁面上溫度分布曲線的變化趨勢，與氣膜對稱線上溫度變化趨勢類似，但二者在圓盤收縮段出口處和氣膜間隙出口處具體溫度值有明顯差異；因此，總的說來，流固耦合壁面上的溫度分布仍然取決于間隙內氣流速度分布，即軸承初始的速度型設計。

(3)共軛傳熱與非共軛傳熱恒溫壁面條件下得到的氣膜對稱線上的馬赫數、靜壓、靜溫以及流固耦合壁面上的熱流密度、溫度均存在明顯差異。非共軛傳熱恒溫壁面條件下，間隙內的氣體只吸熱，流體域耦合壁面上熱流密度均為正值；而共軛傳熱條件下流體域耦合壁面上熱流密度存在正負值，即間隙內氣體在流動過程中，吸熱和放熱同時存在，顯示出軸承圓盤的熱傳導與間隙內氣體的對流換熱具有復雜的共軛作用機制。

(4)分析高壓圓盤氣體軸承的溫度場，推薦采用共軛傳熱方法進行計算，可以得到更為符合實際的結果。