賈晨輝，高靖，邱明

（河南科技大學(xué) 機(jī)電工程學(xué)院，河南 洛陽 471003）

混合氣體軸承結(jié)合了靜壓與動壓軸承的優(yōu)點，利用靜壓供氣方式使軸承轉(zhuǎn)子在啟動階段形成靜壓承載能力，避免發(fā)生干摩擦；隨著轉(zhuǎn)速增加，利用螺旋槽產(chǎn)生的楔形動壓效應(yīng)形成動壓承載能力，從而避免持續(xù)供給軸承所需的高壓氣體［1-3］。動靜壓氣體軸承在同類軸承中具有較高的承載力。

氣體軸承運轉(zhuǎn)時的動態(tài)特性直接影響轉(zhuǎn)子的非線性動力學(xué)行為及穩(wěn)定性，該特性可通過氣膜的動態(tài)剛度和動態(tài)阻尼系數(shù)反映［4］，從而為軸承轉(zhuǎn)子系統(tǒng)的穩(wěn)定性分析提供理論基礎(chǔ)。

通常采用數(shù)值求解的方法研究氣體軸承的靜動態(tài)特性。文獻(xiàn)［4］采用有限元法對動靜壓軸承進(jìn)行優(yōu)化分析；文獻(xiàn)［5］采用有限差分法對靜壓軸承進(jìn)行性能研究?？紤]到數(shù)值求解較為復(fù)雜，文中采用有限差分法計算氣膜的動態(tài)剛度和阻尼系數(shù)，在不同偏心率下研究供氣壓力和供氣切向角度對氣膜動態(tài)特性系數(shù)的影響規(guī)律。

1 球面氣體軸承結(jié)構(gòu)

球面螺旋槽動靜壓氣體軸承結(jié)構(gòu)如圖1所示。軸承由定子與轉(zhuǎn)子組成，供氣孔在定子上，轉(zhuǎn)子上有螺旋槽，供氣孔位置與螺旋槽的分布不相互影響。圖中，β為螺旋角；ω為軸頸轉(zhuǎn)速；φ為供氣切向角；α0為轉(zhuǎn)子小端角；α1為螺旋槽起始端角；α2為轉(zhuǎn)子大端角；ps為小孔供氣壓力；br為臺寬；bg為槽寬；hg為槽區(qū)間隙；h0為臺域內(nèi)氣膜平均間隙。

圖1 球面氣體軸承結(jié)構(gòu)Fig.1 Structure of spherical gas bearing

2 潤滑分析數(shù)學(xué)模型

以氣體潤滑運動方程、連續(xù)性方程、狀態(tài)方程以及節(jié)流孔流量方程為基礎(chǔ)［4，6］，結(jié)合Newton黏性定律和Reynolds方程假設(shè)，在球面坐標(biāo)系（圖2）下，推導(dǎo)供氣切向角可變的球面動靜壓氣體軸承非線性動態(tài)潤滑分析數(shù)學(xué)模型［7-10］。

圖2 球面坐標(biāo)系Fig.2 Spherical coordinate

圖2中旋轉(zhuǎn)面上任意一點M（?，α，r），其中，?為周向坐標(biāo)；α為沿旋轉(zhuǎn)素線的子午線方向坐標(biāo)；r為M與旋轉(zhuǎn)素線的距離。

通過氣體潤滑運動方程推導(dǎo)出球坐標(biāo)系下沿周向和沿子午線方向的質(zhì)量流量方程分別為

式中：ρ為氣體密度；h為氣膜厚度；μ為氣體黏度；p為氣膜壓力；v為節(jié)流孔速度。

球坐標(biāo)系中連續(xù)性方程式為

式中：t為時間。

氣體狀態(tài)方程為

式中：pa為標(biāo)準(zhǔn)大氣壓；ρa(bǔ)為空氣密度。

節(jié)流孔的流量方程為

式中：mi為通過第i個小孔的質(zhì)量流量；?orifice為流量系數(shù)；Aorifice為小孔橫截面積；ψi為流量函數(shù)。

聯(lián)立（1）～（5）式可推導(dǎo)出球面動靜壓氣體軸承非線性量綱一的Reynolds方程為

式中：Q為節(jié)流孔引入的氣體質(zhì)量流量因子；R0為轉(zhuǎn)子半徑；P為量綱一的氣膜壓力；H為量綱一的氣膜厚度；T為量綱一的時間。

為了便于求解，對（6）式進(jìn)行保角變換［11］。變換形式取ξ＝ln［tan（α／2）］，將（6）式由三維求解域轉(zhuǎn)變成求解域為平面矩形的球面動靜壓氣體軸承潤滑分析方程［12］。氣膜的動態(tài)特性表現(xiàn)為轉(zhuǎn)子在外部隨機(jī)擾動下偏離靜平衡位置，在其附近作變位運動時氣膜力的相應(yīng)變化情況，所以應(yīng)以不定常工況的Reynolds方程作為分析計算的基礎(chǔ)［4］。將（6）式變?yōu)?/p>

式中：Ve為徑向偏心距方向的量綱一的速度；ve為沿徑向偏心方向的速度；Vθ為轉(zhuǎn)子量綱一的角速度；vθ為轉(zhuǎn)子中心繞軸承中心轉(zhuǎn)動的線速度；Vz為轉(zhuǎn)子軸向的量綱一的速度；vz為軸心沿軸向方向的速度。

3 動態(tài)特性系數(shù)求解

3.1 擾動壓力方程的推導(dǎo)

采用偏導(dǎo)數(shù)法［12］對動態(tài)控制方程（7）式在位移與速度方向上的擾動變量進(jìn)行求導(dǎo)，從而得到擾動壓力控制方程。求解域分為小孔區(qū)域（連續(xù)區(qū)域）、非小孔區(qū)域（包括連續(xù)與不連續(xù)區(qū)域）。以求解小孔區(qū)域為例，將（7）式對ε求偏導(dǎo)。

等式左邊第1項為

等式左邊第2項為

等式左邊第3項為

等式右邊第1項為

等式右邊第2項為

將5個偏導(dǎo)數(shù)項式相加可得

穩(wěn)態(tài)控制方程為

將（14）式變形得

將（15）式代入（13）式，化簡得出關(guān)于Pε的擾動壓力控制方程為

由于氣膜厚度不連續(xù)，分別計算臺區(qū)和槽區(qū)的氣膜厚度：

臺區(qū)量綱一的氣膜厚度為

槽區(qū)量綱一的氣膜厚度為

式中：Hb為量綱一的槽深；ε1為量綱一的軸向偏心率；e1為軸向偏心距；ε為量綱一的徑向偏心率；e為徑向偏心距。

對氣膜厚度方程求偏導(dǎo)，代入（16）式可得關(guān)于Pε的擾動壓力控制方程為

同理，對其余5個擾動變量求偏導(dǎo)，可得其他5個擾動壓力的控制方程。

3.2 數(shù)值計算擾動壓力

在廣義坐標(biāo)系下用有限差分法對（19）式離散化，推導(dǎo)出擾動壓力的差分表達(dá)式。為了計算的方便與精確，在斜坐標(biāo)系下（x＝?＋ξ／tanβ，y＝-ξ／sinβ）劃分網(wǎng)格，周向網(wǎng)格數(shù)為500（根據(jù)計算調(diào)整設(shè)定），徑向網(wǎng)格數(shù)為90；周向步長0.012 56，徑向步長0.019 23。網(wǎng)格劃分時周向槽臺邊界線到原點之間設(shè)定網(wǎng)格數(shù)50（根據(jù)計算調(diào)整設(shè)定），使該邊界線正好在網(wǎng)格點上；徑向槽臺邊界線計算時對網(wǎng)格數(shù)取整，使邊界線落在網(wǎng)格線上［6］，網(wǎng)格劃分如圖3所示。

圖3 求解域網(wǎng)格劃分Fig.3 Meshing for solution domain

以求解小孔區(qū)域為例，先對網(wǎng)格區(qū)域內(nèi)平行四邊形abcd區(qū)域（圖4）進(jìn)行面積分，再利用格林公式得到相應(yīng)區(qū)域的線積分，最后采用有限差分法對積分公式離散化處理，得到擾動壓力的差分表達(dá)式［4，12］。

圖4 氣膜厚度連續(xù)區(qū)域網(wǎng)格Fig.4 Continuous regional grid of gas film thickness

小孔區(qū)域關(guān)于ε擾動壓力控制方程的差分表達(dá)式為

式中：Ai，j，Bi，j，Ci，j，Di，j，Ei，j，F(xiàn)i，j是與氣膜厚度Hi，j及穩(wěn)態(tài)氣膜壓力Pi－1，j，Pi＋1，j，Pi，j－1，Pi，j＋1有關(guān)的系數(shù)。同理可以求出其他5個擾動壓力的差分表達(dá)式。

求解擾動控制方程的邊界條件為：1）大氣邊界條件為

2）邊界對稱條件為

3）小孔處邊界條件為

式中：Pd為節(jié)流孔后端壓力。

采用超松弛迭代法對（20）式求解可得擾動壓力Pε為

式中：λ為松弛因子，一般取大于1；k為迭代次數(shù)。

同理可以求出小孔區(qū)域其他5個擾動壓力。非小孔區(qū)域擾動壓力的求解與小孔區(qū)域類似，這里不再贅述。

3.3 動態(tài)剛度系數(shù)和動態(tài)阻尼系數(shù)的求解

在轉(zhuǎn)子穩(wěn)定位置給一個小擾動，轉(zhuǎn)子偏離穩(wěn)態(tài)平衡位置作變位運動，氣膜壓力發(fā)生改變。軸承在?，ξ，z方向偏離穩(wěn)態(tài)平衡位置O1時的受力分析如圖5所示。

圖5 軸承受力分析Fig.5 Load analysis of bearing

軸心在穩(wěn)態(tài)平衡位置的壓力分布為

軸心動態(tài)壓力分布為

將動態(tài)壓力分布p在穩(wěn)態(tài)平衡位置的周圍展開為軸心偏離平衡位置的瞬時位移和瞬時變位速度的Taylor級數(shù)（只考慮一階小項）［13］

式中：（pe，pθ，pz，pe，pθ，pz）分別為p隨（Δe，eΔθ，Δz，Δe，eΔθ，Δz）的變化率，即擾動壓力。

對動態(tài)壓力分布進(jìn)行積分，得到軸心作變位運動時各方向的氣膜力為

氣膜力是由氣膜壓力積分所得，通過氣膜力對各方向位移與速度求導(dǎo)可得氣膜的剛度和阻尼系數(shù)為

對于Fe，fe ＝cos?sin2α，

對于Fθ，fθ＝sin?sin2α，

對于Fz，fz ＝cosαsinα，

式中：Sj為某個方向的位移擾動。

采用C＋＋編程計算軸承的剛度和阻尼系數(shù)，流程圖如圖6所示。為了便于計算，采用承載力進(jìn)行判別，而不用氣膜力，保證轉(zhuǎn)子在穩(wěn)定位置附近發(fā)生小擾動。

圖6 軸承剛度和阻尼系數(shù)計算流程圖Fig.6 Calculation flow chart of stiffness and damping coefficients of bearing

4 動態(tài)特性系數(shù)計算與分析

氣體軸承參數(shù)見表1。

轉(zhuǎn)速n＝20 000 r／min，ps＝0.3 MPa下，不同供氣切向角和偏心率對剛度、阻尼系數(shù)的影響（氣膜剛度和阻尼系數(shù)的增大或減小均是相對其絕對值而言）如圖7所示。

由圖7a—圖7c可知：在同一偏心率下，剛度系數(shù)Kee，Kθe，Kθθ隨供氣切向角的增加而減小；供氣切向角不變時，除了Kee外，偏心率增大時，剛度系數(shù)增加；除了Kee，Kθe，Kθθ外，供氣切向角對剛度系數(shù)的影響較小。

由圖7d—圖7f可知：在相同供氣切向角下，除了Bee，Bθe外，偏心率增大時，動態(tài)阻尼系數(shù)增加；阻尼系數(shù)Bee，Bθe受偏心率的影響較弱；當(dāng)偏心率不變時，阻尼系數(shù)Bθθ，Bze，Bez，Bzz隨供氣切向角的增大而減??；除了Bθθ，Bez，Bze，Bzz外，供氣切向角對阻尼系數(shù)的影響較小。

圖7 供氣切向角和偏心率對剛度和阻尼系數(shù)的影響Fig.7 Influence of tangential angle of gas supply and eccentricity on stiffness and damping coefficient

n＝20 000 r／min，φ＝30°下，不同供氣壓力和偏心率對動態(tài)剛度和阻尼系數(shù)的影響如圖8所示。

由圖8a—圖8c可知：除了Kee外，偏心率增大時，剛度系數(shù)增加；除了Kee，Kθe，Kθθ，供氣壓力對量綱一的剛度系數(shù)的影響較?。辉谕黄穆氏?，剛度系數(shù)Kee，Kθe，Kθθ隨供氣壓力的增大而增大。

由圖8d—圖8f可知：供氣壓力不變時，除了Bee，Bθe外，偏心率增大時，阻尼系數(shù)增加；阻尼系數(shù)Bee，Bθe受偏心率的影響較弱；當(dāng)偏心率不變時，供氣壓力的增大，阻尼系數(shù)Bze，Bθθ，Bzz增加；除了Bze，Bθθ，Bzz外，供氣壓力對阻尼系數(shù)的影響較小。

圖8 供氣壓力和偏心率對剛度和阻尼系數(shù)的影響Fig.8 Influence of pressure of gas supply and eccentricity on stiffness and damping coefficient

5 結(jié)論

將半球面氣體軸承的動態(tài)剛度系數(shù)、阻尼系數(shù)與穩(wěn)定性關(guān)聯(lián)起來，分析結(jié)論如下：

1）除了Kee外，偏心率增大，剛度系數(shù)增加；除了Bee，Bθe外，偏心率增大，阻尼系數(shù)增加。

2）供氣切向角增大，動態(tài)剛度系數(shù)Kee，Kθe，Kθθ減小，阻尼系數(shù)Bθθ，Bez，Bze，Bzz減小。

3）供氣壓力增大時，剛度系數(shù)Kee，Kθe，Kθθ增加，阻尼系數(shù)Bze，Bθθ，Bzz增加。