張晨，許亮，孫丹，李雨薇，高飛，荊建平

(1.上海交通大學(xué)機械系統(tǒng)與振動國家重點實驗室，上海 200240；2.中國船舶集團有限公司第七〇三研究所， 黑龍江哈爾濱 150078；3.新鄉(xiāng)航空工業(yè)(集團)有限公司，河南新鄉(xiāng) 453000)

高速風機作為航空、艦船、化工等工業(yè)領(lǐng)域常見的流體機械，其制造產(chǎn)業(yè)的提升對工業(yè)流體機械的發(fā)展尤為重要。對于風機而言，隨著工作時間的增加，其動平衡會被逐漸破壞，軸承部位的振動會不斷加大。一旦軸承的振動大于某一限值，就會對風機造成巨大的影響，必須停機修理[1]。隨著科學(xué)技術(shù)水平的不斷提升，風機應(yīng)用的環(huán)境多變復(fù)雜，風機向著體積小、質(zhì)量輕的方向不斷發(fā)展，這對風機轉(zhuǎn)子系統(tǒng)所使用的軸承的精度提出了更高的要求。氣浮軸承憑借著精度高、發(fā)熱小、噪聲低等各項優(yōu)點，正在逐步替代傳統(tǒng)的機械軸承，成為風機等高精度設(shè)備中的關(guān)鍵部件[2]。如今，氣浮軸承應(yīng)用于風機的轉(zhuǎn)子系統(tǒng)中，針對風機類產(chǎn)品轉(zhuǎn)子系統(tǒng)特性的分析、臨界轉(zhuǎn)速的預(yù)測等研究便是所需要解決的問題。目前市面上常見的商業(yè)計算軟件均缺失較為系統(tǒng)的氣浮軸承特性計算方法，氣浮軸承的研究也在不斷地發(fā)展與完善。

工業(yè)中的高速風機多數(shù)選擇使用高速氣浮軸承，結(jié)構(gòu)如圖1所示，該類軸承的特性計算難度更大[3]。并且，為了提高氣浮軸承的抗沖擊能力與自適應(yīng)能力，在軸與軸承之間加入了具有柔性的金屬箔片(見圖2)，其能夠根據(jù)轉(zhuǎn)速與載荷的變化建立起合適的工作氣膜厚度，同時箔片軸承的結(jié)構(gòu)阻尼效應(yīng)能夠更好地抑制轉(zhuǎn)子的振動，從而有更好的穩(wěn)定性。因此對多葉箔片軸承的特性進行理論分析有更加重要的現(xiàn)實意義[4]。

圖2 氣浮軸承箔片F(xiàn)ig.2 Foil for air bearing

WALOWIT等[5-6]是較早研究波箔型箔片軸承的學(xué)者，他們將起支撐作用的箔片從軸承中分離，分析單個箔片的受力情況，建立對應(yīng)的軸承模型，求解得到箔片軸承的靜態(tài)特性。以HESHMAT為代表的一批學(xué)者對波箔型箔片軸承進行了比較深入的理論和實驗研究，采用有限差分法來計算彈性流體動力潤滑(EHD)問題，形成研究波箔型箔片軸承剛度和阻尼特性的經(jīng)典理論方法[7-9]。

在國內(nèi)，龔煥孫[10]最早開展波箔型箔片軸承剛度和阻尼特性的研究，他基于Krichhoff假設(shè)，施加載荷與邊界條件進行計算，得到軸承靜態(tài)特性，同時搭建實驗臺，對其靜態(tài)特性進行試驗研究，與理論分析的結(jié)果相符。宋來運[11]從理論和實驗驗證2個層面對動靜壓氣浮軸承轉(zhuǎn)子系統(tǒng)的靜態(tài)性能、動態(tài)性能和動態(tài)穩(wěn)定性等綜合性能進行了系統(tǒng)性的研究，為設(shè)計研發(fā)高速高精度動靜壓氣浮主軸系統(tǒng)提供理論基礎(chǔ)和技術(shù)支持。劉佳琦等[12]從彈性箔片軸承的結(jié)構(gòu)特性、涂層特性和建模方法3個方面對氣體彈性箔片軸承的研究進展進行綜述，展望彈性箔片軸承在高精尖等領(lǐng)域的前景，并對彈性箔片軸承的研究方向提出了建議。劉恒等人[13]對箔片含有預(yù)變形的多葉箔片氣體動壓軸承進行研究，將箔片簡化為彎曲梁，推導(dǎo)出多葉箔片的變形方程和考慮預(yù)變形后的氣膜厚度方程，計算研究轉(zhuǎn)速等參數(shù)對軸承承載力的影響。

然而，當前的氣浮軸承的特性計算多是基于一定的假設(shè)條件，通過簡化軸承特性方程來得到軸承特性。郭軍剛等[14]建立多葉油潤滑箔片軸承Reynolds方程、對應(yīng)的邊界條件以及收斂條件，通過傳遞矩陣計算箔片軸承軸系臨界轉(zhuǎn)速，同時，建立對應(yīng)的試驗系統(tǒng)進行評估與驗證。

由于眾多假設(shè)條件的存在，計算結(jié)果往往與實際結(jié)果存在偏差。因此，為了提高高速風機轉(zhuǎn)子動力學(xué)分析精度，減少理論計算與實際情況之間的偏差，本文作者針對現(xiàn)有的一類多箔片氣浮軸承，建立流體域模型，通過使用商業(yè)仿真軟件分析風機產(chǎn)品內(nèi)高速氣浮箔片軸承的靜態(tài)特性與動態(tài)特性隨著轉(zhuǎn)速提升的變化規(guī)律，從而更加精準地預(yù)估風機模態(tài)、臨界轉(zhuǎn)速以及失穩(wěn)轉(zhuǎn)速。在設(shè)計使用氣浮軸承的轉(zhuǎn)子系統(tǒng)時，依據(jù)氣浮軸承自身的動態(tài)特性變化趨勢來設(shè)計對應(yīng)的轉(zhuǎn)子系統(tǒng)結(jié)構(gòu)，從而能夠更大程度地提高轉(zhuǎn)子系統(tǒng)的穩(wěn)定性與安全性，這是開展軸承作為支承部件的轉(zhuǎn)子系統(tǒng)動力學(xué)特性研究的基礎(chǔ)。

1 計算方法與流體域模型

針對使用氣浮軸承的風機轉(zhuǎn)子系統(tǒng)，氣浮軸承中的多葉箔片在轉(zhuǎn)子系統(tǒng)啟動時，起到支撐作用；隨著軸頸轉(zhuǎn)速的不斷提升，帶動氣浮軸承的流體域內(nèi)的流體高速流動，由收斂楔進入發(fā)散楔，產(chǎn)生氣膜；氣膜可視作具有較大剛度的支撐體，具有足夠的承載能力，隨著氣膜的動壓力不斷提升，將軸頸慢慢托起，與箔片中間形成空隙，不產(chǎn)生接觸。當軸頸的載荷與氣膜所能夠提供的承載能力實現(xiàn)相互平衡時，即為氣浮軸承的工作狀態(tài)。

文中提出的軸承特性求解的方法原理如圖3所示，主要過程為：首先建立氣浮軸承的流體域模型，針對流體域模型進行網(wǎng)格劃分，通過商用仿真軟件模擬一定轉(zhuǎn)速下氣浮軸承流體域的分布情況；移動軸頸使得軸頸受到流體域的作用力與所受到的載荷相互平衡，從而得到氣浮軸承在一定載荷下的軸心軌跡圖，即為氣浮軸承的靜特性計算；在靜特性的計算結(jié)果的基礎(chǔ)上，賦予軸頸與氣浮軸承轉(zhuǎn)速相同的同頻擾動，讀取在擾動的作用下流體域?qū)S頸的作用力，并沿軸頸表面進行積分，得到軸頸沿X、Y軸方向的受力變化，從而計算氣浮軸承的剛度與阻尼系數(shù)，即為氣浮軸承的動特性計算。

圖3 軸承特性求解流程Fig.3 Flow for solving bearing characteristics

建立氣浮軸承的流體域模型時，為了更加便于仿真計算與后續(xù)的模型修正，基于流體域沿軸向分布變化較小的情況，建立流體域二維平面模型。使用該模型進行模擬仿真，得到流體域結(jié)果后，通過使用修正系數(shù)，將二維模型計算結(jié)果進行修正，得到三維流體域模型的結(jié)果。

文中研究的氣浮軸承流體域內(nèi)的流體為可壓縮氣體，使用商業(yè)仿真軟件根據(jù)氣體的連續(xù)方程與動量方程等進行求解，通過仿真模擬計算得到氣浮軸承流體域內(nèi)的氣體分布情況以及軸頸表面的壓力分布，列出軸頸的受力方程，從而可以求解得到8個動特性系數(shù)[11]。

在轉(zhuǎn)子軸心附近產(chǎn)生一個小擾動Δx、Δy，流體域?qū)D(zhuǎn)子的作用力ΔFx、ΔFy為

(1)

將位移與受力取諧波分量可得：

(2)

(3)

將諧波分量代入式(1)，可得：

(4)

求解該方程組，即可得到氣浮軸承的剛度系數(shù)與阻尼系數(shù)。

2 模型網(wǎng)格劃分與邊界條件設(shè)定

2.1 模型網(wǎng)格劃分

由于氣浮軸承結(jié)構(gòu)特殊性，為滿足CFD仿真模擬計算得到的流體域結(jié)果的精確度同時提高計算過程效率，選擇使用流體域二維模型進行仿真計算[11]，同時需要對流體域進行高質(zhì)量的網(wǎng)格劃分，得到更合適的網(wǎng)格結(jié)構(gòu)。網(wǎng)格分為結(jié)構(gòu)化網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格，因為該類氣浮軸承的流體域模型比較簡單，且比較規(guī)則，沒有太多復(fù)雜的流體區(qū)域，因此，選擇劃分結(jié)構(gòu)化網(wǎng)格。該類網(wǎng)格生成速度比較快，并且數(shù)據(jù)結(jié)構(gòu)簡單，可以滿足氣浮軸承流體域的仿真計算。網(wǎng)格單元越小的情況下，得到的計算結(jié)果精度更高，但同時會增加仿真模擬的計算時間；且在仿真模擬計算的過程中，需要控制軸頸移動，導(dǎo)致流體域隨之改變，涉及動網(wǎng)格設(shè)置問題，對網(wǎng)格的質(zhì)量要求也隨之提高。因此，需要選取滿足計算精度的網(wǎng)格尺寸與網(wǎng)格數(shù)量。文中采用的結(jié)構(gòu)化網(wǎng)格劃分流體域模型(見圖4)，經(jīng)過網(wǎng)格無關(guān)性驗證，采用網(wǎng)格數(shù)量分別為51 090、239 820、420 000與802 670的流體域模型。在該網(wǎng)格數(shù)下分別移動軸頸至工作轉(zhuǎn)速下的位置，得到的誤差如圖5所示?？梢钥闯?，節(jié)點數(shù)為399 800，網(wǎng)格數(shù)量為420 000的網(wǎng)格劃分得到的模型即可滿足仿真計算求解的精度要求。質(zhì)量較高的流體域模型網(wǎng)格為后續(xù)的氣浮軸承靜特性與動特性的計算求解奠定了基礎(chǔ)。

圖4 氣浮軸承流體域模型及網(wǎng)格劃分示意Fig.4 Fluid domain model and grid generation of air bearing

圖5 網(wǎng)格無關(guān)性驗證Fig.5 Grid independence verification

2.2 模型邊界條件設(shè)定

使用該類氣浮軸承的風機轉(zhuǎn)子系統(tǒng)的工作轉(zhuǎn)速為100 000 r/min，軸頸所受到的載荷為轉(zhuǎn)子自身的重量3 N，流體域內(nèi)的材料為可壓縮氣體，氣體的密度為8.13 kg/m3，氣體的壓縮系數(shù)為1.4，在仿真模擬的計算過程中不考慮氣體黏性的影響，氣體的壓力為800 kPa。模擬仿真的求解選擇Pressure-Based，Velocity Formulation選擇Absolute，Time選擇Transient，2D Space選擇Planar。

模型的設(shè)置中，Energy選擇On，Viscous選擇SSTk-ε，該模型可以較好地處理求解邊界層網(wǎng)格數(shù)據(jù)，在近壁面的位置能夠更好地進行處理。因為氣浮軸承的流體域模型的特點，可以看出，k-ε模型能夠更高精度地處理求解流體域內(nèi)的流動情況；同時，選用SST，相較于BSL，依然保留了在近壁面采用k-ε模型，但SST簡化了模型，減小計算量，對壁面的距離更加敏感，能夠進行更加精確地仿真求解；同時，壁面的設(shè)置中加入軸頸的轉(zhuǎn)速與流體域的壓力邊界條件，以及考慮移動軸頸的動網(wǎng)格，根據(jù)軸頸位置的改變來調(diào)整網(wǎng)格結(jié)構(gòu)。最后，根據(jù)網(wǎng)格的尺寸與軸頸的移動設(shè)定合適的時間步長與時間步數(shù)來得到更加精確的結(jié)果。得到靜特性的計算結(jié)果之后，賦予軸頸同頻擾動，求解得到軸頸表面的壓力分布，對軸頸表面進行積分得到軸頸沿X、Y軸方向受到的作用力，從而求解得到氣浮軸承的剛度系數(shù)與阻尼系數(shù)。完成100 000 r/min工作轉(zhuǎn)速下的氣浮軸承的特性計算之后，繼續(xù)求解60 000、80 000 r/min轉(zhuǎn)速下的氣浮軸承的特性，從而得到氣浮軸承從開始啟動到轉(zhuǎn)速提升至工作轉(zhuǎn)速這一過程中軸心位置的變化以及氣浮軸承的靜特性與動特性的變化趨勢。

3 仿真結(jié)果及分析

對氣浮軸承流體域模型進行計算流體力學(xué)求解并進行仿真計算，可以得到在不同的轉(zhuǎn)速下，氣浮軸承流體域內(nèi)的氣體速度分布、湍動能分布以及氣體壓力分布情況。圖6、7、8所示為100 000 r/min轉(zhuǎn)速下氣浮軸承流體域內(nèi)氣體的速度、壓力分布情況以及流體域內(nèi)的湍動能分布。

從圖6中可以看出，在流體域內(nèi)靠近軸頸邊界的氣體速度較大，隨之速度沿徑向方向不斷下降，在箔片邊界附近的氣體速度最低，與實際情況相比較為符合。從圖7中可以看出，流體域內(nèi)氣體的壓力最大值發(fā)生在流體域厚度最小附近，該現(xiàn)象產(chǎn)生的原因可由流體動壓潤滑理論解釋。流體動壓潤滑是指軸頸在旋轉(zhuǎn)時，在軸頸與軸承的摩擦表面之間會充滿流體，由于流體自身具有黏性，當軸頸達到足夠高的旋轉(zhuǎn)速度時，流體會被軸頸帶入軸頸與軸承表面之間的空間。對于文中所研究的氣浮軸承，即為將流體帶入軸頸與箔片之間的流體域內(nèi)，在軸頸與箔片之間形成楔形間隙，從而形成流體動壓效應(yīng)；被帶入楔形間隙的流體會產(chǎn)生壓力，該壓力與軸頸賦予的載荷(即為外載荷)平衡時，軸頸與箔片之間形成穩(wěn)定的氣膜，從而實現(xiàn)流體的動壓潤滑。因此，流體域內(nèi)氣體的壓力分布情況與理論分析結(jié)果相符合。如圖8所示，軸承的流體域內(nèi)變化較大的湍動能主要發(fā)生在箔片搭接的位置，即為流體域內(nèi)氣體厚度較大的位置，因為該處為箔片搭接，流體域邊界發(fā)生突變，邊界的變化會導(dǎo)致流體域厚度隨之改變，從而在此處更易發(fā)生湍流。

圖6 100 000 r/min轉(zhuǎn)速下流體域速度分布Fig.6 Velocity distribution in fluid domain at speed of 100 000 r/min

圖7 100 000 r/min轉(zhuǎn)速下流體域壓力分布Fig.7 Fluid domain pressure distribution at speed of 100 000 r/min

圖8 100 000 r/min轉(zhuǎn)速下流體域湍動能分布Fig.8 Turbulent kinetic energy distribution in fluid domain at speed of 100 000 r/min

圖9—11所示為氣浮軸承在80 000和60 000 r/min轉(zhuǎn)速下流體域內(nèi)的氣體速度分布、湍流分布以及氣體壓力分布情況。

圖9 不同轉(zhuǎn)速下氣體速度分布Fig.9 Speed distribution of gas at different speeds：(a) 80 000 r/min；(b)60 000 r/min

圖10 不同轉(zhuǎn)速下氣體壓力分布Fig.10 Pressure distribution of gas at different speeds： (a) 80 000 r/min；(b)60 000 r/min

圖11 不同轉(zhuǎn)速下氣體湍動能分布Fig.11 Turbulent kinetic energy distribution at different speeds： (a) 80 000 r/min；(b)60 000 r/min

從圖9—11中可以看出不同轉(zhuǎn)速下，速度分布、壓力分布以及湍動能分布情況以及變化趨勢都大致相同，不同之處在于轉(zhuǎn)速的變化，軸心位置會隨之改變，承載相同壓力的氣膜形狀發(fā)生變化，但氣膜對軸頸的壓力合力仍然保持與軸承的載荷相平衡，從而保證氣膜能夠穩(wěn)定存在，軸頸能夠正常運轉(zhuǎn)。

根據(jù)賦予軸頸同頻擾動后求解得到軸頸表面的壓力分布情況，對軸頸表面進行積分得到軸頸沿X、Y軸方向受到的作用力，圖12與圖13所示為賦予軸頸同頻擾動后，軸頸X、Y軸方向的位置變化與軸頸所受到的2個方向的壓力變化情況。

圖13 軸頸表面壓力變化Fig.13 Change of journal surface pressure

從圖12、13中可以看出，軸頸表面所受到的壓力隨著軸心位置的變化呈現(xiàn)周期性變化的規(guī)律。因此，需要通過數(shù)據(jù)處理，從壓力變化曲線中找到與軸心位置變化頻率相同的曲線；同時求解軸心位置變化曲線與壓力變化曲線之間的相位差，根據(jù)頻率與相位差值即可得到求解剛度系數(shù)與阻尼系數(shù)所需要的壓力變化表達式。使用相同的方法，賦予同一轉(zhuǎn)速2個不同的同頻擾動，處理數(shù)據(jù)得到對應(yīng)的壓力變化表達式，將表達式代入方程組內(nèi)，從而求解得到氣浮軸承的剛度系數(shù)與阻尼系數(shù)。圖14與圖15所示為剛度系數(shù)與阻尼系數(shù)隨著轉(zhuǎn)速提升的變化情況。

圖15 阻尼系數(shù)隨轉(zhuǎn)速的變化Fig.15 Variation of damping coefficient with speed

從圖14、15中可以看出，X軸與Y軸方向的直接剛度系數(shù)隨著轉(zhuǎn)速的提升不斷下降，直接阻尼系數(shù)的變化趨勢則是先下降而后上升；交叉剛度系數(shù)隨著轉(zhuǎn)速的提升先下降而后上升，X軸與Y軸的交叉阻尼系數(shù)的變化趨勢則相反。這是因為根據(jù)動壓潤滑理論與方程可知，任意一點的氣膜壓力與該點的速度梯度的導(dǎo)數(shù)有關(guān)。因此，對于某一固定位置，當軸頸的轉(zhuǎn)速較低時，軸頸偏心率較大；隨著轉(zhuǎn)速的提升，軸頸偏心率會逐漸減小，在此過程中，該點的速度梯度對應(yīng)發(fā)生變化，從而產(chǎn)生的氣膜壓力也隨之發(fā)生改變，代入方程求解得到的剛度系數(shù)與阻尼系數(shù)呈現(xiàn)不斷下降的趨勢。隨著軸頸轉(zhuǎn)速的不斷提升，軸頸的偏心率變化會不斷變小，該位置的速度梯度會隨之增大，產(chǎn)生的氣膜壓力也不斷提高，從而會使得剛度系數(shù)與阻尼系數(shù)隨著轉(zhuǎn)速的不斷提升呈現(xiàn)前后不同的變化趨勢。

圖14 剛度系數(shù)隨轉(zhuǎn)速的變化Fig.14 Variation of stiffness coefficient with speed

4 結(jié)論

建立了高速氣浮軸承流體域模型，通過仿真計算方法求解了軸承特性，分析了氣浮軸承工作狀態(tài)下穩(wěn)定時的軸心位置隨著轉(zhuǎn)速提升對應(yīng)的變化情況，以及氣浮軸承的動態(tài)特性系數(shù)與轉(zhuǎn)速的關(guān)系，發(fā)現(xiàn)隨著轉(zhuǎn)速的不斷提升，穩(wěn)定時的軸心位置會逐漸向氣浮箔片軸承的圓心靠近；同時，軸承的直接剛度系數(shù)在一定的工作轉(zhuǎn)速范圍內(nèi)隨著轉(zhuǎn)速的提升而下降，直接阻尼系數(shù)則在工作轉(zhuǎn)速范圍內(nèi)下降到達低值后上升。研究結(jié)果可為預(yù)估風機的失穩(wěn)轉(zhuǎn)速與臨界轉(zhuǎn)速提供理論支撐。