計算流體潤滑頻變動力學特性的一種完整變量擾動偏導數法

畢春曉， 韓東江， 楊金福

(1. 中國科學院 工程熱物理研究所，北京 100190；2. 中國科學院大學，北京 100049)

高速旋轉機械應用廣泛，在能源動力、國防軍工、航空航天、車輛、艦船等具有高性能、極端參數、長壽命要求的領域發揮著重要作用[1]。高效率和結構緊湊的特點使旋轉機械在保證可靠性的同時向高轉速發展，因此高速軸承和密封的性能分析、優化設計及其軸系穩定性控制已成為關鍵問題[2]。

使用工質流體作為潤滑介質的軸承既能徹底避免工質被油污染、降低密封技術難度，同時縮短轉子以提高渦輪機軸系的動力學性能[3]，是國內外相關領域的重要發展方向之一[4]，比如自潤滑軸承用于低溫工程[5]和有機朗肯循環(organic rankine cycle, ORC)的工質泵[6]，液體火箭發動機的高速渦輪泵[7-8]。1970年美國的空間Brayton循環項目中的發電裝置最終采用He-Xe工質潤滑箔片軸承[9]。超臨界二氧化碳(S-CO2)閉式Brayton循環在太陽能光熱發電[10]、船舶和飛行器余熱利用[11-12]等領域有著廣闊前景。目前S-CO2自潤滑軸承已用于多個高速渦輪機樣機[13-14]，其結構形式的研究和設計已引起國內外越來越多學者的重視[15-17]。

特殊參數下的流體通常具有復雜的熱物性，從而帶來潤滑膜流動的多種附加效應，如表1所列舉。具體到數學模型即潤滑雷諾方程中不只含有膜厚和壓力兩個顯函數變量，還包含熱物性變量以及表示附加效應的非線性函數。

表1 特種參數流體潤滑的附加效應Tab.1 Additional effects of fluid lubrication with special parameters

上述文獻對于軸承或密封的潤滑流動附加效應研究集中于穩態特性，然而這些復雜效應更加導致動力學特性的求解難度明顯增大。Lund[26]提出計算滑動軸承剛度阻尼系數的攝動法并結合“瓦塊裝配法”求解可傾瓦軸承同步動力學系數，即軸頸和瓦塊的擾動頻率等于軸頸旋轉的圓頻率。雖然該方法在推導過程中通過歸并消掉了位移和速度擾動量，從而天然滿足臨界特性與線性失穩分析所要求的平衡位置鄰域內擾動[27]。然而大量試驗結果表明，無論是油潤滑可傾瓦軸承[28]還是氣體軸承[29-30]，在相同的旋轉頻率(轉速)下，動力學剛度和阻尼系數隨擾動頻率變化。文獻[31]通過含擾動頻率的復指數描述小擾動下的膜厚和壓力動態變化，從而能夠處理空氣軸承雷諾方程中的時滯壓力項，并揭示了可壓縮潤滑軸承動力學特性的頻率效應。而且頻率擾動方法能夠納入軸瓦與彈性支承的運動和動態變形引起的頻率效應，實現箔片軸承或可傾瓦軸承的動態耦合分析[32-33]。比如Yang等[34]將偏導數法用于可傾瓦空氣軸承，求解包含瓦塊和支點擾動的全局動力學系數，并在此基礎上提出一種穩定性主動控制方法。李長林[35]將接觸約束引入結構擾動模型，分析了多滑動梁空氣箔片軸承的線性穩定性，其中的動力學系數反映了箔片結構之間庫倫摩擦阻尼起到的作用。

然而對于真實氣體，密度-壓力呈非線性關系，沿用空氣軸承的處理方式(將狀態方程代入雷諾方程消去熱物性變量)將會使方程變得非常復雜，用來求解頻變動力學特性更加困難。而且對于熱力學行為嚴重偏離理想氣體的超臨界流體，需要采用基于數值計算的模型或數據庫才能得到較為準確的熱物性結果[36]。溫建全基于小擾動法分析了S-CO2箔片軸承的動特性，并與位移速度增量法進行對比，由于在推導擾動方程時只考慮了壓力擾動，從而文中對兩種方法結果差異的解釋為：熱物性和湍流系數也會隨著轉子擾動發生動態變化[37]。Bi等[38]對偏導數法做出改進，通過建立動態映射引入包含擾動密度、黏度、擾動雷諾數和湍流系數的完整變量擾動，首次給出任意擾動頻率下考慮真實氣體和湍流動態影響的超臨界二氧化碳潤滑軸承動力學剛度阻尼系數的求解方法。江錦波等[39]進一步拓展了Bi等的方法，將動態離心慣性效應引入完整變量擾動，用于超臨界二氧化碳干氣密封的動力學特性研究。

本文給出考慮完整變量擾動的偏導數法，推導可壓縮湍流潤滑的擾動雷諾方程，介紹動力學系數求解方法，適用于具有附加效應的軸承與密封的頻變動力學特性計算?；跍熟o態處理的可壓縮湍流潤滑雷諾方程，通過對比偏導數法和有限增量法求解的S-CO2圓柱瓦軸承剛度阻尼系數，驗證了本文的偏導數法可確保動力學系數包含附加效應動態變化的影響。同時澄清基于頻率擾動的和有限位移速度增量的兩類方法各自的特點，并以S-CO2和油潤滑高速可傾瓦軸承為案例展現本文方法的應用。

1 考慮完整變量動態變化的偏導數法

1.1 瞬態雷諾方程的頻率擾動分析

為了使偏位角迭代易于處理，本文采用始于豎直位置的潤滑膜圓周方向角度θ，如圖1所示。一般形式的可壓縮湍流潤滑雷諾方程為

(1)

式中：R為軸承半徑；θ為環向角度坐標；z為軸向坐標；ρ為密度；μ為動力黏度；h為潤滑膜厚；p為潤滑膜壓力；ω為旋轉圓頻率；t為時間；kθ和kz分別為圓周方向和軸向湍流潤滑系數，可通過局部雷諾數Re的非線性函數描述湍流效應。

kθ=12+α1Reβ1=12+0.013 6Re0.9,
kz=12+α2Reβ2=12+0.004 3Re0.96

(2)

式中： 常數α1，β1，α2，β2的值由Ng-Pan湍流潤滑模型確定，該模型在來流馬赫數小于5的情形對可壓縮和不可壓縮潤滑均適用[40]；層流時上述所有常數均等于零，由于潤滑膜同時存在壓差流和剪切流，以Re<1 000為層流判據。

式(1)不依賴密度和黏度模型，因而不局限于特定的潤滑流體，基于該方程可發展出通用的潤滑特性分析方法。

(3)

式中：ψ=C/R為間隙比；C為名義半徑間隙；Λ=μaω/2paψ2為軸承數，是可壓縮潤滑問題特有的無量綱參數；μa為環境黏度；ρa為環境密度；pa為環境壓力。

(4)

在圖1中，ε和φ為軸頸在擾動狀態下某一瞬時所在位置的偏心率和偏位角； 對應的ε0和φ0為軸頸穩態平衡位置的偏心率和偏位角；Ob為軸承中心；Oj和Oj0為瞬時位置和穩態平衡位置的軸頸中心。

在擾動狀態下，可壓縮湍流雷諾方程中的密度、黏度和湍流系數同樣都發生與壓力擾動和膜厚擾動耦合的動態變化。本文給出雷諾方程中完整變量的擾動形式并建立與壓力和膜厚擾動的動態映射關系，進而可將基于數據庫和數值方法的熱物性模型納入該求解體系。

(5)

(6)

局部雷諾數Re的動態變化反映了湍流對動力學特性的影響，引入由名義半徑間隙、環境密度和環境黏度定義的平均雷諾數Rea=ρaωRC/μa，可得由無量綱變量表示的動態局部雷諾數。

(7)

由冪函數的泰勒公式，結合式(6)有

(8)

將式(4)、式(6)和式(8)代入式(7)，略去二階及以上的擾動項，得到動態局部雷諾數Re的擾動形式。

(9)

式中，Re0為對應于軸頸平衡位置的穩態雷諾數，擾動雷諾數Red是擾動密度、擾動黏度和擾動膜厚的聯合響應。

利用冪函數的泰勒公式不難推出Reβ。

(10)

(11)

1.2 穩態和擾動雷諾方程

(12)

(13)

1.3 動力學系數的求取方法

(14)

式中：L為軸承長；D為軸承直徑。

式(15)可將動力學系數從耦合坐標系變換到x-y直角坐標系。

(15)

式中， [A]為旋轉矩陣，其表達式為

(16)

數值積分公式和旋轉矩陣的具體表達式與坐標軸指向和軸頸旋轉方向有關，式(14)和式(16)對應圖1的坐標系。對于可傾瓦軸承，擾動膜厚表達式中還包含瓦塊擺動和支點運動等相關的擾動變量，并需要借助瓦塊-支點動力學方程推導折合剛度和阻尼系數，具體過程同Yang等的研究。

2 完整變量擾動偏導數法的驗證

2.1 有限增量法(僅用于驗證)

(17)

(18)

求解各方向位移和速度增量對應的五個類似穩態雷諾方程形式的式(17)，通過增量壓力分布與穩態壓力分布的差值可提取剛度阻尼系數。

(19)

其中，

(20)

動力學系數從ε-φ坐標系轉換到圖1的x-y坐標系公式如式(21)所示，其中[A]與式(16)相同。

(21)

這類方法的命名并不統一，表達形式也不盡相同，文獻[41]以載荷增量方式實現剛度阻尼系數的提取，方法名稱中文直譯為數值微分法；國內大部分資料[42]稱為有限差分法，是通過在不同位置和速度下積分得到油膜力的差值實現的。溫建全提出一種新的實現方式，并根據方法的特點稱為位移速度增量法。本質上都是通過給定有限的軸頸位移和速度實現剛度阻尼系數求解，因而為了避免與數值方法中的差分或微分的概念造成混亂，本文稱之為有限增量法，制造增量的方式可借鑒溫建全的研究。

對于可壓縮潤滑來講，該方法得到的是靜態剛度和阻尼系數，即擾動頻率無限趨近于零Ω→0+的極端情況。但是能夠在準靜態處理(去掉可壓縮時滯項?ρ/?t)的前提下，自動包含密度、黏度和湍流系數的動態變化效應，因此可用來驗證考慮完整變量動態變化的偏導數法。

2.2 對比驗證結果

表2 S-CO2 圓柱瓦軸承的結構和運行參數Tab.2 Parameters of S-CO2 cylindrical bearing

3 案例計算結果與討論

3.1 案例1—超臨界二氧化碳高速可傾瓦軸承

表3 S-CO2高速可傾瓦軸承的參數Tab.3 Parameters of S-CO2 tilting pad bearing

需要說明的是，本章兩節內容中可傾瓦軸承的偏心率指的是以裝配半徑間隙Cb為參照的裝配偏心率，記作ε′0，以區分潤滑膜厚度公式中的偏心率ε0，后者以名義半徑間隙(瓦塊半徑間隙)Cp為參照。兩個偏心率的關系ε0=ε′0(1-m)為，其中m=1-Cb/Cp為可傾瓦軸承的預負荷系數。本章兩個案例的結果與討論中所提及的“偏心率”指裝配偏心率。

圖5顯示了偏心率0.1，轉速55 000 r/min的可傾瓦軸承動力學剛度和阻尼系數隨著無量綱擾動頻率Ω的變化。在Ω較低時，準靜態假設的剛度系數結果接近完整變量擾動，而只考慮壓力擾動導致剛度結果明顯偏大。隨著Ω增加，剛度系數明顯增大，在Ω<1.5范圍內只有壓力擾動的結果均偏大。

隨著的增加，阻尼系數先緩慢減小，大約從Ω>0.5開始明顯減小，只考慮壓力擾動或準靜態假設都導致阻尼系數結果明顯偏大。

圖6為偏心率0.8，轉速20 000 r/min的結果，剛度系數的變化趨勢與圖5高轉速小偏心情形基本一致，只是增長較緩慢。阻尼系數隨著Ω的增加而緩慢增大。完整變量擾動對阻尼系數結果的影響規律與圖5相同。

對于可傾瓦軸承，可壓縮時滯效應增大了剛度系數而減小了阻尼系數。忽略熱物性和湍流的動態變化效應將導致剛度和阻尼系數明顯偏大。

3.2 案例2—高速可傾瓦油潤滑軸承

由于不可壓縮潤滑的雷諾方程中不存在壓力或密度的時滯項，油潤滑可傾瓦軸承的頻率效應與潤滑膜自身特性無關。因而也有一類方法采用Lund攝動法處理油膜擾動并與引入頻率的瓦塊和支點動力學模型結合，可求解不同頻率的可傾瓦軸承動力學系數[43]。本文的偏導數法能夠統一處理潤滑膜和結構的頻率擾動。首先將式(3)各項中表征可壓縮性的密度變量去掉，然后用偏導數法計算高速可傾瓦油潤滑軸承的折合剛度和阻尼系數，求解方程時采用雷諾邊界條件處理油膜破裂。軸承的結構形式和坐標系選自手冊《Journal Bearing Data Book》[44]第二部分No.49數據的五瓦瓦上加載可傾瓦動壓軸承。Someya采用Lund攝動法計算可傾瓦軸承同步動力學系數，因此在本節的偏導數法計算中，無量綱擾動頻率Ω=1。軸承的結構和運行參數列于表4。

表4 高速可傾瓦油潤滑軸承的參數Tab.4 Parameters of tilting pad oil bearing

圖8顯示了可傾瓦軸承的無量綱折合剛度和阻尼系數隨偏心率的變化，圖8中縱坐標為無量綱動力學系數與無量綱承載力的比值，與手冊《Journal Bearing Data Book》的無量綱方式一致。當偏心率大于0.6時，對于x方向的剛度阻尼，偏導數法的層流結果比手冊中數據偏大，這是由于本節的計算沒有采用質量守恒邊界條件處理油膜破裂。

對比考慮湍流系數擾動和只有壓力擾動的湍流結果，由于擾動項不影響靜特性因而兩種求解條件下無量綱承載力相同，可以看出忽略湍流系數的動態變化效應將使剛度系數結果偏大，偏心率越小差異越明顯；而湍流系數擾動對阻尼系數幾乎沒有影響。

4 結 論

本文給出考慮完整變量擾動的偏導數法，從一般形式的可壓縮湍流潤滑雷諾方程出發，推導擾動變量的復合動態映射表達式以及擾動雷諾方程，并求解軸承動力學剛度與阻尼系數。該方法適用于具有真實氣體、湍流等附加效應的軸承與密封流體膜頻變動力學特性計算。對方法進行驗證并給出應用案例，計算分析的結果可得如下結論。

(1) 在準靜態假設下(可壓縮湍流潤滑雷諾方程去掉密度時滯項)用偏導數法求解S-CO2圓柱瓦軸承的剛度與阻尼系數，與有限增量法結果對比，驗證了本文所提出的偏導數法求解動力學系數的準確性與有效性。

(2) 本文的偏導數法能夠統一處理結構擾動和潤滑膜可壓縮性共同作用產生的動力學系數的頻率效應。S-CO2高速可傾瓦軸承的分析結果表明，忽略壓力和膜厚以外的擾動變量將對剛度與阻尼系數結果帶來很大偏差。

(3) 高速可傾瓦油潤滑軸承在潤滑膜湍流效應明顯時，若忽略湍流系數的動態變化效應將導致小偏心率下的剛度系數結果偏大。