基于流體域仿真的高速風(fēng)機(jī)氣浮箔片軸承動(dòng)靜特性研究

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

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

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

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

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

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

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

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

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

1 計(jì)算方法與流體域模型

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

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

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

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

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

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

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

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

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

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

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

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

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

從圖6中可以看出，在流體域內(nèi)靠近軸頸邊界的氣體速度較大，隨之速度沿徑向方向不斷下降，在箔片邊界附近的氣體速度最低，與實(shí)際情況相比較為符合。從圖7中可以看出，流體域內(nèi)氣體的壓力最大值發(fā)生在流體域厚度最小附近，該現(xiàn)象產(chǎn)生的原因可由流體動(dòng)壓潤滑理論解釋。流體動(dòng)壓潤滑是指軸頸在旋轉(zhuǎn)時(shí)，在軸頸與軸承的摩擦表面之間會(huì)充滿流體，由于流體自身具有黏性，當(dāng)軸頸達(dá)到足夠高的旋轉(zhuǎn)速度時(shí)，流體會(huì)被軸頸帶入軸頸與軸承表面之間的空間。對(duì)于文中所研究的氣浮軸承，即為將流體帶入軸頸與箔片之間的流體域內(nèi)，在軸頸與箔片之間形成楔形間隙，從而形成流體動(dòng)壓效應(yīng)；被帶入楔形間隙的流體會(huì)產(chǎn)生壓力，該壓力與軸頸賦予的載荷(即為外載荷)平衡時(shí)，軸頸與箔片之間形成穩(wěn)定的氣膜，從而實(shí)現(xiàn)流體的動(dòng)壓潤滑。因此，流體域內(nèi)氣體的壓力分布情況與理論分析結(jié)果相符合。如圖8所示，軸承的流體域內(nèi)變化較大的湍動(dòng)能主要發(fā)生在箔片搭接的位置，即為流體域內(nèi)氣體厚度較大的位置，因?yàn)樵撎帪椴罱?，流體域邊界發(fā)生突變，邊界的變化會(huì)導(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)速下流體域湍動(dòng)能分布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)速下氣體湍動(dòng)能分布Fig.11 Turbulent kinetic energy distribution at different speeds： (a) 80 000 r/min；(b)60 000 r/min

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

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

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

從圖12、13中可以看出，軸頸表面所受到的壓力隨著軸心位置的變化呈現(xiàn)周期性變化的規(guī)律。因此，需要通過數(shù)據(jù)處理，從壓力變化曲線中找到與軸心位置變化頻率相同的曲線；同時(shí)求解軸心位置變化曲線與壓力變化曲線之間的相位差，根據(jù)頻率與相位差值即可得到求解剛度系數(shù)與阻尼系數(shù)所需要的壓力變化表達(dá)式。使用相同的方法，賦予同一轉(zhuǎn)速2個(gè)不同的同頻擾動(dòng)，處理數(shù)據(jù)得到對(duì)應(yīng)的壓力變化表達(dá)式，將表達(dá)式代入方程組內(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ì)則是先下降而后上升；交叉剛度系數(shù)隨著轉(zhuǎn)速的提升先下降而后上升，X軸與Y軸的交叉阻尼系數(shù)的變化趨勢(shì)則相反。這是因?yàn)楦鶕?jù)動(dòng)壓潤滑理論與方程可知，任意一點(diǎn)的氣膜壓力與該點(diǎn)的速度梯度的導(dǎo)數(shù)有關(guān)。因此，對(duì)于某一固定位置，當(dāng)軸頸的轉(zhuǎn)速較低時(shí)，軸頸偏心率較大；隨著轉(zhuǎn)速的提升，軸頸偏心率會(huì)逐漸減小，在此過程中，該點(diǎn)的速度梯度對(duì)應(yīng)發(fā)生變化，從而產(chǎn)生的氣膜壓力也隨之發(fā)生改變，代入方程求解得到的剛度系數(shù)與阻尼系數(shù)呈現(xiàn)不斷下降的趨勢(shì)。隨著軸頸轉(zhuǎn)速的不斷提升，軸頸的偏心率變化會(huì)不斷變小，該位置的速度梯度會(huì)隨之增大，產(chǎn)生的氣膜壓力也不斷提高，從而會(huì)使得剛度系數(shù)與阻尼系數(shù)隨著轉(zhuǎn)速的不斷提升呈現(xiàn)前后不同的變化趨勢(shì)。

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

4 結(jié)論

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