韓中合，侯棟楠，趙建宏，朱霄珣

(華北電力大學(xué) 能源動(dòng)力與機(jī)械工程學(xué)院，河北 保定 071003)

隨著氣體潤(rùn)滑技術(shù)在20世紀(jì)中期迅速發(fā)展，氣體軸承作為該技術(shù)的核心產(chǎn)品被開發(fā)出來。氣體軸承與滾動(dòng)軸承和液體滑動(dòng)軸系相比具有高精度、高轉(zhuǎn)速、振動(dòng)和噪聲小、壽命長(zhǎng)和摩擦功耗小等優(yōu)點(diǎn)，在高速旋轉(zhuǎn)機(jī)械上得到了廣泛應(yīng)用[1]。同時(shí)氣體軸承存在著承載能力低，穩(wěn)定性差等缺點(diǎn)。氣體軸承作為氣體軸承-轉(zhuǎn)子系統(tǒng)的支撐部分，直接影響轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)特性[2]。因此需要將氣體軸承和轉(zhuǎn)子耦合起來進(jìn)行分析，考慮氣體軸承對(duì)整個(gè)系統(tǒng)的影響。

1886年,Reynolds推導(dǎo)出解決流體潤(rùn)滑軸承理論計(jì)算的Reynolds方程，為解決流體潤(rùn)滑問題奠定了理論基礎(chǔ)。隨著計(jì)算機(jī)技術(shù)的進(jìn)步，通常采用數(shù)值計(jì)算求解Reynolds方程，得到氣體軸承的性能參數(shù)[3]。隨著計(jì)算流體力學(xué)(CFD)技術(shù)的發(fā)展，采用CFD軟件對(duì)氣體軸承內(nèi)流場(chǎng)進(jìn)行模擬，也可以得到氣體軸承的性能參數(shù)。對(duì)比采用有限差分法求解靜壓氣體Reynolds方程和采用CFD計(jì)算的結(jié)果，可知兩者具有較好的一致性[4]。同時(shí)利用CFD軟件對(duì)氣體軸承的性能進(jìn)行計(jì)算，并根據(jù)計(jì)算結(jié)果研究氣體軸承-轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)特性，也能得到與試驗(yàn)相對(duì)一致的結(jié)論[5]。在試驗(yàn)中，發(fā)現(xiàn)通過升高軸承供氣壓力，能夠提高氣體軸承剛度，消除和抑制轉(zhuǎn)子低頻振動(dòng)并增大氣體軸承-轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速[6]。綜上所述，采用CFD和求解Reynolds方程的方法計(jì)算靜壓氣體軸承的性能，都能得到相對(duì)準(zhǔn)確的結(jié)果。且氣體軸承的供氣壓力、軸頸轉(zhuǎn)速等因素都對(duì)氣體軸承-轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)特性產(chǎn)生影響。

采用CFD方法計(jì)算氣體軸承性能雖然較為精確，但使用CFD計(jì)算方法需要對(duì)不同工況下流場(chǎng)的邊界條件進(jìn)行逐一設(shè)置并計(jì)算。在研究氣體軸承的性能對(duì)氣體軸承-轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)特性影響時(shí)，需要得到供氣壓力、轉(zhuǎn)速變化范圍內(nèi)任意工況的氣體軸承的性能。而采用CFD計(jì)算方法對(duì)所需工況進(jìn)行計(jì)算耗時(shí)漫長(zhǎng)、效率較低，不利于為系統(tǒng)動(dòng)力學(xué)研究提供支撐參數(shù)。

針對(duì)該問題本文提出了采用支持向量機(jī)回歸算法(SVR)建立氣體軸承性能函數(shù)，根據(jù)此函數(shù)對(duì)氣體軸承-轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速計(jì)算方法。該方法只需對(duì)多個(gè)由CFD計(jì)算得到的氣體軸承的性能參數(shù)進(jìn)行訓(xùn)練，就能得到準(zhǔn)確的氣體軸承性能函數(shù)，此函數(shù)能快速、準(zhǔn)確地預(yù)測(cè)一定壓力、轉(zhuǎn)速范圍內(nèi)任意工況下的軸承性能，利用該模型對(duì)供氣壓力變化時(shí)氣體軸承-轉(zhuǎn)子系統(tǒng)進(jìn)行動(dòng)力學(xué)分析，能夠得到和試驗(yàn)接近的結(jié)果，為選擇合理的供氣方案提供參考。

1 基于fluent的靜壓氣體軸承性能計(jì)算

1.1 靜壓氣體軸承的結(jié)構(gòu)尺寸

本文研究靜壓氣體軸承結(jié)構(gòu)如圖1所示，結(jié)構(gòu)參數(shù)如表1所示。

圖1 靜壓氣體軸承結(jié)構(gòu)Fig.1 Static pressure gas bearing structure

1.2 模型的建立與假設(shè)

根據(jù)圖1與表1的結(jié)構(gòu)尺寸，利用前處理軟件GAMBIT建立靜壓氣體軸承的流道模型(圖2)。

氣體軸承的供氣孔區(qū)是計(jì)算的關(guān)鍵部位，需要對(duì)其進(jìn)行網(wǎng)格加密，因此將流場(chǎng)劃分為供氣孔區(qū)和無供氣孔區(qū)。在供氣孔區(qū)采用非結(jié)構(gòu)化網(wǎng)格(圖3)，在無供氣孔區(qū)采用結(jié)構(gòu)化網(wǎng)格。為提高流場(chǎng)計(jì)算精度，在氣膜的高度方向?qū)⒕W(wǎng)格分為10層。

對(duì)模型設(shè)置3類邊界條件：① 壓力入口，供氣壓力pin=0.25 MPa、0.35 MPa、0.45 MPa、0.6 MPa;② 壓力出口，兩側(cè)出氣口壓力等于大氣壓力，相對(duì)壓力pout=0 MPa;③ 旋轉(zhuǎn)面，選取軸頸狀態(tài)為靜止、10 000 r/min、20 000 r/min、…、60 000 r/min。

表1 靜壓氣體軸承結(jié)構(gòu)參數(shù)Tab.1 Static pressure gas bearing structural parameters

圖2 氣體軸承網(wǎng)格劃分Fig.2 Mesh generation of gas bearings圖3 供氣孔區(qū)網(wǎng)格Fig.3 Air supply hole area grid

實(shí)際的靜壓氣體軸承流場(chǎng)較為復(fù)雜，所以在計(jì)算時(shí)做如下假設(shè)：① 軸承間隙內(nèi)的氣體是理想氣體，氣體與壁面無熱量交換，即為絕熱過程;② 軸承與軸頸工作面為光滑曲面，氣膜厚度為常數(shù)，不考慮粗糙度對(duì)軸承特性的影響[7];③ 軸頸面沿y軸負(fù)方向偏移，且不考慮偏心角的影響。

本文采用Fluent軟件對(duì)流場(chǎng)數(shù)值仿真，選用基于壓力的隱式求解器。由于軸承氣膜中存在著不同尺度的漩渦[8]，因此采用湍流模型中的模型，該模型能夠較準(zhǔn)確模擬軸承間隙內(nèi)的流場(chǎng)[9]。選擇工質(zhì)為空氣，方程采用SMPLE算法求解。

1.3 計(jì)算結(jié)果與試驗(yàn)對(duì)比

采用與上一節(jié)相同的網(wǎng)格劃分、邊界條件和計(jì)算假設(shè)條件，建立文獻(xiàn)[10]中試驗(yàn)所采用的氣體軸承流道模型，得到CFD計(jì)算與試驗(yàn)結(jié)果如圖4所示。

結(jié)果表明，采用Fluent對(duì)氣體軸承的承載力計(jì)算的結(jié)果與實(shí)際較為吻合。隨著偏心率增大到0.6，仿真結(jié)果與試驗(yàn)出現(xiàn)偏差。這是由于隨著偏心率的增大，氣體軸承一側(cè)間隙較小，實(shí)際流場(chǎng)和層流模型存在差異導(dǎo)致的。

根據(jù)文獻(xiàn)[11]，本文研究的氣體軸承-轉(zhuǎn)子系統(tǒng)中，轉(zhuǎn)子正常運(yùn)行時(shí)的最大偏心率小于0.4，因此采用Fluent對(duì)氣體軸承的性能進(jìn)行計(jì)算可以得到較準(zhǔn)確結(jié)果。

圖4 仿真結(jié)果和實(shí)際對(duì)比Fig.4 Simulation results and the actual comparison

1.4 氣體軸承承載力計(jì)算

考慮到氣體的供氣壓力會(huì)對(duì)氣體軸承的剛度產(chǎn)生影響，設(shè)置0.25 MPa、0.35 MPa、0.45 MPa、0.55 MPa的供氣壓力。設(shè)軸頸沿y軸負(fù)方向偏移，偏心率為0～0.4。轉(zhuǎn)速設(shè)置為0～60 000 r/min。分別計(jì)算軸頸x、y方向上的合力。再由式(1)得到軸承的承載力

(1)

式中:f為承載力，fx為軸頸x方向上的合力，fy為軸頸y方向上的合力。

繪制不同供氣壓力下軸承承載力隨轉(zhuǎn)速和偏心率變化的曲線圖。(如圖5～圖8所示)。

圖5 0.25 MPa時(shí)軸承承載力Fig.5 Bearing capacity of bearing at 0.25 MPa

圖6 0.35 MPa時(shí)軸承承載力Fig.6 Bearing capacity of bearing at 0.35 MPa

圖7 0.45 MPa時(shí)軸承承載力Fig.7 Bearing capacity of bearing at 0.45 MPa

圖8 0.6 MPa時(shí)軸承承載力Fig.8 Bearing capacity of bearing at 0.6 MPa

分析承載力變化規(guī)律得到，相同的供氣壓力下，承載力隨著偏心率的增大近似線性增大。隨著轉(zhuǎn)速的升高，承載力也增大。在供氣壓力升高時(shí)，承載力和轉(zhuǎn)速與偏心率的變化規(guī)律不變，但承載力整體增大，這是由于供氣壓力升高，在相同的偏心率下氣膜兩側(cè)的壓差增大，產(chǎn)生的承載力也增大。

2 基于SVM的靜壓氣體軸承性能計(jì)算

支持向量機(jī)(Support Vector Machine,SVM)是一種分類算法，通過尋求結(jié)構(gòu)化風(fēng)險(xiǎn)最小來提高學(xué)習(xí)機(jī)泛化能力，實(shí)現(xiàn)經(jīng)驗(yàn)風(fēng)險(xiǎn)和置信范圍的最小化，從而達(dá)到在樣本量較少的情況下，亦能獲得良好統(tǒng)計(jì)規(guī)律的目的[12]。通俗來講，它是一種二類分類模型，其基本模型定義為特征空間上的間隔最大的線性分類器，即支持向量機(jī)的學(xué)習(xí)策略便是間隔最大化，最終可轉(zhuǎn)化為一個(gè)凸二次規(guī)劃問題的求解。

2.1 氣體軸承承載力計(jì)算

對(duì)于本文模型，首先通過CFD計(jì)算，得到不同壓力p、轉(zhuǎn)速n、偏心率e下的軸承承載力，并將其作為SVM的訓(xùn)練樣本集：T={(x1,y1),(x2,y2),…,(xn,yn)}，其中，x=(p，n，e)，y為承載力F限定回歸問題所需要尋找的函數(shù)f(x)為線性函數(shù)

y=f(x)=ωx+b

(2)

然后構(gòu)造并求解最優(yōu)化問題

(3)

(4)

(5)

式中：xi為樣本輸入；yi為樣本輸出；ω為權(quán)重向量；b為偏置項(xiàng)；C為懲罰系數(shù)；αi、αi*為拉格朗日乘子；ε為不敏感損失值。最后構(gòu)造出最優(yōu)線性回歸函數(shù)[13]。

在實(shí)際情況下對(duì)線性問題進(jìn)行拓廣后可得到非線性情況下的回歸函數(shù)，其關(guān)鍵是引進(jìn)核函數(shù)K(xi，xj)，把數(shù)據(jù)從低維空間映射到高維空間[14]，優(yōu)化問題轉(zhuǎn)化為

(6)

最終可得到最優(yōu)的非線性回歸函數(shù)。

在選取合適的懲罰系數(shù)以及核函數(shù)后利用訓(xùn)練樣本建立支持向量機(jī)(SVM)預(yù)測(cè)模型[15]。模型建立好以后，利用預(yù)測(cè)模型分別預(yù)測(cè)已經(jīng)輸入到SVM的訓(xùn)練集數(shù)據(jù)和未輸入到SVM的數(shù)據(jù)。再利用CFD方法計(jì)算進(jìn)行預(yù)測(cè)工況條件下的結(jié)果，并進(jìn)行對(duì)比。結(jié)果見圖9、圖10。

圖9 SVM預(yù)測(cè)訓(xùn)練集數(shù)據(jù)Fig.9 SVM Prediction of Training Set Data

圖10 SVM預(yù)測(cè)未輸入數(shù)據(jù)Fig.10 SVM prediction of not entered data

從圖9、圖10可以看出預(yù)測(cè)模型計(jì)算出的數(shù)據(jù)和輸入模型的已知數(shù)據(jù)以及未輸入模型的數(shù)據(jù)對(duì)比準(zhǔn)確率都很高。根據(jù)計(jì)算，采用SVM預(yù)測(cè)訓(xùn)練集的最大誤差為0.4%，預(yù)測(cè)未輸入到SVM的數(shù)據(jù)最大誤差為0.8%。該結(jié)果滿足對(duì)軸承性能計(jì)算的要求。

2.2 氣體軸承剛度計(jì)算

在圖1所示的靜壓氣體軸承中，轉(zhuǎn)子到受氣體軸承承載力實(shí)現(xiàn)無接觸懸浮，所以氣體軸承的剛度直接影響轉(zhuǎn)子的支撐剛度，為了方便計(jì)算轉(zhuǎn)子臨界轉(zhuǎn)速，氣體軸承-轉(zhuǎn)子系統(tǒng)的支撐剛度由式(7)得到[16]。

(7)

式中:k為系統(tǒng)支撐剛度，F(xiàn)為軸承承載力，h為軸頸偏移量。

通過CFD計(jì)算無法得到承載力F關(guān)于偏移量h的函數(shù)，而通過SVM建立的模型也難以直接進(jìn)行求導(dǎo)。故在實(shí)際過程中采用差分代替微分計(jì)算支撐剛度k(如式(3))

(8)

當(dāng)Δh值越小，式(8)的計(jì)算結(jié)果越接近與式(7)。但當(dāng)Δh值過小時(shí)，會(huì)因承載力F的擬合誤差導(dǎo)致計(jì)算精度下降，因此本文選取Δh值為1 μm。

通過由SVM建立的模型，利用式(8)計(jì)算出轉(zhuǎn)子的剛度隨轉(zhuǎn)速偏心率、轉(zhuǎn)速、供氣壓力的變化規(guī)律如圖11～圖13所示。

根據(jù)圖11可得，在偏心率小于0.4的情況下，轉(zhuǎn)子的支撐剛度受偏心率的影響較小。

由圖12可得，當(dāng)轉(zhuǎn)子的轉(zhuǎn)速小于30 000 r/min時(shí)，轉(zhuǎn)速對(duì)轉(zhuǎn)子的支撐剛度影響也很小，當(dāng)轉(zhuǎn)速大于30 000 r/min時(shí)，隨著轉(zhuǎn)速升高轉(zhuǎn)子支撐剛度增大，而且供氣壓力越大，轉(zhuǎn)子支撐剛度受轉(zhuǎn)速影響越大。

由圖13可得，相對(duì)于偏心率、轉(zhuǎn)速，供氣壓力對(duì)轉(zhuǎn)子支撐剛度的影響最大。供氣壓力在0.25～0.55 MPa內(nèi)轉(zhuǎn)子支撐剛度隨供氣壓力近似線性增大，當(dāng)供氣壓力大于0.55 MPa后，轉(zhuǎn)子支撐剛度對(duì)供氣壓力不敏感。

圖11 60 000 r/min時(shí)轉(zhuǎn)子支撐剛度隨偏心率變化
Fig.11 Rotor support stiffness changes with eccentricity at 60 000 r/min

圖12 偏心率0.3時(shí)轉(zhuǎn)子支撐剛度隨轉(zhuǎn)速變化Fig.12 Rotor support stiffness changes with speed at eccentricity 0.3

圖13 偏心率0.3時(shí)轉(zhuǎn)子支撐剛度隨供氣壓力變化Fig.13 Rotor support stiffness changes with air supply pressure at eccentricity 0.3

3 供氣壓力對(duì)轉(zhuǎn)子臨界轉(zhuǎn)速的影響

3.1 模態(tài)分析的理論基礎(chǔ)

模態(tài)分析是求得系統(tǒng)的各階模態(tài)參數(shù)，對(duì)一個(gè)N自由度的線性振動(dòng)系統(tǒng)，根據(jù)彈性力學(xué)有限元理論，得到其運(yùn)動(dòng)微分方程

(9)

在對(duì)系統(tǒng)進(jìn)行模態(tài)分析時(shí)，若無外力作用設(shè)f(t)=0,得到系統(tǒng)的自由振動(dòng)方程。同時(shí)在求解系統(tǒng)的固有頻率時(shí)，阻尼對(duì)氣體軸承影響不大，忽略阻尼對(duì)系統(tǒng)的影響得到

(10)

設(shè)式(10)的解為簡(jiǎn)諧振動(dòng)

x=X0sin(ωt+θ)

(11)

將式(11)代入得到

([K]-ω2[M])X0=0

(12)

式中:ω為轉(zhuǎn)子的固有頻率,X0為轉(zhuǎn)子的振幅矩陣。求解此方程即可得到轉(zhuǎn)子的固有頻率。

3.2 利用Ansys建立氣體軸承-轉(zhuǎn)子系統(tǒng)的模型

氣體軸承-轉(zhuǎn)子系統(tǒng)如圖14所示，主要包括1徑向氣體軸承；2止推軸承。長(zhǎng)度L為582 mm，直徑D為50 mm。

圖14 氣體軸承-轉(zhuǎn)子結(jié)構(gòu)Fig.14 Gas bearing-rotor structure

根據(jù)圖12所示的結(jié)構(gòu)尺寸，利用ANSYS前處理器建立模型。定義單元的類型和材料屬性。選取solid185單元對(duì)主軸進(jìn)行劃分，選用combin14彈簧單元對(duì)徑向氣體軸承進(jìn)行模擬。其中combin14單元可以設(shè)置彈性剛度和阻尼系數(shù)，由于氣體軸承阻尼系數(shù)較小，所以忽略阻尼系數(shù)的影響，根據(jù)支持向量機(jī)計(jì)算的結(jié)果得到不同工況下combin14元的彈性剛度。本文不考慮止推軸承的性能對(duì)轉(zhuǎn)子固有頻率的影響，在止推盤處添加位移約束模擬止推軸承。

轉(zhuǎn)子的主體采用的材料為1Cr13結(jié)構(gòu)，可以用于各種旋轉(zhuǎn)機(jī)械。其為密度7 850 kg/m3，彈性模量210 GPa。泊松比為0.3。為提高計(jì)算效率，在劃分網(wǎng)格時(shí)忽略倒角和圓角，并采用“自由”網(wǎng)格劃分方法對(duì)網(wǎng)格進(jìn)行劃分。最終在Ansys中建立如圖15所示模型。

圖15 氣體軸承-轉(zhuǎn)子模型Fig.15 Gas bearing - rotor model

3.3 計(jì)算結(jié)果與分析

本文采用Block Lanczos法提取轉(zhuǎn)子前三階模態(tài)，并對(duì)轉(zhuǎn)子設(shè)置0、10 000、…、60 000 r/min的旋轉(zhuǎn)速度，根據(jù)計(jì)算的氣體軸承特性，設(shè)置不同轉(zhuǎn)速下combin14單元的剛度值。通過固有頻率和等速線的交點(diǎn)計(jì)算出轉(zhuǎn)子的臨界轉(zhuǎn)速。

通過改變氣體軸承的供氣壓力，設(shè)置不同工況下combin14單元的剛度值，分別計(jì)算出供氣壓力在0.25～0.6 MPa間變化時(shí)的轉(zhuǎn)子前三階臨界轉(zhuǎn)速。結(jié)果如圖16所示。

圖16 臨界轉(zhuǎn)速隨供氣壓力變化Fig.16 The critical speed varies with the supply pressure

根據(jù)文獻(xiàn)[5]中用相同轉(zhuǎn)子結(jié)構(gòu)在供氣壓力0.6 MPa下試驗(yàn)得到的轉(zhuǎn)子升速幅頻特性曲線，可知轉(zhuǎn)子一階臨界轉(zhuǎn)速在16 580～17 360 r/min，二階臨界轉(zhuǎn)速在23 400～24 300 r/min，三階臨界轉(zhuǎn)速發(fā)生在43 300 r/min之后。

本文通過計(jì)算得到在0.6 MPa供氣壓力下轉(zhuǎn)子一階臨界轉(zhuǎn)速18 130 r/min，二階臨界轉(zhuǎn)速24 102 r/min，三階臨界轉(zhuǎn)速50 273 r/min。計(jì)算結(jié)果和試驗(yàn)數(shù)據(jù)間的存在差距，這可能是由于試驗(yàn)采用的軸承和計(jì)算采用軸承結(jié)構(gòu)存在的差異造成的。

計(jì)算結(jié)果表明，隨著供氣壓力的提高，轉(zhuǎn)子的一、二階臨界轉(zhuǎn)速發(fā)生升高，三階臨界轉(zhuǎn)速近似不變。同時(shí)當(dāng)供氣壓力增大到0.5 MPa以上時(shí)，供氣壓力對(duì)臨界轉(zhuǎn)速的影響較小。

4 結(jié) 論

(1)本文利用CFD商業(yè)軟件Fluent建立靜壓氣體軸承流場(chǎng)模型，在相同供氣壓力下，軸承承載力隨偏心率增大近似線性增大，也隨轉(zhuǎn)速升高而增大；再供氣壓力增大情況下，承載力隨著轉(zhuǎn)速和偏心變化規(guī)律不變，但承載力整體增大。

(2)采用SVM對(duì)由CFD得到的氣體軸承性能數(shù)據(jù)進(jìn)行訓(xùn)練，得到的氣體軸承性能參數(shù)預(yù)測(cè)模型具有很高的精度，能夠滿足計(jì)算要求。

(3)利用SVM得到的模型進(jìn)行計(jì)算，得到轉(zhuǎn)子支撐剛度隨偏心率、轉(zhuǎn)速、供氣壓力的變化規(guī)律，轉(zhuǎn)子支撐剛度在0.4偏心率以下變化較小，當(dāng)轉(zhuǎn)速在30 000 r/min以下時(shí)變化較小，在30 000 r/min以上時(shí)轉(zhuǎn)子支撐剛度隨轉(zhuǎn)速升高而升高。當(dāng)供氣壓力從0.25 MPa增大到0.5 MPa，轉(zhuǎn)子支撐剛度隨供氣壓力線性增大，當(dāng)供氣壓力增大到0.5 MPa以上時(shí)，轉(zhuǎn)子支撐剛度受供氣壓力影響較小。

(4)利用ANSYS建立氣體軸承-轉(zhuǎn)子系統(tǒng)的有限元模型得到，當(dāng)供氣壓力在0.25～6 MPa區(qū)間內(nèi)升高，氣體軸承-轉(zhuǎn)子系統(tǒng)的一、二階臨界轉(zhuǎn)速增大，三階臨界轉(zhuǎn)速近似不變。研究結(jié)果表明通過改變供氣壓力能夠調(diào)節(jié)系統(tǒng)的臨界轉(zhuǎn)速，使系統(tǒng)的臨界轉(zhuǎn)速避開工作轉(zhuǎn)速。