計算流體潤滑頻變動力學(xué)特性的一種完整變量擾動偏導(dǎo)數(shù)法

畢春曉， 韓東江， 楊金福

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

高速旋轉(zhuǎn)機(jī)械應(yīng)用廣泛，在能源動力、國防軍工、航空航天、車輛、艦船等具有高性能、極端參數(shù)、長壽命要求的領(lǐng)域發(fā)揮著重要作用[1]。高效率和結(jié)構(gòu)緊湊的特點使旋轉(zhuǎn)機(jī)械在保證可靠性的同時向高轉(zhuǎn)速發(fā)展，因此高速軸承和密封的性能分析、優(yōu)化設(shè)計及其軸系穩(wěn)定性控制已成為關(guān)鍵問題[2]。

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

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

表1 特種參數(shù)流體潤滑的附加效應(yīng)Tab.1 Additional effects of fluid lubrication with special parameters

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

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

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

1 考慮完整變量動態(tài)變化的偏導(dǎo)數(shù)法

1.1 瞬態(tài)雷諾方程的頻率擾動分析

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

(1)

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

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

(2)

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

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

(3)

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

(4)

在圖1中，ε和φ為軸頸在擾動狀態(tài)下某一瞬時所在位置的偏心率和偏位角； 對應(yīng)的ε0和φ0為軸頸穩(wěn)態(tài)平衡位置的偏心率和偏位角；Ob為軸承中心；Oj和Oj0為瞬時位置和穩(wěn)態(tài)平衡位置的軸頸中心。

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

(5)

(6)

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

(7)

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

(8)

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

(9)

式中，Re0為對應(yīng)于軸頸平衡位置的穩(wěn)態(tài)雷諾數(shù)，擾動雷諾數(shù)Red是擾動密度、擾動黏度和擾動膜厚的聯(lián)合響應(yīng)。

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

(10)

(11)

1.2 穩(wěn)態(tài)和擾動雷諾方程

(12)

(13)

1.3 動力學(xué)系數(shù)的求取方法

(14)

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

式(15)可將動力學(xué)系數(shù)從耦合坐標(biāo)系變換到x-y直角坐標(biāo)系。

(15)

式中， [A]為旋轉(zhuǎn)矩陣，其表達(dá)式為

(16)

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

2 完整變量擾動偏導(dǎo)數(shù)法的驗證

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

(17)

(18)

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

(19)

其中，

(20)

動力學(xué)系數(shù)從ε-φ坐標(biāo)系轉(zhuǎn)換到圖1的x-y坐標(biāo)系公式如式(21)所示，其中[A]與式(16)相同。

(21)

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

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

2.2 對比驗證結(jié)果

表2 S-CO2 圓柱瓦軸承的結(jié)構(gòu)和運(yùn)行參數(shù)Tab.2 Parameters of S-CO2 cylindrical bearing

3 案例計算結(jié)果與討論

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

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

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

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

隨著的增加，阻尼系數(shù)先緩慢減小，大約從Ω>0.5開始明顯減小，只考慮壓力擾動或準(zhǔn)靜態(tài)假設(shè)都導(dǎo)致阻尼系數(shù)結(jié)果明顯偏大。

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

對于可傾瓦軸承，可壓縮時滯效應(yīng)增大了剛度系數(shù)而減小了阻尼系數(shù)。忽略熱物性和湍流的動態(tài)變化效應(yīng)將導(dǎo)致剛度和阻尼系數(shù)明顯偏大。

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

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

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

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

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

4 結(jié) 論

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

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

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

(3) 高速可傾瓦油潤滑軸承在潤滑膜湍流效應(yīng)明顯時，若忽略湍流系數(shù)的動態(tài)變化效應(yīng)將導(dǎo)致小偏心率下的剛度系數(shù)結(jié)果偏大。