動壓軸承流固界面非一致滑移流動分析方法及分布特征

張鏡洋，陳哲，呂元偉，孫義建，張靖周，陳衛東，3，羅欣洋

1.南京航空航天大學 航天學院，南京 211106

2.南京航空航天大學 能源與動力學院，南京 210016

3.南京航空航天大學 機械結構力學及控制國家重點實驗室，南京 210016

4.南京機電液壓工程研究中心 航空機電系統綜合航空科技重點實驗室，南京 211106

非接觸式動壓氣體軸承具有輕質、低摩阻、無油污染等優勢，作為轉子支撐機構可大幅提高高速旋轉機械的極限轉速和功重比，在制冷透平、沖壓渦輪、冷卻風扇、小型渦輪發動機等航空航天旋轉機械上具有突出的應用優勢［1-4］。以大預緊力裝配和高溫涂層為特征的第四代多葉波箔動壓氣體軸承顯著改善了承載性能和極端工況下的運行穩定性［5］。該類軸承利用轉子與流體間粘性剪切力驅動氣體周向流動，在變截面微間隙內形成動壓壓縮效應而提供承載。動壓氣體軸承轉靜間隙氣膜流動為典型的“泰勒-庫特”流動［6］。由于氣膜流動間隙尺度小（10?6m 量級）、轉子轉速高（105r/min 量級）、剪切效應強，氣膜流動達到滑移流范疇［7］。在稀薄效應和強剪切復合作用下，微尺度間隙內會發生近壁面氣體界面滑移現象。由于界面滑移與環楔形間隙承載氣膜壓力分布強相關，軸承承載等核心性能對其極其敏感，動壓氣體軸承間隙內氣固界面滑移機制及其對軸承性能的影響成為該領域當前的研究熱點。

眾多研究者相繼利用滑移速度模型和極限剪應力模型等對動壓氣體軸承間隙氣膜滑移流動開展研究。學者認為滑移速度模型下界面滑移速度與速度梯度成正比，逐步發展出1.0 階［8］、2.0 階［9］、1.5 階滑移速度模型［10］，并基于上述模型推演出修正的雷諾方程，系統地研究了滑移對動壓氣體軸承間隙內壓力分布和氣膜承載力的影響。譬如，馮凱［11］和黃海［12］等分別基于1.0 階和2.0 階滑移速度模型，分析了軸承間隙氣膜壓力分布以及承載力變化規律。在此基礎上，Wu［13］推演了基于分子動力學理論的滑移速度模型，解決了克努森數Kn>1 時氣膜間隙流動處于分子流狀態下的界面滑移模型問題，后續蘭正義和伍林［14］、燕震雷和伍林［15］利用Wu 滑移模型分析了界面滑移對不同類型動壓氣體軸承性能的影響。極限剪應力模型則認為當流體界面剪應力達到流固界面切應力的極限值時才發生滑移［16］。Spikes和Granick［17］提出了一種極限剪應力為滑移判據復合滑移速度模型的軸承微間隙流動界面滑移模型。馬國軍［18］則通過計算和實驗發現滑移速度模型適合低剪切界面滑移，極限剪應力模型更適合高剪切率的界面滑移。王麗麗等［19］以徑向動壓軸承為研究對象，基于極限剪應力推導了無滑移、轉子側滑移、靜子側滑移、兩側均滑移4 種不同滑移狀態雷諾方程，其分析結果表明了軸承間隙流動處于4 種不同滑移狀態時承載壓力存在明顯差異。李旺［5］、張鏡洋等［20］以2 種類型的波箔型動壓氣體軸承為研究對象，基于極限剪應力模型分析了不同界面滑移狀態對軸承靜特性的影響，結果表明箔片側滑移對靜特性有利，軸面側滑移對靜特性不利。近年來，Bhattacharya［21］、Rao［22］等設計了動壓氣體軸承局部滑移速度梯度分布的氣固界面，發現不同局部滑移分布形態能夠顯著影響軸承性能。實際上，為了形成提供承載的氣膜壓力，實現在負載作用下動壓軸承轉子的偏心轉動，軸承楔形間隙周向流動尺度原本就存在差異，有研究表明其周向克努森數Kn相差可達一個數量級［23］，這會導致近壁面氣體界面滑移狀態存在非一致性。而前述研究中，所涉及的滑移速度模型多認為近壁面氣體周向各局部同步產生滑移；所涉及的極限剪應力模型雖然考慮了轉靜兩側不一致滑移狀態，但是在周向上多依據近壁平均剪應力關系來計算一致的滑移狀態。滑移狀態一致假設或人為設定滑移速度分布，與軸承間隙周向非一致的流動尺度及可壓縮氣體狀態參數相矛盾。尤其在飛行器大過載間隙大偏心下這種矛盾會被加劇，而且第四代多葉波箔型動壓軸承具有多楔型間隙特征、間隙內流動尺度分布更為復雜［5］。因此，工程應用中亟待發展一種適用的動壓氣體軸承間隙周向非一致滑移流動模型，并能掌握其復雜滑移狀態的演化機制。

本文以多葉波箔型動壓氣體徑向軸承為研究對象，推演局部界面剪應力與極限剪應力計算局部滑移速度的方法，建立基于剪應力方程、滑移修正雷諾方程、氣膜厚度方程耦合迭代的流固界面非一致滑移流動分析方法；通過與傳統滑移模型流場分布特征以及實驗值的分析比較驗證模型有效性；進而數值研究旋轉及幾何參數變化下的軸承間隙滑移流動演化規律。以期為飛行器高速旋轉機械非接觸轉子支撐結構設計提供技術支撐。

1 物理模型和數值方法

1.1 物理模型

多葉波箔動壓氣體軸承由軸承套、平箔和波箔組成，波箔作為平箔的支撐元件，一端固定在軸承套上，如圖1（a）所示。其中，O為軸承套圓心，O′ 為軸頸圓心，E0為偏心距，偏心率e=E0/c，c為軸承間隙尺度，ω為旋轉角速度，θ為轉動角，P0為環境壓力。在轉子偏心和裝配后的多楔形箔片結構下流道寬度即氣膜厚度h沿周向呈現明顯變化，并且周向壓力P多峰非均勻分布，這是引起周向克努森數強不均的原因，如圖1（b）所示。軸承結構參數如表1 所示。

表1 軸承模型結構參數Table 1 Structural parameters of bearing model

圖1 軸承結構和間隙參數示意圖Fig.1 Schematic diagram of bearing structure and pa‐rameters in clearance

1.2 滑移速度模型與極限剪應力模型

1）滑移速度模型

滑移速度模型由牛頓內摩擦定律和分子動力學方程導出，適用于描述微尺度下滑移流動邊界條件，基本表達為

式中：us為滑移速度；u為周向線速度；y為徑向位置；λ為平均分子自由程；α為動量協調系數；b為速度偏導二階項系數。

式（1）中關聯了近壁法向速度梯度，反映了剪應力對滑移的影響；平均分子自由程λ、動量調節系數α組成速度梯度項系數，反映微尺度下分子與壁面作用對于界面滑移的影響程度；利用速度偏導二階項系數b修正了速度梯度高階影響程度。為適應不同特征下的滑移流動模擬，學者通常將滑移速度模型分為1.0 階、2.0 階、1.5 階，三者的區別主要在于二階項系數b，分別為0、1/2、2/9。動壓軸承領域學者也分別采用過上述3 種滑移速度模型［8-10］對其間隙流動界面滑移進行模擬，利用滑移速度修正原無滑移邊界速度。軸面側、箔片側近壁氣體線速度u1′、u2′的計算公式分別為

式中：u0為軸頸旋轉線速度。

2）極限剪應力模型

式中：m、n分別為周向和軸向網格數；x為周向位置；Pa為近壁面氣體的平均壓力；τ0為初始剪應力；β為調節系數。

3）傳統滑移模型的問題分析

動壓氣體軸承間隙尺度為μm 級，其間隙流動克努森數Kn范圍為10?5～10?2［15，23］，間隙內局部處于滑移流范疇。從前述動壓氣體軸承界面滑移模型應用研究中可以看出，以往多認為近壁面氣體周向滑移狀態一致，即認為當流固界面出現滑移時則整面各局部均滑移，未細致考慮流固界面會出現局部滑移與局部非滑移區的混合狀態。且滑移速度模型運用流固界面各處相同的動量協調系數α經驗值，極限剪應力模型則通過比較面平均的界面剪應力、極限剪應力，判斷其有無界面滑移或計算滑移速度大小。然而，動壓軸承轉子的偏心轉動和間隙高度的周向非一致，會誘發周向強不勻且非線性的徑向速度梯度分布（見圖2），并且周向不均勻壓力會使得局部可壓縮氣體密度等物性變化較大，表征局部流動尺度的克努森數Kn也會有較大差異。此情形下動壓氣體軸承間隙周向應具有不同的滑移狀態（見圖2）或尺度。但是，在傳統滑移模型中這種影響滑移的局部非一致性因素并未被充分考慮。圖2 中：u1、u2分別為無滑移時軸面側和箔片側線速度；us1、us2分別為軸面側和箔片側滑移速度；下標A、B、C 分別代表圖1（a）中不同截面。

圖2 不同截面界面滑移示意圖Fig.2 Schematic diagram of interface slip in different sections

1.3 界面非一致滑移流動模型

本文考慮近壁面氣體周向局部滑移狀態的非一致性，建立滑移修正雷諾方程，在此基礎上推演局部剪應力計算局部滑移狀態和滑移速度的方法。

無滑移時的無量綱雷諾方程［11-12］為

考慮了軸面與箔片兩側不同的滑移速度，引入無量綱近壁氣體速度對式（9）進行修正，則得到滑移修正的無量綱雷諾方程［7-15］為

式中：±取舍原則為保證τ為正值。

根據牛頓內摩擦定律，軸面側、箔片側界面剪應力τ1、τ2分別為

借鑒極限剪應力經驗關系［16］，獲得關聯局部氣膜壓力的極限剪應力方程

式中：τs為局部極限剪應力；τ0為局部初始剪應力；β為局部調節系數。因兩側材料及壓力相同，故兩側具有相同的τs，根據文獻［18］，取τ0=6 Pa，β=0.000 7。

根據前述界面剪應力和極限剪應力公式，結合式（17），可求出后文軸面側和箔片側局部滑移速度。

根據局部剪應力和極限剪應力方程［24］，可求出軸面側和箔片側局部滑移速度us1、us2，計算公式為

式（18）、式（19）說明若局部近壁氣體流固界面剪應力大于極限剪應力，則判斷局部發生滑移，反之則局部不發生滑移，即滑移速度為0。如此，則利用局部界面剪應力和極限剪應力關聯了局部流動狀態，來計算周向非一致的滑移速度。

1.4 耦合箔片變形的氣膜厚度方程

無量綱氣膜厚度的表達式為

式中：hf為箔片變形量。

箔片變形矩陣根據柔度和壓力矩陣［5］求出

式中：U為柔度矩陣；P為氣膜壓力矩陣。

1.5 網格劃分和邊界條件

對軸承間隙內氣體進行周向和軸向網格劃分，將圓柱面沿周向展開，如圖3 所示。將周向和軸向分別劃分為m、n個單元，形成m+1、n+1個節點；將軸承進口和出口設置為大氣壓；周向展開時軸向兩側為循環邊界，具體邊界條件為

1.6 計算流程和工況

根據前述模型建立滑移修正雷諾方程、局部剪應力方程、彈流潤滑氣膜厚度方程耦合迭代的界面非一致滑移流動分析方法，計算流程如圖4所示。將雷諾方程耦合氣膜厚度方程，對其差分方程采用超松弛迭代法求解，直至氣膜壓力收斂，壓力殘差收斂公式為

圖4 計算流程圖Fig.4 Calculation flow chart

式中：i、j分別為周向和軸向網格點坐標，取值范圍為［0，m］和［0，n］；k為雷諾方程迭代周期數；k′為滑移速度迭代周期數。然后通過兩側剪應力計算兩側局部滑移速度，并結合計算出的兩側局部滑移速度計算滑移速度殘差，軸面側和箔片側滑移速度殘差可分別表示為

判斷軸面側和箔片側局部各節點滑移速度殘差是否收斂，若不收斂則重新計算氣膜厚度和氣膜壓力，經多次循環迭代后可獲得計算域內所有節點的局部滑移速度、滑移狀態、氣膜壓力、氣膜厚度以及承載力，且承載力計算公式為

本文計算工況參數范圍為：軸承數Λ∈[0.15，0.90]、偏心率e∈[0.1，0.8]。

1.7 計算方法驗證

對本文提出的數值計算方法進行驗證。首先，保持軸承參數與文獻［25］中實驗一致，并將計算結果與文獻［25］實驗結果進行對比，如圖5（a）所示。結果表明：不同滑移模型計算的承載力隨轉速變化規律一致，轉速較低時，不同滑移模型計算結果接近，隨轉速升高，不同滑移模型計算結果均高于實驗值。與其他模型相比，非一致滑移模型的計算結果與實驗值最為接近，最大誤差為8.5%。此外，構建軸承間隙壓力分布實驗系統，實驗裝置如圖5（b）所示。實驗測試了軸承間隙c=600 μm、偏心率e=0.6 時軸承Z=L/2 處周向壓力分布，將實驗結果與計算和結果進行對比，如圖5（c）所示。結果表明：n0=6×103，9×103r/min 工況下，數值計算結果與實驗測得的氣膜壓力分布規律一致，且除高壓區和低壓區附近外結果接近，實驗與數值模擬平均表面壓力的誤差在11%以內。與n0=6×103r/min工況相比，轉速提升至n0=9×103r/min 工況，高轉速下軸自身的動力學變形效應，使得轉軸變得不穩定，氣膜厚度處于變化之中，尤其在軸承氣膜壓力高壓區與低壓區位置附近處（位于最小和最大氣膜厚度附近）對氣膜厚度變化極度敏感，這使得理論值與實驗值相差較大。這說明了本文數值計算方法的有效性。

圖5 有效性驗證Fig.5 Validity verification

2 計算結果與分析

2.1 不同滑移模型對壓力和滑移速度分布影響

圖6 為偏心率e=0.6 時不同工況下局部克努森數Kn分布，從圖中可以看出，軸承間隙內局部克努森數Kn沿周向分布變化明顯，最大可達1.8×10?3、最小為4×10?4，相差幾近一個數量級。尤其在周向θ=0.8π～1.5π 位置最小氣膜厚度附近，局部克努森數Kn的范圍為10?3～10?1，處于滑移流區。

圖6 不同工況下局部克努森數Kn 分布（e=0.6）Fig.6 Kn distribution under different working condi‐tions（e=0.6）

圖7、圖8 分別為Λ=0.45，0.90 工況下滑移速度分布。可以看出，在兩工況下，一階滑移速度模型下軸面側和箔片側流固界面均發生滑移，局部滑移速度數值差異較小。極限剪應力模型下Λ=0.45 工況兩側流固界面均不發生滑移；而Λ=0.90 工況軸面側流固界面均發生滑移，箔片側流固界面均無滑移。而兩工況中，非一致滑移模型下兩側局部發生滑移，滑移速度呈現多峰值分布，與氣膜厚度分布狀態相似，局部克努森數Kn較大區域滑移速度較大、較小區域無滑移或者滑移速度較小。本文數值方法獲得的滑移速度分布與局部克努森數Kn分布更符合。

圖7 Λ=0.45 時滑移速度分布Fig.7 Slip velocity distribution at Λ=0.45

圖8 Λ=0.90 時滑移速度分布Fig.8 Slip velocity distribution at Λ=0.90

圖9 為不同滑移模型下軸承Z=L/2 處壓力分布。與無滑移模型相比，一階滑移速度模型下Λ=0.45 工況的高壓區壓力峰值降低、低壓區壓力無明顯變化，最大影響小于1.5%。Λ=0.90 工況的高壓區壓力峰值降低、低壓區壓力峰值升高，但影響較小，最大影響小于5%；極限剪應力模型下Λ=0.45 工況的壓力分布與無滑移模型完全相同。Λ=0.90 工況的高壓區和低壓區壓力峰值均降低，對高壓區影響較大，可達20%，對低壓區影響較小，最大為6%；非一致滑移模型下Λ=0.45 工況的高壓區壓力峰值升高、低壓區壓力峰值降低，最大影響小于2.3%。Λ=0.90 工況的高壓區壓力峰值降低、低壓區壓力峰值升高，最大影響分別為10.6%、5.0%。不同滑移模型與無滑移相比，壓力分布趨勢相同，但對決定軸承承載能力的壓力峰值影響較大，尤其是軸承數越大影響程度越大，相較于無滑移模型分析結果，非一致滑移模型下壓力峰值變化幅度位于滑移速度模型和極限剪應力模型之間。產生上述現象的具體原因：滑移速度模型和極限剪應力模型認為軸承間隙各處流動尺度相當，采用間隙各處相同的滑移模型參數，導致過高或過低估計邊界滑移速度的影響；局部流動狀態與動量協調系數α及面平均的剪應力-τ、-τs強相關，而傳統滑移模型中并未考慮滑移速度與局部流動狀態的耦合影響，在迭代收斂過程中僅進行了單次賦值，因此均出現了滑移速度分布不合理的情況。而本文提出的非一致滑移模型，考慮了偏心間隙尺度及局部壓力分布對剪應力的影響，并將剪應力計算與雷諾方程進行耦合迭代，關聯了局部流動狀態和滑移速度的相互關系。從圖10 中可以看出，經多次迭代后流固界面剪應力變化幅度較大，這證明了剪應力計算與雷諾方程耦合迭代的必要性。經耦合迭代后局部界面剪應力收斂至小于等于局部極限剪應力的合理范圍。因此，關聯滑移與局部流動狀態進行動壓氣體軸承間隙滑移流動分析，可以有效解決周向非一致滑移速度的準確預計問題。

圖9 不同滑移模型下氣膜壓力分布Fig.9 Pressure distribution at different slip models

圖10 壓力和剪應力隨迭代次數變化Fig.10 Variation of pressure and shear stress distribu‐tion with iteration times

2.2 不同軸承參數對界面非一致滑移特征影響

2.2.1 軸承數

軸承數反映了轉靜子相對運動速度和間隙尺度的平方之比的大小，對軸承動壓壓縮效應及周向滑移速度有著顯著的影響。圖11 為偏心率e=0.6 時非一致滑移模型下滑移速度隨軸承數變化規律。隨軸承數增大，箔片側與軸面側依次出現滑移。Λ=0.30 工況下，僅箔片側發生滑移；Λ=0.45 工況下，軸面側也發生滑移。隨軸承數增大，軸面側和箔片側滑移區域沿周向從軸承最小氣膜厚度處逐漸向兩側擴張且滑移速度不斷增大，箔片側局部滑移區域小于軸面側，軸面側滑移區域面積最高約為箔片側的2 倍。圖12為不同工況下Z=L/2 處剪應力分布，從圖中可以看出，兩側剪應力峰值隨軸承數增大而增大，易發生滑移。如圖12（a）所示，Λ=0.30 工況下，僅箔片側部分區域的氣體界面剪應力達到極限剪應力；如圖12（b）所示，Λ=0.60 工況下，軸面側部分區域的氣體界面剪應力也達到極限剪應力，故箔片側和軸面側依次發生滑移。從圖12（c）中可以看出，隨軸承數的進一步增大，兩側氣體界面剪應力達到極限剪應力的區域沿周向從最小氣膜厚度處向兩側逐漸擴張，且剪應力大處滑移速度較大。其中軸面側剪應力更大，故軸面側滑移區域擴張范圍更廣，這是兩側滑移速度分布隨軸承數變化規律的產生內因。

圖11 滑移速度分布隨軸承數變化規律Fig.11 Variation of slip velocity distribution with bearing number

圖12 不同工況下Z=L/2 處剪應力分布Fig.12 Shear stress distribution under different working conditions at Z=L/2

2.2.2 偏心率

偏心率反映了軸承間隙周向尺度變化幅度大小，同樣對周向滑移速度分布有顯著影響。圖13 為軸承數Λ=0.60 時非一致滑移模型下滑移速度隨偏心率變化規律。與軸承數增大時相同，隨偏心率增大，箔片側與軸面側依次出現滑移。e=0.2 工況下，僅箔片側發生滑移；e=0.5 工況下，軸面側也發生滑移。隨偏心率增大，軸面側和箔片側滑移區域沿周向從最小氣膜厚度處逐漸向兩側擴張且滑移速度均不斷增大，箔片側滑移區域小于軸面側。圖14為不同工況下Z=L/2 處剪應力分布，從圖中可以看出，兩側剪應力峰值隨偏心率增大而增大，易發生滑移。如圖14（a）所示，e=0.2 工況下，僅箔片側部分區域的氣體界面剪應力達到極限剪應力；如圖14（b）所示，e=0.5 工況下，軸面側部分區域的氣體界面剪應力也達到極限剪應力。與軸承數增大時不同的是，偏心率增大時兩側滑移區域擴張范圍較小，沿周向遠離最小氣膜厚度處的界面剪應力逐漸降低。其內因是增大偏心率使得間隙周向通道高度變化更為劇烈，偏心率增大僅使得最小氣膜厚度附近徑向速度梯度增大，而沿周向從最小氣膜厚度處向兩側界面剪應力迅速降低，限制了滑移區域的擴張，這不同于軸承數提高整體界面剪應力的情況。

圖13 滑移速度分布隨偏心率變化規律Fig.13 Variation of slip velocity distribution with eccentricity

圖14 不同工況下Z=L/2 處剪應力分布Fig.14 Shear stress distribution under different working conditions at Z=L/2

3 結論

為獲得動壓氣體軸承變截面微間隙局部滑移特征，開展了界面非一致滑移分析方法及滑移速度演化規律研究，得到如下主要結論：

1）針對傳統滑移模型未考慮變截面間隙流固界面局部流動尺度差異問題，推導了局部剪應力與滑移速度的關聯方程，建立了局部剪應力方程、滑移修正雷諾方程、彈流潤滑氣膜厚度方程耦合迭代的流固界面非一致滑移流動分析方法，有效解決了周向非一致滑移速度的預計問題。

2）在本文研究參數范圍內，與傳統滑移模型相比，非一致滑移流動分析方法獲得的滑移速度分布與局部克努森數Kn分布一致性更好。與無滑移相比，不同滑移模型壓力分布趨勢相同，但對決定軸承承載能力的壓力峰值影響較大。一階滑移速度模型和極限剪應力模型對壓力幅值最大影響分別為5%、20%，而非一致滑移模型介于兩者之間，最大影響為10.6%。

3）非一致滑移模型下，隨軸承數和偏心率的增大，箔片側與軸面側依次出現滑移，且兩側滑移區域沿周向從軸承最小氣膜厚度處逐漸向兩側擴張，軸面側滑移區域面積最大為箔片側的2倍。但軸承數增大時滑移區域面積迅速增大，而偏心率增大時引起局部流動尺度差異變大，僅最小氣膜厚度處的滑移面積增大。