多孔質靜壓徑向軸承的理論建模與數值計算

顧延東，成 立，B?HLE Martin，SCHIMPF Artur，袁壽其

(1.揚州大學 水利科學與工程學院，江蘇 揚州 225009；2.江蘇大學 國家水泵及系統工程技術研究中心，江蘇 鎮江 212013；3.凱澤斯勞滕工業大學 流體力學和流體機械研究所，德國凱澤斯勞滕 67663)

靜壓徑向軸承主要依靠供壓源和節流器產生承載能力，具有承載能力大、啟停性能好等諸多優點，常應用于旋轉機械中。節流器是其關鍵零件之一，有多孔質式、小孔式等。其中，利用多孔質材料加工出具有節流降壓作用的軸襯，再設計供壓和封裝部分，就可制造出多孔質軸承。與小孔式軸承相比，多孔質軸承具有更大的承載能力、剛度等。

由于微米級潤滑間隙、軸系對中、軸系旋轉等原因，很難開展多孔質靜壓徑向軸承實驗研究，尤其是測試偏心率-承載能力等靜特性[1-2]。目前，建立理論模型并數值求解是獲得軸承特性的主要方法[3-4]。Sneck等[5]假設多孔質軸襯為一維徑向可壓流動，對多孔質達西方程徑向積分并代入雷諾潤滑方程，得到二維非線性偏微分方程，采用攝動法求解小偏心率下的靜特性。Sun[6]先用Newton-Raphson法對上述方程線性化處理，再用逐次超松弛法(successive over relaxation,SOR)加速求解，指出多孔質滲透率較大時，不僅降低承載能力，而且加劇供氣功耗。Majumder等[7]建立了多孔質三維流動方程并與氣膜雷諾潤滑方程耦合，分析了多孔質-氣膜交界面速度滑移對氣動不穩定性及臨界質量的影響。Lee等[8]假設多孔質-氣膜交界面上無速度滑移，使用低松弛因子保證可壓計算的穩定性，給出了寬徑比等參數的設計建議。馮凱等[9]考慮溫度對氣體物性及材料變形的影響，建立了帶限制層多孔質軸承的溫度模型，指出轉速對溫升的影響比載荷大。楊國棟等[10]忽略溫度影響，使用非結構網格和有限體積法準確描述帶槽液體潤滑軸承的幾何特征，獲得了寬徑比、偏心率對液膜壓力的影響規律。商業通用CFD(computational fluid dynamics)軟件是預測多孔質軸承靜特性的另一方法[11-12]，但較難完成嚴謹合格的數值模擬，而且版權費高、計算時間長。在工程設計中，往往需要快速查找或計算零部件靜特性，定制CFD程序恰好滿足該需求。

本研究的主要目的是開發一種準確、快速、穩定的多孔質軸承靜特性計算程序。首先，建立多孔質、多孔質-潤滑膜交界面及動壓潤滑膜的理論模型。然后，采用數值方法求解理論模型，使用C++編寫求解器PBS，并與ANSYS Fluent對比。最后，利用PBS分析供給壓差和潤滑劑黏度對多孔質軸承靜特性的影響。

1 多孔質軸承參數

針對超臨界氣體旋轉機械等特殊裝備，采用碳纖維增強碳基復合材料[13]制造多孔質軸承并利用超臨界氣體潤滑，如圖1所示。超臨界氣體的密度與液體相近、黏度又明顯小于液體，是值得探索的潤滑劑。設計了3個多孔質軸承，相關參數如表1和圖2所示。以下計算均基于表1參數。

圖1 多孔質軸承Fig.1 Porous bearing

表1 軸承幾何參數及運行工況Tab.1 Geometric parameters and working conditions

2 理論建模與數值計算

2.1 雷諾潤滑方程

圖2是多孔質軸承-轉子系統示意圖。采用雷諾潤滑方程(reynolds lubrication equation)描述多孔質軸承潤滑膜流動，該方程可以從簡化Navier-Stokes(N-S)方程、再代入連續方程積分得到，或直接從黏性流體本構方程和連續方程推得[14-15]，其一般完整形式為

圖2 多孔質軸承-轉子系統示意圖Fig.2 Diagrammatic representation for porous bearing-rotor system

(1)

式中：x，y和z分別是周向、軸向和徑向坐標；h是徑向間隙函數[16]；μ是動力黏度；u，v和w分別是周向、軸向和徑向速度，下標a和b表示軸頸和軸襯上的分量。接下來簡化式(1)。

假設軸頸為無滑移壁面，考慮剪切和擠壓效應，則軸頸上的速度邊界條件是

ua=R1ωcosβ≈R1ω

(2)

va=0

(3)

(4)

在等溫、不可壓及層流條件下，潤滑劑物性參數不變，結合上述邊界條件式(2)～式(4)，式(1)可簡化為

(5)

2.2 多孔質軸襯流動模型

基于多孔質達西方程和連續方程，建立多孔質軸襯的流動模型。假設多孔質軸襯是三維層流，多孔質各項同性并忽略慣性阻力，采用達西方程描述其內部流動

(6)

(7)

式中：α是滲透率；i是自由標；j是啞標。將式(6)代入不可壓連續方程式(7)并應用到圖3坐標系，得多孔質軸襯控制方程，即圓柱坐標系下的壓力Laplace方程

圖3 控制方程及邊界條件作用域Fig.3 Scopes for equations and boundary conditions

(8)

需要強調的是，式(8)也適用于多孔質-固體交界面。

2.3 潤滑膜流動模型

基于雷諾潤滑方程和多孔質達西方程，建立多孔質-潤滑膜交界面及動壓潤滑膜的流動模型。對于多孔質-流體交界面，既有無滑移壁面處理方法，也有多種滑移壁面處理方法[17]。假設該交界面為滑移壁面，周向和軸向速度由達西方程控制，即

(9)

(10)

徑向速度為

(11)

將式(9)～式(11)代入式(5)，得多孔質-潤滑膜交界面控制方程

(12)

非多孔質區域的壁面條件是

ub_non-porous=0

(13)

vb_non-porous=0

(14)

wb_non-porous=0

(15)

將式(13)～式(15)代入式(5)，得動壓潤滑膜控制方程

(16)

式(8)、式(12)和式(16)構成了多孔質靜壓徑向軸承的流動控制方程組，作用域如圖3所示。

2.4 邊界條件

上述橢圓型偏微分方程的未知量是壓力，因而在邊界條件設置中，需給定壓力值或壓力梯度或兩者組合，即Dirichlet(Ⅰ)、Neumann(Ⅱ)及Robin(Ⅲ)邊界條件。實際應用中，在多孔質軸襯進口供給一定壓力的潤滑劑，再從軸承兩端排出到工作環境，因而設置進出口為壓力Dirichlet邊界條件。在多孔質-固體交界面上無工作介質通過，交界面法向速度(軸向速度)為0，即

(17)

則有

(18)

式(17)為壓力Neumann邊界條件。也就是說，多孔質-固體交界面上的壓力被式(8)和式(18)同時控制。多孔質軸承的各邊界條件類型如圖3所示。

2.5 網格劃分及離散格式

采用數值方法對式(8)、式(12)和式(16)同步求解。基于有限差分法，對計算域劃分如圖4所示的網格。在潤滑膜和多孔質的周向及軸向上，采用均勻節點分布。在多孔質徑向上，采用間距等比增長的節點分布，即Δzk+1=qΔzk。

圖4 PBS計算網格Fig.4 Mesh for PBS

對控制方程中的周向和軸向一階、二階導數采用中心差分格式

(19)

(20)

(21)

(22)

對式(12)中的徑向一階導數采用三節點向前差分格式[18]

(23)

對式(8)中的徑向一階、二階導數采用中心差分格式

(24)

(25)

式(19)～式(24)均為二階精度，式(25)的精度在一階到二階之間，取決于公比q。對比了網格增長比率1(均勻分布，二階精度)和1.1(近似二階精度)的計算結果，壓力幾乎沒有區別，因而在后續計算中都采用了1.1增長率的徑向網格。交界面處細化網格對共軛傳熱計算非常重要，這在以后的研究中作詳細說明。

多孔質-固體交界面上的壓力只需式(18)的單向差分格式就可計算，這也是現有文獻的處理方法。但式(18)的二階精度、單向差分格式收斂較慢，潤滑膜壓力分布不連續，使用逐次超松弛法時易發散，制約了整個迭代計算的收斂速度及穩定性。如2.4節所述，多孔質-固體交界面上的壓力被式(8)(Laplace方程)和式(18)(Neumann邊界條件)同時控制，這里提出一種Laplace-Neumann虛擬節點法：①采用虛擬節點法處理Neumann邊界條件[19]，其節點構造如圖4(c)所示，利用虛擬節點pvirtual、內部節點pi,j-1,k及中心差分格式離散式(18)，得式(26)及式(27)，該處理方法為二階精度；②采用式(21)、式(22)、式(24)及式(25)并結合虛擬節點對式(8)進行離散，然后用式(27)替換虛擬節點，則式(8)的軸向二階導數離散形式寫成式(28)并代入迭代矩陣。這個方法收斂速度極快、殘差極小，更重要的是壓力分布連續，符合預期的物理現象。基于上述方法，使用C++編寫了快速穩定的求解器PBS.exe。

(26)

pvirtual=pi,j-1,k

(27)

(28)

3 PBS與Fluent結果對比

3.1 Fluent計算設置

對比商業通用CFD軟件ANSYS Fluent 2019 R1的模擬結果，初步驗證PBS理論模型和數值方法的準確性。在ANSYS Workbench中，首先使用Design Modeler建立計算域的幾何模型。然后，使用Mesh模塊的MultiZone策略劃分六面體網格，如圖5所示。在偏心方向進行局部加密，交界面網格節點完全對齊，保證了良好的網格正交性、長寬比等。通過網格無關性分析，確定3個軸承的網格節點數均為15 949 060。最后，使用Fluent對不同偏心率(E/h0)下的多孔質軸承進行模擬。在Fluent設置中，使用層流模型，采用壓力進口、壓力出口等邊界條件，使用SIMPLE算法，各項離散精度均為二階，收斂殘差為1.0×10-6。

圖5 Fluent計算網格Fig.5 Mesh for Fluent

PBS與Fluent主要差別為：①Fluent是有限體積法，PBS是有限差分法；②Fluent潤滑膜的控制方程是三維N-S方程，PBS潤滑膜是二維雷諾潤滑方程；③Fluent中的多孔質-潤滑膜交界面不能進行滑移處理，而一定條件下滑移對動壓效應影響顯著。在PBS中可以植入多種滑移模型，本文采用的是達西方程控制的滑移邊界條件；④Fluent是將多孔質達西方程作為阻力源項，直接代入完整的N-S方程求解，滲透率大到一定程度，擴散項和達西阻力項同等顯著，多孔質流動甚至從Stokes流轉變為對流-擴散流動。而在PBS中，只用達西方程和連續方程描述多孔質流動，適用于Stokes流。因此，在大滲透率下Fluent和PBS不再有可比性。由于多孔質軸承是小滲透率，可對比兩者。

3.2 PBS計算設置

在PBS設置中，采用迭代矩陣余量的均方根作為收斂準則，收斂精度為1.0×10-6。利用逐次超松弛法加速收斂，松弛因子為1.92。選取偏心率0.9下的承載能力(L)和偏位角(A)做網格無關性分析(坐標系如圖3所示)，計算公式為

(29)

(30)

(31)

(32)

通過網格無關性分析，確定3個軸承的網格節點數均為192 000。

3.3 PBS與Fluent結果對比

圖6對比了PBS與Fluent的承載能力(L)及偏位角(A)。總體而言，PBS與Fluent具有很好的一致性。在承載能力方面，PBS比Fluent高，大偏心率下(e>0.5)差值維持在2.6 N內，產生偏差的主要原因如3.1節所述。在偏位角方面，PBS比Fluent略大。動壓相對于靜壓較小，造成偏位角較低。更重要的是，與商業通用CFD軟件相比，定制CFD程序PBS具有即時仿真(instant simulation)、計算硬件要求低、可操作性強等優點，滿足設計人員“one click，one result”需求。以Intel E5-2696V4 CPU為例，PBS 1核計算多孔質軸承靜特性曲線需10 s左右，而Fluent 20核需4.68×104s左右。

圖6 PBS與Fluent的承載能力及偏位角對比Fig.6 Comparisons of load capacity and azimuthal angle between PBS and Fluent

為進一步驗證PBS，以CaseA偏心率0.9為例，對比PBS和Fluent的潤滑膜壓力場，如圖7所示。為方便對比，在ANSYS CFD-Post中輸出Fluent結果到MATLAB作圖。PBS和Fluent的潤滑膜壓力分布高度相似。PBS壓力最大值是4.006 8 bar，Fluent是4.005 2 bar，兩者相差很小。由于動壓效應，壓力最大值高于供給壓力。PBS和Fluent壓力最小值均在出口，即壓力邊界條件1 bar。PBS高壓區面積略大于Fluent，這是圖6中PBS承載能力高于Fluent的主要原因之一。

圖8對比了PBS與Fluent的多孔質軸襯中截面上的壓力，兩者非常相似。壓力最大值出現在偏心方向上，PBS壓力最大值是4.006 8 bar，Fluent是4.005 2 bar，與圖7一致。壓力最小值出現在偏心反方向上，PBS壓力最小值是2.737 5 bar，Fluent是2.725 2 bar，兩者相差很小。

圖7 PBS與Fluent的潤滑膜壓力對比Fig.7 Comparisons of pressure in lubricating film between PBS and Fluent

圖9對比了PBS與Fluent的多孔質軸襯中截面上的徑向速度，與上述對比一樣，徑向速度分布也很相似。高速區出現在偏心反方向內層處，PBS和Fluent徑向速度最大值均為0.159 7 m/s。在偏心方向上，徑向速度趨于0，而且最低速度區略微順時針偏移，與圖8高壓區偏移一致，這是由動壓效應造成的。

圖8 PBS與Fluent的多孔質壓力對比Fig.8 Comparisons of pressure in porous restrictor between PBS and Fluent

圖9 PBS與Fluent的多孔質徑向速度對比Fig.9 Comparisons of radial velocity in porous restrictor between PBS and Fluent

通過對比PBS與Fluent，初步驗證了多孔質徑向靜壓軸承的理論模型與數值方法是準確可靠的。

4 結果與討論

4.1 供給壓差對靜特性的影響

以CaseA為例，圖10是供給壓差對承載能力的影響。在同一偏心率下，承載能力隨著供給壓差增大而線性增大，比如供給壓差2 bar的承載能力是1 bar的2倍左右，4 bar是1 bar的4倍左右。在同一供給壓差下，承載能力隨著偏心率增大而增大。圖11是供給壓差對偏位角的影響。在同一偏心率下，偏位角隨著供給壓差增大而減小。這是因為偏位角由動壓和靜壓共同決定，靜壓隨著供給壓差增大，而動壓幾乎不變，因而動壓與靜壓比值減小、偏位角減小。在同一供給壓差下，偏心率小于0.7時偏位角幾乎不變，偏心率大于0.7時偏位角略微增大。

圖10 不同供給壓差下的承載能力Fig.10 Load capacity under different feeding pressure differences

圖11 不同供給壓差下的偏位角Fig.11 Azimuthal angle under different feeding pressure differences

圖12和13分別是供給壓差對供給流量(F)和供給功耗(P)的影響。計算公式為

圖12 不同供給壓差下的流量Fig.12 Flowrate under different feeding pressure differences

(33)

P=FΔp

(34)

在同一偏心率下，隨著供給壓差增大，流量線性增大，功耗平方增大，比如供給壓差4 bar與1 bar相比，壓差比是4倍，4 bar的流量是1 bar的4倍左右，功耗是16倍左右。在同一供給壓差下，流量和功耗最小值都出現在零偏心下。隨著偏心率增大，流量和功耗逐漸增大。

圖13 不同供給壓差下的功耗Fig.13 Power consumption under different feeding pressure differences

4.2 潤滑劑黏度對靜特性的影響

以CaseA為例，黏度設置為2.20×10-5Pa·s、2.99×10-5Pa·s、3.66×10-5Pa·s及4.27×10-5Pa·s。圖14是潤滑劑黏度對承載能力的影響。在同一偏心率下，隨著黏度增大，承載能力略微增大，偏位角近似線性增大。這是因為：①多孔質壓力控制方程式(8)是Laplace方程，各項約去了黏度和滲透率，造成多孔質壓力分布完全依賴于邊界條件；②將多孔質-潤滑膜交界面控制方程式(12)的兩邊乘以常數黏度，得

圖14 不同潤滑劑黏度下的承載能力Fig.14 Load capacity under different lubricant viscosity

(35)

圖15是潤滑劑黏度對偏位角的影響。對于該多孔質軸承，靜壓在承載能力中絕對主導，式(32)tan(FX/FY)很小時，有

圖15 不同潤滑劑黏度下的偏位角Fig.15 Azimuthal angle under different lubricant viscosity

(36)

所以偏位角隨著黏度增大而近似線性增大。

盡管黏度對承載能力影響很小，但對功耗和流量影響很大。圖16和17分別是潤滑劑黏度對供給流量和功耗的影響。在同一偏心率下，隨著黏度增大，流量和功耗線性減小。比如黏度4.27×10-5Pa·s與2.20×10-5Pa·s相比，黏度比是1.94，流量比和功耗比是0.52。上文從控制方程的角度分析了黏度對承載能力和偏位角的影響，即黏度對徑向壓力梯度影響較小，再結合式(33)和式(34)可知，流量和功耗均與黏度成反比，因此隨著黏度增大，流量和功耗線性減小。與供給壓差對流量和功耗的影響一樣，在同一黏度下，流量和功耗最小值都出現在零偏心下。隨著偏心率增大，流量和功耗逐漸增大。

圖16 不同潤滑劑黏度下的流量Fig.16 Flowrate under different lubricant viscosity

圖17 不同潤滑劑黏度下的功耗Fig.17 Power consumption under different lubricant viscosity

5 結 論

(1)針對超臨界氣體旋轉機械等特殊裝備，采用碳纖維增強碳基復合材料制造多孔質軸承并探索超臨界氣體潤滑。基于雷諾潤滑方程、達西方程和連續方程，考慮多孔質-潤滑膜交界面速度滑移，建立了多孔質軸承的理論模型。引入有限差分法、變步長差分格式、逐次超松弛法，提出Laplace-Neumann虛擬節點法，建立了該理論模型的數值求解方法。使用C++編寫了多孔質軸承求解器PBS，PBS與Fluent結果相近，但PBS實現了即時仿真。

(2)PBS結果表明：隨著供給壓差增大，承載能力線性增大，偏位角減小，流量線性增大，功耗平方增大；隨著潤滑劑黏度增大，承載能力略微增大，偏位角增大，流量和功耗線性減小；隨著偏心率增大，承載能力、流量及功耗增大。