波箔動壓氣體軸承建模中典型摩擦模型的對比分析

周如能，顧永鵬，任革學，于溯源

(清華大學 a.能源與動力工程系；b.核能與新能源技術研究院；c.航天航空學院，北京 100084)

箔片動壓氣體軸承具有無油、高溫和高速性能好等優點，被廣泛應用于飛機的空氣循環機和燃料電池的壓縮機等領域[1-3]，其中波箔式箔片軸承應用最為廣泛[4-5]。由于波箔結構的特殊性，波箔動壓氣體軸承系統一部分非線性來源于箔片結構的接觸摩擦行為[6]，系統動力學建模中箔片結構的接觸摩擦模型是準確預測波箔動壓氣體軸承靜、動態特性的關鍵[7]。

文獻[8]最早提出了波箔氣體軸承的理論模型——簡單彈性基礎模型，后來被拓展到動力學計算中，并引入了黏性阻尼近似模擬結構耗散特性[9-10]，但結構耗散特性與軸承-轉子系統的運動狀態關系較大，根據某一工況估算的等效黏性阻尼不能準確表征實際的能量耗散機制。這一局限性使很多學者提出了考慮波箔結構及波箔之間相互作用的模型耦合摩擦模型，進而對波箔動壓氣體軸承系統進行靜、動力學分析[11-13]。

摩擦模型較多，其選擇受系統求解方法的限制。傳統的顯式分別迭代傳參求解方法可以采用最經典的庫侖摩擦模型，在每個時間步內迭代判斷每個摩擦節點的狀態并得到摩擦力[12]，但各子系統變量之間存在時間滯后，需極小的時間步長來保證收斂性，計算效率低，故文獻[14-15]提出了全耦合求解方法，系統所有變量可同時求解，大幅提高了計算效率。但由于庫侖摩擦模型在零速度時的不連續性而無法直接使用，需要改進摩擦模型。按是否采用微分方程庫侖摩擦模型的改進可分為無需微分方程的代數摩擦模型和采用微分方程的動態摩擦模型[16]。最典型的代數摩擦模型是在經典庫侖摩擦的間斷處引入近似函數進行平滑處理，文獻[17]采用代數摩擦模型實現了波箔動壓氣體軸承的全耦合建模與求解，但平滑模型無法準確描述靜摩擦，即無法準確預測箔片節點的黏著狀態，進而影響軸承性能的預測。波箔動壓氣體軸承分析中使用的動態摩擦模型主要包括Petrov和Ewins提出的動態摩擦模型[18]和LuGre動態摩擦模型[19-20]，前者同樣是基于符號函數的近似衍生得到，理論上無法準確描述靜摩擦，后者可以很好地描述斯特里貝克效應以及黏著-滑移現象。文獻[19]最早將LuGre動態摩擦模型引入波箔動壓氣體軸承，并與link-spring模型進行聯合仿真，但其分析局限于箔片內部的力學行為，未考慮轉子和氣膜的行為，且求解方法仍為分別迭代傳參求解。

本文在上述研究的基礎上，將LuGre動態摩擦模型引入全耦合的波箔動壓氣體軸承模型，提出了考慮轉子、動壓氣膜、箔片結構、波箔之間的相互作用以及箔片節點之間黏著-滑移接觸摩擦行為的綜合動力學模型[20]。基于該模型，進一步對比了LuGre動態摩擦模型和平滑近似代數摩擦模型在程序實現、穩態特性以及軸承性能預測方面的差異。

1 波箔動壓氣體軸承-轉子系統動力學建模及求解

波箔動壓氣體軸承結構如圖1所示，以軸承中心O為原點建立坐標系Oxyz，θ為位置角(與x軸負方向的夾角)。工作時轉子高速旋轉，在轉子與頂箔之間形成一層流體動壓氣膜，流體動壓力作用于轉子上使轉子懸浮的同時也作用于頂箔上，而頂箔和波箔的受力變形會導致波箔與頂箔及軸承座之間出現復雜的接觸摩擦行為。

1—氣膜；2—焊接端；3—頂箔；4—波箔；5—軸承座。

為了求解這一包含可壓縮流體潤滑、接觸摩擦等復雜非線性行為的流固耦合問題，氣膜壓力和箔片變形均采用標準的有限單元法進行空間離散，如圖2所示。氣膜壓力沿氣膜厚度方向的變化可忽略，不同位置的氣膜壓力可表示為p(θ,z)。對于第一代波箔動壓氣體軸承，一般假設箔片的曲率效應及變形沿軸承軸向的變化可忽略，為分析變形，以軸承周向為xf軸，徑向為yf軸建立平面坐標系。為提高計算效率，箔片采用二維的梁單元建模。在標準的有限元框架下，氣膜壓力p(θ,z)和箔片變形u可以分別用單元節點向量插值得到，即

(1)

式中：Np，Nq分別為氣膜壓力和箔片變形的形函數矩陣；pe和qe分別為單元節點的壓力向量和變形向量。氣膜壓力和箔片變形的廣義坐標向量分別標記為p和q。

1.1 轉子的控制方程

因2套軸承對稱支承一個剛性轉子時，2套軸承位置轉子響應相同，在此僅建立一套軸承模型，轉子等效為質點，記xr為轉子的位移矢量， 則轉子

圖2 有限單元法離散示意圖

的控制方程為

(2)

式中：Mr為轉子質量矩陣；Wr為軸承靜載荷，本文指轉子重力；Fr為氣膜壓力對轉子的合力；Fub為轉子不平衡力。

Fr可通過氣膜壓力p(θ,z)積分得到，即

z∈[-L/2,L/2]，

(3)

Aθ=[cosθ,sinθ]T，

式中：R為軸承半徑；L為軸承長度；pa為環境壓力；Ωe為單元面積。

不平衡力矢量Fub為

(4)

式中：u為不平衡量；ω為轉子角速度；θ0為不平衡力的初始相位。

1.2 氣膜的控制方程

氣膜壓力p(θ,z)的控制方程為等溫可壓縮理想流體的雷諾方程，即

(5)

U=[ωR/2,0]T，

式中：h為氣膜厚度；μ為空氣動力黏度。

忽略轉子和軸承的熱膨脹以及離心膨脹，假設軸承間隙C為常數，氣膜厚度h可表示為

h=C(1+excosθ+eysinθ)+hf，

(6)

[ex,ey]T=xr/C，

hf=-Nq,yfqe，

式中：ex，ey分別為在x，y方向量綱一的轉子位移；hf為箔片變形導致氣膜厚度的變化量；Nq,yf為頂箔單元對應于yf方向變形的形函數矩陣。

由于氣膜與外界環境相通，氣膜控制方程的邊界條件為在所有邊界上p=pa。

氣膜壓力采用雙線性四節點單元進行劃分，利用布勃諾夫-伽遼金法得到單元雷諾方程的弱形式為

(7)

Bp=?Np，

氣膜的整體控制方程可以通過組裝所有單元方程得到

(8)

1.3 箔片的控制方程

采用二維歐拉梁單元建立箔片模型，單元的控制方程為

(9)

Ce=βKe，

式中：Me，Ce，Ke分別為單元的質量矩陣、阻尼矩陣(此處采用比例阻尼)和剛度矩陣；Fe為單元的外部力矢量；β為比例阻尼系數。此外，箔片材料彈性模量需根據平面應力假設修正。

由于波箔大部分節點不受力，且波箔動力學特性可以忽略，模型引入Guyan縮減法縮減波箔的自由度[20]，在保證計算精度的同時可大幅提高計算效率，且縮減后的波箔子結構可以視為超級單元。

組裝所有的頂箔單元方程及縮減后的波箔子結構方程得到箔片整體的控制方程為

(10)

對于平滑近似代數摩擦模型，(10)式可簡化為

(11)

(10)，(11)式的邊界條件為軸承尾端的波箔和頂箔均焊接固定于軸承座。

1.4 接觸摩擦的控制方程

根據波箔縮減策略，采用點對二維梁單元的接觸模型來表征箔片結構的接觸摩擦行為。接觸力包括法向接觸力fn和切向摩擦力ft，基于赫茲接觸模型得到法向接觸力fn為

(12)

式中：k為法向接觸剛度系數；δ為接觸變形；c為法向接觸阻尼系數。

切向摩擦力ft分別采用LuGre動態摩擦模型和平滑近似代數摩擦模型表征。

1.4.1 LuGre動態摩擦模型

LuGre動態摩擦模型可以認為是Dahl模型的改進模型，可以表征斯特里貝克效應以及相應的黏著-滑移運動[21]。對于單個節點摩擦，LuGre動態摩擦模型的控制方程為

(13)

(14)

式中：σ0，σ1，σ2分別為摩擦剛度系數、微觀阻尼系數和黏性阻尼系數；zt為接觸點的內部摩擦狀態變量，描述了接觸表面上的微觀變形；vt為兩接觸體的相對速度；vc為轉換速度，應盡可能小，以保證摩擦模型具有較好的微觀阻尼[22]；g(vt)為表征庫侖摩擦和斯特里貝克效應的函數。

g(vt)一般可表示為

g(vt)=fsl+(fst-fsl)e-(vt/vs)2，

(15)

fsl=μffn，

式中：fsl，fst分別為滑動摩擦力和靜摩擦力；vs為斯特里貝克速度；μf為摩擦因數。

定義黏著系數γ=fst/fsl表征斯特里貝克效應。當滑動速度為常數時，穩態切向摩擦力ft,s可以表示為

ft,s=g(vt)sgn(vt)+σ2vt，

(16)

這與大多數代數摩擦模型一致。

將所有摩擦節點的狀態變量定義為矢量z，得到LuGre動態摩擦的控制方程為

(17)

1.4.2 平滑近似代數摩擦模型

采用最常用的雙曲正切函數近似模型[17]，切向摩擦力為

ft=μffntanh(rvt)，

(18)

式中：r為平滑系數，取值越大越接近經典庫侖摩擦模型，但也將導致嚴重的數值困難，平滑系數的取值直接決定了模型的計算精度和效率。

1.5 系統耦合方程的建立及求解

當采用LuGre動態摩擦模型時，聯立轉子、氣膜、箔片和接觸摩擦的控制方程，得到系統的控制方程為

(19)

定義系統狀態變量為

(20)

將(19)式整理為一般形式為

(21)

當采用平滑近似代數摩擦模型時，系統的控制方程為

(22)

定義系統狀態變量為

(23)

同樣可將(22)式整理為(21)式的一般形式。

基于MATLAB的變階變步長的隱式積分器ode15s[23]對(21)式的所有狀態變量同時積分求解，隱式積分器需要獲得系統方程的雅克比矩陣，仿真程序中采用了解析的雅克比矩陣?F/?y提高計算效率。此外，本模型還忽略了頂箔與波箔之間的分離行為，并引入Gümbel條件修正，即在(3)式中忽略p

2 2種摩擦模型對比分析

2.1 程序實現

在氣體軸承幾何發散區會存在低壓區，該區域頂箔與波箔之間的法向接觸力趨近于0。對于LuGre動態摩擦模型，法向接觸力過小會使(14)式中的g(vt)趨近于0，進而導致數值困難。需定義一個極小的正實數作為接觸容差，δ<的節點對應的摩擦狀態自由度和方程需要被刪除掉，并在δ≥時添加回系統方程。在仿真中還需要執行一個檢測函數，當檢測到接觸狀態變化時終止仿真程序，并在修改系統方程自由度后從斷點繼續仿真。此外，LuGre控制方程對系統雅克比矩陣的貢獻項也需要進行較為繁瑣的推導與編程。故LuGre動態摩擦模型會增加系統自由度，且程序實現復雜。

而平滑近似代數摩擦模型既不會增加系統自由度，也無需在仿真中調整系統方程，程序實現簡單。

2.2 穩態特性

2種摩擦模型的切向摩擦力隨兩接觸體的相對速度變化如圖3所示，LuGre動態摩擦模型未考慮黏性阻尼項，且摩擦力ft基于滑動摩擦力μffn進行了歸一化。由圖3可知：LuGre動態摩擦模型能完全表征靜摩擦(節點黏著狀態)，平滑近似代數摩擦模型理論上靜態時的摩擦力為0，僅能依靠較低的滑移速度來近似表征靜摩擦。

圖3 2種摩擦模型的切向摩擦力隨兩接觸體的相對速度的變化趨勢

2.3 性能預測

以某波箔動壓氣體軸承為研究對象，其主要結構參數見表1。箔片材料參數為：密度7 800 kg/m3，彈性模量207 GPa，泊松比0.29。運行條件：pa=101 325 Pa，μ=1.95×10-5Pa· s。LuGre動態摩擦模型基本參數見表2。氣膜周向和軸向的單元數分別為50和20，頂箔被劃分為100個梁單元。此外，k=1×106N/mm，c=1×102N· s/mm，=1×10-12mm，β=1×10-4s。ode15s函數的相對容差設為1×10-6。為方便對比，黏性阻尼和黏著系數的取值使LuGre動態摩擦模型的穩態摩擦力和庫侖摩擦完全一致。代數摩擦模型中的平滑系數r取10 s/mm[17]，這一取值已經導致了嚴重的數值困難，平均計算時間約為LuGre摩擦模型的5倍，且摩擦因數越大計算效率越低。系統建模及求解方法的準確性已在文獻[20]中得到驗證。

表1 波箔動壓氣體軸承主要結構參數

表2 LuGre動態摩擦模型基本參數

2.3.1 轉子靜態性能

本節所研究的靜態性能為在不同靜載荷下轉子的靜態偏心率，轉子轉速為30 000 r/min，摩擦因數為0.1，主要任務是求解(21)式的靜平衡方程F(y0)=0得到其靜態解y0。由于箔片的摩擦節點較多，導致系統在同一工況下存在多個靜態平衡點，且無法直接將方程在其平衡點附近進行線化迭代求解；此外，平滑近似代數摩擦模型直接求靜態解無意義，因為在理論上該模型的靜態摩擦力為0：故計算的平衡點為通過時域仿真在固定轉速下從坐標原點準靜態逐步加載得到。

2種摩擦模型轉子的靜態偏心率隨靜載荷的變化曲線如圖4所示(為方便對比，圖中還給出了無摩擦的計算結果， 即平滑近似代數摩擦模型的理論結果)：1)隨靜載荷增加，2種摩擦模型的轉子靜態偏心率差異逐漸增大；2)緩慢加載情況下的平滑近似代數摩擦模型處于LuGre動態摩擦模型和無摩擦模型之間，即在非零的低速滑移情況下平滑近似代數摩擦模型可以在一定程度上表現出接近黏著特性的行為。

圖4 2種摩擦模型的轉子靜態偏心率隨靜載荷的變化曲線

2.3.2 轉子穩定性

所研究的穩定性包括不同摩擦因數下系統的局部穩定性(系統在平衡位置處受小擾動后回到靜態平衡位置穩定運行的能力)和全局穩定性(系統受任意非破壞性的沖擊后回到靜態平衡位置穩定運行的能力)[24]。轉子質量為10.2 kg，靜載荷為100 N，不平衡量為0。

局部穩定性和全局穩定性的關系取決于Hopf分岔行為(系統平衡點局部失穩時，根據系統參數的不同所形成的空間分岔形式包括超臨界Hopf分岔和亞臨界Hopf分岔)，如圖5所示：1)對于超臨界Hopf分岔，在穩定極限環存在的轉速區域內，任意初始條件的轉子軌跡都將收斂于該穩定極限環， 系統的全局穩定性和局部穩定性一致； 2)對于亞臨界Hopf分岔， 在不穩定極限環存在的轉速區域內， 擾動較小的轉子軌跡將收斂于穩定的平衡點，而擾動較大的轉子軌跡將向外發散或收斂于外側可能存在的穩定極限環，系統局部穩定，全局不穩定。

(a)超臨界Hopf分岔

局部穩定性的計算需要確定轉子的平衡位置并施加小擾動。平衡位置通過時域仿真在固定轉速下從坐標原點準靜態逐步加載得到。隨后，給予轉子1‰的擾動并仿真足夠長的時間來判斷系統的局部穩定性。若轉子軌跡收斂，系統在該轉速下局部穩定；反之，則為局部不穩定。隨轉子轉速升高，系統在某一轉速下徹底失去局部穩定性，這一閾值轉速即為局部失穩轉速。

全局穩定性的計算主要關心的是受各種非破壞性沖擊后的響應，沖擊效果可以等效為轉子位置主導的初始條件。對于一般的非破壞性沖擊，轉子初始位置基本都在軸承間隙圓內。轉子初始時在固定轉速下被約束到間隙圓內的不同位置，并根據此時軸承力、靜載荷和約束力達到靜平衡的條件，計算得到初始的氣膜壓力、箔片變形以及摩擦狀態。隨后，移除約束力釋放轉子并仿真足夠長的時間來判斷系統的全局穩定性。若從所有初始位置釋放的轉子都回到了靜態平衡位置，系統在該轉速下全局穩定；反之，則為全局不穩定。類似地，隨轉子轉速升高，系統在某一轉速下徹底失去全局穩定性，這一閾值轉速即為全局失穩轉速。

2種摩擦模型局部失穩轉速隨摩擦因數的變化曲線如圖6所示：1)在低摩擦因數區，2種模型局部失穩轉速一致，這是因為箔片摩擦節點在低摩擦因數下主要處于滑移狀態，2種摩擦模型均能很好地表征動力學行為；2)隨摩擦因數增大，平滑近似代數摩擦模型預測的局部失穩轉速先大于LuGre動態摩擦模型， 然后小于LuGre動態摩擦模型，這是由于隨摩擦因數增大， 處于黏著狀態的摩擦節點增多，平滑近似代數摩擦模型無法準確描述黏著狀態。上述分析說明平滑近似代數摩擦模型在低摩擦因數區可準確預估局部失穩轉速，在中摩擦因數區會高估局部失穩轉速，在高摩擦因數區會低估局部失穩轉速。

圖6 2種摩擦模型局部失穩轉速隨摩擦因數的變化曲線

2種摩擦模型全局失穩轉速隨摩擦因數的變化曲線如圖7所示：1)在低摩擦因數區，2種模型全局失穩轉速一致，這是由于幾乎所有箔片摩擦節點均處于滑移狀態；2)在中摩擦因數區，2種模型全局失穩轉速出現差異，但沒有局部失穩轉速大，這是因為與局部的小擾動相比，全局穩定性計算中轉子的擾動幅度較大，箔片摩擦節點大多處于滑移狀態；3)在高摩擦因數區，即使是大幅擾動下很多節點也主要處于黏著狀態，平滑近似代數摩擦模型的全局失穩轉速明顯低于LuGre動態摩擦模型。

圖7 2種摩擦模型的全局失穩轉速隨摩擦因數的變化曲線

3 結論

對比分析了波箔動壓氣體軸承建模中LuGre動態摩擦模型和平滑近似代數摩擦模型在程序實現、穩態特性以及性能預測方面的差異，得到以下結論：

1)LuGre動態摩擦模型可以準確表征黏著狀態，在任何仿真工況下預測結果均比較準確，且計算效率較高；但該模型會增加系統自由度，且需要對零接觸力進行特殊處理并執行檢測函數，導致其程序實現較復雜。

2)平滑近似代數摩擦模型無需增加系統自由度，程序實現簡單，且在低速準靜態或低摩擦因數、大擾動等以滑移為主的工況下預測結果比較準確；但該模型無法準確表征黏著狀態，在高摩擦因數、小擾動等以黏著為主的工況下預測結果偏差較大，且計算效率較低。