顧延?xùn)|，成 立，B?HLE Martin，SCHIMPF Artur，袁壽其

(1.揚(yáng)州大學(xué) 水利科學(xué)與工程學(xué)院，江蘇 揚(yáng)州 225009；2.江蘇大學(xué) 國(guó)家水泵及系統(tǒng)工程技術(shù)研究中心，江蘇 鎮(zhèn)江 212013；3.凱澤斯勞滕工業(yè)大學(xué) 流體力學(xué)和流體機(jī)械研究所，德國(guó)凱澤斯勞滕 67663)

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

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

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

1 多孔質(zhì)軸承參數(shù)

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

圖1 多孔質(zhì)軸承Fig.1 Porous bearing

表1 軸承幾何參數(shù)及運(yùn)行工況Tab.1 Geometric parameters and working conditions

2 理論建模與數(shù)值計(jì)算

2.1 雷諾潤(rùn)滑方程

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

圖2 多孔質(zhì)軸承-轉(zhuǎn)子系統(tǒng)示意圖Fig.2 Diagrammatic representation for porous bearing-rotor system

(1)

式中：x，y和z分別是周向、軸向和徑向坐標(biāo)；h是徑向間隙函數(shù)[16]；μ是動(dòng)力黏度；u，v和w分別是周向、軸向和徑向速度，下標(biāo)a和b表示軸頸和軸襯上的分量。接下來(lái)簡(jiǎn)化式(1)。

假設(shè)軸頸為無(wú)滑移壁面，考慮剪切和擠壓效應(yīng)，則軸頸上的速度邊界條件是

ua=R1ωcosβ≈R1ω

(2)

va=0

(3)

(4)

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

(5)

2.2 多孔質(zhì)軸襯流動(dòng)模型

基于多孔質(zhì)達(dá)西方程和連續(xù)方程，建立多孔質(zhì)軸襯的流動(dòng)模型。假設(shè)多孔質(zhì)軸襯是三維層流，多孔質(zhì)各項(xiàng)同性并忽略慣性阻力，采用達(dá)西方程描述其內(nèi)部流動(dòng)

(6)

(7)

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

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

(8)

需要強(qiáng)調(diào)的是，式(8)也適用于多孔質(zhì)-固體交界面。

2.3 潤(rùn)滑膜流動(dòng)模型

基于雷諾潤(rùn)滑方程和多孔質(zhì)達(dá)西方程，建立多孔質(zhì)-潤(rùn)滑膜交界面及動(dòng)壓潤(rùn)滑膜的流動(dòng)模型。對(duì)于多孔質(zhì)-流體交界面，既有無(wú)滑移壁面處理方法，也有多種滑移壁面處理方法[17]。假設(shè)該交界面為滑移壁面，周向和軸向速度由達(dá)西方程控制，即

(9)

(10)

徑向速度為

(11)

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

(12)

非多孔質(zhì)區(qū)域的壁面條件是

ub_non-porous=0

(13)

vb_non-porous=0

(14)

wb_non-porous=0

(15)

將式(13)～式(15)代入式(5)，得動(dòng)壓潤(rùn)滑膜控制方程

(16)

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

2.4 邊界條件

上述橢圓型偏微分方程的未知量是壓力，因而在邊界條件設(shè)置中，需給定壓力值或壓力梯度或兩者組合，即Dirichlet(Ⅰ)、Neumann(Ⅱ)及Robin(Ⅲ)邊界條件。實(shí)際應(yīng)用中，在多孔質(zhì)軸襯進(jìn)口供給一定壓力的潤(rùn)滑劑，再?gòu)妮S承兩端排出到工作環(huán)境，因而設(shè)置進(jìn)出口為壓力Dirichlet邊界條件。在多孔質(zhì)-固體交界面上無(wú)工作介質(zhì)通過(guò)，交界面法向速度(軸向速度)為0，即

(17)

則有

(18)

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

2.5 網(wǎng)格劃分及離散格式

采用數(shù)值方法對(duì)式(8)、式(12)和式(16)同步求解?；谟邢薏罘址ǎ瑢?duì)計(jì)算域劃分如圖4所示的網(wǎng)格。在潤(rùn)滑膜和多孔質(zhì)的周向及軸向上，采用均勻節(jié)點(diǎn)分布。在多孔質(zhì)徑向上，采用間距等比增長(zhǎng)的節(jié)點(diǎn)分布，即Δzk+1=qΔzk。

圖4 PBS計(jì)算網(wǎng)格Fig.4 Mesh for PBS

對(duì)控制方程中的周向和軸向一階、二階導(dǎo)數(shù)采用中心差分格式

(19)

(20)

(21)

(22)

對(duì)式(12)中的徑向一階導(dǎo)數(shù)采用三節(jié)點(diǎn)向前差分格式[18]

(23)

對(duì)式(8)中的徑向一階、二階導(dǎo)數(shù)采用中心差分格式

(24)

(25)

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

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

(26)

pvirtual=pi,j-1,k

(27)

(28)

3 PBS與Fluent結(jié)果對(duì)比

3.1 Fluent計(jì)算設(shè)置

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

圖5 Fluent計(jì)算網(wǎng)格Fig.5 Mesh for Fluent

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

3.2 PBS計(jì)算設(shè)置

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

(29)

(30)

(31)

(32)

通過(guò)網(wǎng)格無(wú)關(guān)性分析，確定3個(gè)軸承的網(wǎng)格節(jié)點(diǎn)數(shù)均為192 000。

3.3 PBS與Fluent結(jié)果對(duì)比

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

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

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

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

圖7 PBS與Fluent的潤(rùn)滑膜壓力對(duì)比Fig.7 Comparisons of pressure in lubricating film between PBS and Fluent

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

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

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

通過(guò)對(duì)比PBS與Fluent，初步驗(yàn)證了多孔質(zhì)徑向靜壓軸承的理論模型與數(shù)值方法是準(zhǔn)確可靠的。

4 結(jié)果與討論

4.1 供給壓差對(duì)靜特性的影響

以CaseA為例，圖10是供給壓差對(duì)承載能力的影響。在同一偏心率下，承載能力隨著供給壓差增大而線性增大，比如供給壓差2 bar的承載能力是1 bar的2倍左右，4 bar是1 bar的4倍左右。在同一供給壓差下，承載能力隨著偏心率增大而增大。圖11是供給壓差對(duì)偏位角的影響。在同一偏心率下，偏位角隨著供給壓差增大而減小。這是因?yàn)槠唤怯蓜?dòng)壓和靜壓共同決定，靜壓隨著供給壓差增大，而動(dòng)壓幾乎不變，因而動(dòng)壓與靜壓比值減小、偏位角減小。在同一供給壓差下，偏心率小于0.7時(shí)偏位角幾乎不變，偏心率大于0.7時(shí)偏位角略微增大。

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

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

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

圖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倍左右。在同一供給壓差下，流量和功耗最小值都出現(xiàn)在零偏心下。隨著偏心率增大，流量和功耗逐漸增大。

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

4.2 潤(rùn)滑劑黏度對(duì)靜特性的影響

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

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

(35)

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

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

(36)

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

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

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

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

5 結(jié) 論

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

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