雷 群 宋宜強(qiáng) 杜建軍

(1.廣州市昊志機(jī)電股份有限公司 廣東廣州 511356;2.哈爾濱工業(yè)大學(xué)(深圳) 廣東深圳 518055)

近年來，隨著精密加工技術(shù)和超高速旋轉(zhuǎn)機(jī)械的應(yīng)用，人們對其中關(guān)鍵的部件——軸承的性能要求越來越高。從20世紀(jì)起，氣體潤滑技術(shù)發(fā)展迅速，氣體軸承開始被廣泛應(yīng)用。氣體軸承是用氣體代替潤滑油進(jìn)行潤滑、使用氣膜提供支撐力的機(jī)械構(gòu)件，特別適合于輕型、高速的渦輪機(jī)械，例如渦輪壓縮機(jī)和空氣循環(huán)機(jī)等[1-3]。氣體箔片軸承是氣體動壓軸承的一種，由于其出色的性能，近年來越來越受到人們的重視。相比于傳統(tǒng)的剛性氣體動壓軸承，氣體箔片軸承具有壽命長、承載能力大、可靠性高和穩(wěn)定性好的優(yōu)點(diǎn)，并且可以在高溫高速情況下工作。

1968年，LUND[4]在研究氣體動壓軸承時提出了攝動法，即用4個剛度系數(shù)和4個阻尼系數(shù)來計(jì)算氣膜力。這種方法將油膜或氣膜線性化表示，優(yōu)點(diǎn)是建模簡單、計(jì)算速度快，因此獲得了廣泛的應(yīng)用[5-9]。但該方法僅限于小擾動及低速轉(zhuǎn)動情況。軌跡法能夠彌補(bǔ)攝動法這種不足，該方法通過聯(lián)立求解雷諾方程、轉(zhuǎn)子動力學(xué)方程和箔片變形方程，能夠獲得系統(tǒng)在各時刻的運(yùn)動情況，因此有助于更好地研究系統(tǒng)的動力學(xué)特性，從而得到了廣泛應(yīng)用[10-13]。

在使用軌跡法時，一種常用的做法是將雷諾方程、轉(zhuǎn)子動力學(xué)方程和箔片變形方程在每一個時間步內(nèi)分開求解，由于各狀態(tài)量不能同步更新，需要保持積分步長很小才能保證算法收斂，計(jì)算量很大。為了克服弱耦合解法的弱點(diǎn)，本文作者建立了箔片軸承轉(zhuǎn)子系統(tǒng)的強(qiáng)耦合數(shù)學(xué)模型，實(shí)現(xiàn)雷諾方程、轉(zhuǎn)子動力學(xué)方程和箔片變形方程的同步求解，從而使計(jì)算效率得到較大的提高。

1 非線性動力學(xué)耦合模型

圖1所示是第一代氣體箔片軸承的結(jié)構(gòu)示意圖。當(dāng)轉(zhuǎn)子發(fā)生偏移并高速旋轉(zhuǎn)時，轉(zhuǎn)軸與平箔片之間產(chǎn)生楔形間隙，從而為軸承提供承載力。氣膜壓力作用在平箔片上會使平箔片和波箔片產(chǎn)生彈性變形，文中將箔片結(jié)構(gòu)視為普通彈性元件進(jìn)行分析。在受到氣膜力后，箔片結(jié)構(gòu)的變形能夠動態(tài)地調(diào)整氣膜間隙。在求解系統(tǒng)的動態(tài)響應(yīng)時，主要需要考慮3個方面：瞬態(tài)氣膜壓力的產(chǎn)生、轉(zhuǎn)子受氣膜力引起的渦動以及箔片結(jié)構(gòu)受力變形。在建立相應(yīng)的方程時，需要將一些相關(guān)的狀態(tài)量表示出來。

圖1 氣體箔片軸承結(jié)構(gòu)示意

1.1 氣膜壓力求解

在研究氣體軸承轉(zhuǎn)子系統(tǒng)的動力學(xué)特性時，通過求解可壓縮時變雷諾方程得到氣膜壓力。在等溫狀態(tài)下，氣體的時變雷諾方程為

(1)

式中：h為氣膜厚度；p為氣膜壓力；z為軸向坐標(biāo)；μ為氣體黏度；t為時間；ω為轉(zhuǎn)子轉(zhuǎn)速；R為轉(zhuǎn)子半徑。

(2)

式中：θ為周向角度；pa為環(huán)境大氣壓；c為名義間隙；Λ為軸承數(shù)。

將式(1)變換為式(2)是為了方便后面建立狀態(tài)方程。設(shè)箔片變形為q(θ)， 取向外為正，考慮箔片的變形，則氣膜厚度為

h=c+xsinθ-ycosθ+q

(3)

式中：x為轉(zhuǎn)子在x方向的位移；y為轉(zhuǎn)子在y方向的位移；坐標(biāo)軸方向在圖1中顯示。

(4)

使用有限差分法求解式(2)的雷諾方程，網(wǎng)格劃分如圖2所示。

圖2 有限差分法的網(wǎng)格劃分

計(jì)算域采用矩形網(wǎng)格離散為Nx×Nz個點(diǎn)，軸承的兩端作為邊界點(diǎn)，此處的氣壓等于外界大氣壓，不作為狀態(tài)變量。對于氣膜承載力，在計(jì)算中使用量綱一化形式，在考慮方向后積分得到計(jì)算公式：

(5)

(6)

式中:Fx和Fy分別為在x方向和y方向的氣膜承載力;L為箔片軸承寬度。

1.2 轉(zhuǎn)子動力學(xué)方程

轉(zhuǎn)子動力學(xué)方程由牛頓第二定律獲得。系統(tǒng)運(yùn)行時，轉(zhuǎn)子除了受到氣膜力外，還需要考慮由于轉(zhuǎn)子質(zhì)心與形心不重合導(dǎo)致的自轉(zhuǎn)慣性力。因此轉(zhuǎn)子動力學(xué)方程表示為

(7)

(8)

式中:m為轉(zhuǎn)子質(zhì)量;Fub為不平衡力;S為重力;下標(biāo)x和y分別表示在x和y方向。

對式(7)和(8)進(jìn)行量綱一化，得到：

(9)

(10)

1.3 箔片變形方程

對于第一代波箔軸承，箔片沿軸向的變形可忽略不計(jì)[14]，箔片按照簡單彈簧振動模型進(jìn)行分析。根據(jù)文獻(xiàn)[15]，箔片變形方程表示為

(11)

(12)

箔片的剛度計(jì)算通常使用HESHMAT等[16]在1983年提出的公式：

(13)

式中:E為彈性模量；tb為箔片厚度；ν為泊松比；l0為箔片波紋長度的1/2。

式(13)計(jì)算得到的箔片剛度是線剛度，進(jìn)行量綱一化，式(11)對應(yīng)的量綱一箔片面剛度為

(14)

式中：s為波箔節(jié)距。

1.4 系統(tǒng)的狀態(tài)方程與耦合求解

為了建立整個氣體箔片軸承轉(zhuǎn)子系統(tǒng)的狀態(tài)方程，實(shí)現(xiàn)強(qiáng)耦合求解，需要將時變雷諾方程、轉(zhuǎn)子動力學(xué)方程和箔片變形方程以狀態(tài)方程的形式表示。首先對時變雷諾方程進(jìn)行處理，式(2)變換為

(15)

式中：°為哈達(dá)瑪積，表示向量的各項(xiàng)逐項(xiàng)相乘。

式(9)和(10)的轉(zhuǎn)子動力學(xué)方程以向量形式表示為

(16)

式(11)的箔片變形方程同樣以向量形式表示為

(17)

將式(15)、(16)和(17)聯(lián)立起來，便得到箔片軸承轉(zhuǎn)子系統(tǒng)的總體方程，以狀態(tài)方程的形式表示為

(18)

其中，狀態(tài)變量y為

(19)

求解式(18)的系統(tǒng)狀態(tài)方程，便能同時得到瞬態(tài)氣膜壓力、轉(zhuǎn)子位移和箔片變形量，實(shí)現(xiàn)強(qiáng)耦合求解，計(jì)算效率相比弱耦合解法有了較大提升。

2 仿真結(jié)果與分析

如無特別指出，文中仿真所選用的箔片軸承和轉(zhuǎn)子的基本參數(shù)如表1所示。

表1 箔片軸承和轉(zhuǎn)子的基本參數(shù)

2.1 模型驗(yàn)證

設(shè)轉(zhuǎn)子初始位置的中心線與軸承的軸心重合，初始平動速度為0，轉(zhuǎn)速為10 000 r/min，計(jì)算從量綱一時間為0～10π的轉(zhuǎn)子運(yùn)動軌跡與箔片變形，并與文獻(xiàn)[17]中的數(shù)據(jù)比較來驗(yàn)證模型的正確性，如圖3所示。

圖3 轉(zhuǎn)速為10 000 r/min時仿真結(jié)果與文獻(xiàn)結(jié)果比較

可以看到，文中仿真得到的轉(zhuǎn)子軌跡、箔片變形與文獻(xiàn)[17]中的結(jié)果基本一致，說明文中建立的箔片軸承轉(zhuǎn)子系統(tǒng)動力學(xué)模型是可行的。

2.2 轉(zhuǎn)子轉(zhuǎn)速對系統(tǒng)動力學(xué)特性的影響

轉(zhuǎn)子轉(zhuǎn)速對系統(tǒng)的穩(wěn)定性和其他動力學(xué)特性產(chǎn)生重要影響。以轉(zhuǎn)速為研究對象，使用不同轉(zhuǎn)速，計(jì)算系統(tǒng)的運(yùn)動狀態(tài)。分別取轉(zhuǎn)速為9 000、10 500和12 000 r/min，轉(zhuǎn)子不平衡偏心量為1 μm，繪制系統(tǒng)進(jìn)入穩(wěn)定狀態(tài)后，轉(zhuǎn)子的運(yùn)動軌跡和x方向的頻譜圖，如圖4—6所示。

圖4(a)—6(a)所示是轉(zhuǎn)子的運(yùn)動軌跡圖，其橫坐標(biāo)為轉(zhuǎn)子軸心的水平偏心率，縱坐標(biāo)為軸心的豎直偏心率?？梢钥闯?，當(dāng)轉(zhuǎn)子轉(zhuǎn)速為9 000 r/min時，轉(zhuǎn)子的運(yùn)動軌跡是一個規(guī)則的橢圓形；而在轉(zhuǎn)速提高為10 500 r/min時，轉(zhuǎn)子軌跡略微紊亂，但仍有跡可循；但當(dāng)轉(zhuǎn)速提升到12 000 r/min時，轉(zhuǎn)子軌跡圖變得復(fù)雜混亂，難以對其走向進(jìn)行估計(jì)。

圖4(b)—6(b)所示是對轉(zhuǎn)子軸心在x方向的運(yùn)動做快速傅里葉變換得到的頻譜圖。可以看出，當(dāng)轉(zhuǎn)速為9 000 r/min時，轉(zhuǎn)子渦動與轉(zhuǎn)速同頻渦動；當(dāng)提高轉(zhuǎn)速到10 500 r/min時，除了1倍頻外，還出現(xiàn)了0.3倍頻和0.7倍頻；當(dāng)轉(zhuǎn)速為12 000 r/min時，同樣出現(xiàn)了多種頻率，頻率分量包含0.35倍頻、0.65倍頻和1倍頻，且0.65倍頻分量的振幅接近于1倍頻分量的振幅，此時系統(tǒng)接近失穩(wěn)。據(jù)頻譜圖可以知道，在轉(zhuǎn)速較低時，轉(zhuǎn)子做穩(wěn)定的周期運(yùn)動，但是當(dāng)轉(zhuǎn)速提高后，頻譜圖會出現(xiàn)其他頻率分量，轉(zhuǎn)子可能做擬周期運(yùn)動。

圖4 轉(zhuǎn)速為9 000 r/min時轉(zhuǎn)子的運(yùn)動狀態(tài)

圖5 轉(zhuǎn)速為10 500 r/min時轉(zhuǎn)子的運(yùn)動狀態(tài)

圖6 轉(zhuǎn)速為12 000 r/min時轉(zhuǎn)子的運(yùn)動狀態(tài)

圖7所示是以轉(zhuǎn)速為分岔參數(shù)繪制的轉(zhuǎn)子在x方向和y方向的分岔圖。分岔圖是描繪狀態(tài)變量隨分岔參數(shù)變化的關(guān)系?？梢钥闯?，當(dāng)轉(zhuǎn)速達(dá)到9 500 r/min時，轉(zhuǎn)子的運(yùn)動出現(xiàn)分叉，說明此時轉(zhuǎn)子進(jìn)入擬周期運(yùn)動狀態(tài)。

圖7 以轉(zhuǎn)速為分岔參數(shù)的分岔圖

2.3 轉(zhuǎn)子質(zhì)量對系統(tǒng)動力學(xué)特性的影響

改變轉(zhuǎn)子質(zhì)量進(jìn)行數(shù)值仿真計(jì)算。取轉(zhuǎn)子轉(zhuǎn)速為10 000 r/min，轉(zhuǎn)子質(zhì)量分別為2、4和6 kg，繪制其運(yùn)動軌跡圖和頻譜圖進(jìn)行分析，如圖8—10所示。

圖8 轉(zhuǎn)子質(zhì)量為2 kg時轉(zhuǎn)子的運(yùn)動狀態(tài)

圖9 轉(zhuǎn)子質(zhì)量為4 kg時轉(zhuǎn)子的運(yùn)動狀態(tài)

圖10 轉(zhuǎn)子質(zhì)量為6 kg時轉(zhuǎn)子的運(yùn)動狀態(tài)

圖8(a)—10(a)顯示了轉(zhuǎn)子的軌跡?？芍?，當(dāng)轉(zhuǎn)子質(zhì)量為2 kg時，轉(zhuǎn)子軌跡為規(guī)則的橢圓形；而在轉(zhuǎn)子質(zhì)量為4 kg時，軌跡變得比較復(fù)雜，且轉(zhuǎn)子的振幅增大；當(dāng)轉(zhuǎn)子質(zhì)量增大到6 kg時，轉(zhuǎn)子軌跡進(jìn)一步紊亂。這個規(guī)律與增加轉(zhuǎn)子轉(zhuǎn)速情況相似。

圖8(b)—10(b)所示是對轉(zhuǎn)子軸心在x方向的運(yùn)動做快速傅里葉變換得到的頻譜圖。可以看出，當(dāng)轉(zhuǎn)子質(zhì)量為2 kg時，轉(zhuǎn)子渦動與轉(zhuǎn)速同頻渦動；在轉(zhuǎn)子質(zhì)量為4 kg時，除了基頻，頻譜圖上還出現(xiàn)0.35倍頻和0.65倍頻；當(dāng)轉(zhuǎn)子質(zhì)量為6 kg時，頻譜圖上出現(xiàn)0.45倍頻、0.55倍頻和1倍頻。據(jù)此規(guī)律分析，當(dāng)轉(zhuǎn)子質(zhì)量進(jìn)一步增加時，將會出現(xiàn)0.5倍頻的渦動，即為半速渦動。

圖11所示是以轉(zhuǎn)子質(zhì)量為分岔參數(shù)繪制的分岔圖。可知，當(dāng)轉(zhuǎn)子質(zhì)量小于等于3 kg時，轉(zhuǎn)子做周期運(yùn)動；但是當(dāng)轉(zhuǎn)子質(zhì)量大于3 kg時，開始出現(xiàn)其他頻率的振動；隨著轉(zhuǎn)子質(zhì)量的進(jìn)一步增加，系統(tǒng)最終將失穩(wěn)。

圖11 以轉(zhuǎn)子質(zhì)量為分岔參數(shù)的分岔圖

2.4 名義間隙對系統(tǒng)動力學(xué)特性的影響

以名義間隙為研究對象，分別取名義間隙為25、35和50 μm進(jìn)行仿真，計(jì)算結(jié)果如圖12—14所示。

圖12(a)—14(a)顯示了轉(zhuǎn)子的軌跡?？芍x間隙為25 μm時，系統(tǒng)不穩(wěn)定，轉(zhuǎn)子運(yùn)動情況較混亂；而在名義間隙為35 μm時，轉(zhuǎn)子的運(yùn)動狀態(tài)得到改善，運(yùn)動軌跡偏移范圍較??；當(dāng)名義間隙為50 μm時，運(yùn)動振幅進(jìn)一步減小。

圖12(b)—14(b)所示是對轉(zhuǎn)子軸心在x方向的運(yùn)動做快速傅里葉變換得到的頻譜圖?？梢钥闯觯?dāng)名義間隙為25 μm時，頻譜圖上除同頻分量外，還出現(xiàn)0.2倍頻和0.8倍頻；在名義間隙為35 μm時，頻譜圖上主要頻率為基頻分量，其他另有0.7倍頻和0.75倍頻，但僅占很小部分；當(dāng)名義間隙為50 μm時，主要頻率成分仍為基頻，還出現(xiàn)0.2倍頻和0.8倍頻?？梢灾溃?dāng)名義間隙增大到某一值時，再增加名義間隙對系統(tǒng)穩(wěn)定性影響不大。

圖12 名義間隙為25 μm時轉(zhuǎn)子的運(yùn)動狀態(tài)

圖13 名義間隙為35 μm時轉(zhuǎn)子的運(yùn)動狀態(tài)

圖14 名義間隙為50 μm時轉(zhuǎn)子的運(yùn)動狀態(tài)

3 結(jié)論

(1)針對氣體箔片軸承轉(zhuǎn)子系統(tǒng)的動態(tài)響應(yīng)進(jìn)行研究，建立同步求解雷諾方程、轉(zhuǎn)子動力學(xué)方程和箔片變形方程的強(qiáng)耦合數(shù)學(xué)模型，相比傳統(tǒng)的弱耦合求解方法，能夠進(jìn)一步提升計(jì)算速度和提高計(jì)算精度，為箔片軸承的設(shè)計(jì)和性能預(yù)測提供了研究工具。

(2)探究不同設(shè)計(jì)參數(shù)對軸承轉(zhuǎn)子系統(tǒng)性能的影響。通過改變轉(zhuǎn)速、轉(zhuǎn)子質(zhì)量和名義間隙等參數(shù)進(jìn)行數(shù)值仿真，并繪制轉(zhuǎn)子運(yùn)動軌跡圖及頻譜圖，分析系統(tǒng)的穩(wěn)定性。結(jié)果顯示，當(dāng)轉(zhuǎn)子的質(zhì)量增加或轉(zhuǎn)速提高時，系統(tǒng)的穩(wěn)定性呈下降趨勢；反之，系統(tǒng)穩(wěn)定性增強(qiáng)。另外，當(dāng)名義間隙過小時，系統(tǒng)容易出現(xiàn)失穩(wěn)狀況，當(dāng)名義間隙增大到某一值時，再增加名義間隙對系統(tǒng)穩(wěn)定性影響不大。