基于流固耦合的箔片氣體軸承動態特性分析

摘要

基于有限元軟件ANSYS Workbench，建立箔片氣體軸承在可壓縮流體介質中運動的有限元模型，采用6DOF動網格計算方法對軸承的運動狀態進行流固耦合數值模擬，探討了不同轉速和波箔片結構參數（波箔片的長度比、高度以及厚度）對軸承動態特性的影響規律。仿真結果表明：轉速增加，軸承的承載能力增加，但穩定性有所下降，更容易發生失穩現象；選取波箔片的長度比在1～1.5之間、厚度為0.16 mm，既可以保證軸承具有較高的剛度，同時又能獲得較大的阻尼；波箔片高度與軸承動態特性成反比關系。將仿真結果與試驗結果進行對比，驗證了仿真計算方法的正確性和有效性。同時本文的研究為優化波箔片結構，改善軸承動態特性，提高軸承運行穩定性提供理論依據。

關鍵詞

箔片氣體軸承; 流固耦合; 波箔片結構參數; 動態特性; 仿真

引 言

箔片氣體軸承是目前研究和使用最多的動壓軸承，其具有顯著的性能優勢，例如：自作用，安全運行，無外部加壓要求。由于其以空氣作為潤滑劑從而減少了潤滑裝置，且節省了維護成本，因此在航空、航天等領域得到了廣泛應用。隨著研究的不斷深入，其在醫療器械和船舶等領域的應用也有重大突破［1?2］。

李昊等［3］采用有限差分法和小擾動法，分析了箔片氣體軸承中箔片安裝位置對軸承的靜/動態特性的影響，認為安裝位置存在一“敏感”區間，在該區間內軸承的穩定性會受到嚴重影響，且該位置區間在數值上與軸承的偏位角近似互補，因此箔片安裝位置在偏位角對側時，能有效提高軸承的穩定性。梁波等［4］應用有限元理論提出了小柔度步計算方法，分析了不同平箔片厚度和不同波箔片彈性對軸承承載性能的影響，結果表明，波寬是影響軸承承載性能的重要參數，過小的平箔片厚度會引起氣膜振蕩。Zhao等［5］基于有限元理論使用梁單元將平箔片和波箔片進行離散，建立了考慮摩擦接觸情況下的箔片結構二維非線性模型，發現摩擦會導致箔片結構的變形發生遲滯效應，且會影響軸承能量的耗散。隨著研究的深入，馮凱等［6］利用有限元法對新型三瓣式箔片軸承展開研究，分析在接觸面間庫侖摩擦力條件下預載和安裝角度對軸承靜/動態規律的影響，發現安裝角度對軸承的靜/動態特性的影響呈正/余弦變化，同時發現預載可以提高軸承的負載能力。Xu等［7］在可壓縮氣體的雷諾方程的基礎上，提出了流體結構相互作用的算法，分析了箔片氣體軸承的靜/動態特性，認為波箔片厚度與軸承的承載力成反比，與偏位角成正比，而波箔片間距恰恰相反，且平箔片厚度的增加會導致動態剛度、阻尼系數增加。吳垚等［8］利用MATLAB中的偏微分方程工具箱，探討了軸承間隙、氣體稀薄效應和軸承參數對軸承動態特性的影響，結果表明，軸承間隙在納米級別時，軸承參數對動態特性影響較小，間隙減小、氣體稀薄程度增加會導致動態剛度系數減小，阻尼系數增加。目前的理論研究一般以雷諾方程為基礎，應用有限差分法和有限元法等方法將箔片軸承的波箔片簡化為線性彈簧或梁單元等形式以構建箔片軸承的彈性模型，高估了波拱柔度，且在計算過程中認為波箔片與平箔片、軸承座的接觸方式是線接觸，并忽略了波箔片的周向變形，因此計算結果與實際情況存在一定的誤差。隨著計算流體力學理論的逐步完善和相關商業計算軟件的不斷成熟，一些學者成功采用ANSYS Workbench中的流固耦合仿真方法和6DOF動網格計算方法［9］模擬了研究對象的材料屬性、接觸和變形情況［10?11］，并取得了很好的研究成果，而國內少有采用該方法對箔片氣體軸承進行的研究。

基于上述情況，本文以箔片氣體軸承為研究對象，將流固耦合理論與流體潤滑理論相結合，基于ANSYS Workbench仿真平臺，借助流體求解器中的6DOF動網格計算方法，對氣膜流場和彈性箔片間的相互作用進行流固耦合仿真。計算得到軸承的軸徑氣膜力、位移、速度等參數，進而對軸承的動態特性系數進行數值求解，分析軸承在不同轉速和波箔片結構參數條件下的動態特性規律，為該軸承的設計與應用提供理論支撐。

1 箔片氣體軸承的結構

箔片氣體軸承由平箔片、波箔片與軸承殼體三部分構成。其中頂部光滑的平箔片形成軸承的內表面，與轉子一同形成楔形間隙，為動壓氣膜的形成提供條件，該平箔片靠在由波紋狀波箔制成的彈性支撐結構上。波箔片上每個凸起的波箔塊都相當于一個小的彈簧，共同形成支撐結構，該結構為軸承提供結構剛度與阻尼，并會影響軸承的流體力學性能和最大負載能力。平箔片與波箔片的一端通過插槽固定在軸承殼體上，另一端自由延伸。

圖1（a）為軸承結構示意圖，圖1（b）為波箔片結構示意圖，其中，hb為波箔片高度，tb為波箔片厚度，lo為波箔片半波長度，ls為兩凸起波箔塊之間的距離。箔片氣體軸承的結構參數如表1所示。在后期對多種波箔片設計參數的箔片氣體軸承動態特性進行仿真分析時，為得出長度比ls/lo、波箔片高度hb、波箔片厚度tb與動態特性的關系，采用控制變量法將參數設置為ls/lo=1/2，hb=0.5 mm，tb

=0.1 mm。

2 仿真計算

2.1 網格無關性驗證

在有限元仿真計算中，網格的疏密程度會直接影響到計算結果的精確性。如果網格劃分得很細密，雖然可以提高計算結果的精確度，但也會大大增加計算時間，而且當網格達到一定數量時，計算精度并不會明顯提高，因此需要做網格無關性驗證。本文以穩定后流體域輸出的最大壓力作為評價標準，從圖2中可以看出，當流體域網格數量達到29萬后，流體域輸出最大壓力隨網格數量的變化基本不再改變，故流體域的網格劃分數量為294196。因流體域外表面與平箔片內表面是流固耦合面，這兩個面間的網格大小應盡可能一致，所以固體域的網格劃分數量為50888。流體域和固體域網格劃分如圖3所示。

2.2 邊界條件

設定符合實際情況的邊界條件對準確求解實際問題來說是極其重要的。邊界條件的選定要充分考慮模型所處的環境狀況、流體的類型以及相關控制方程等。針對本文研究的箔片氣體軸承，各邊界條件的設定如下：

（1）氣膜內表面即轉子外表面設置為旋轉壁面，用來模擬轉子的旋轉運動。

（2）氣膜外表面設置為無滑移壁面，其與平箔片的內表面構成流固耦合交界面。

（3）氣膜兩端面一端設置為壓力入口邊界條件，另一端設置為壓力出口邊界條件，壓力數值大小為當地大氣壓值。

2.3 動態特性仿真計算

軸承的動態仿真計算是箔片氣體軸承在特定的運行參數和結構參數條件下，通過流體求解器仿真求解楔形間隙內的流體壓力分布，進而得到軸徑的氣膜力、位移、速度值。軸承的動態特性計算流程圖如圖4所示，其中Fluent（流體求解器）中的動網格可以將計算過程中發生變形的網格進行重構，使網格質量始終保持在較高水平，而有限差分法和有限元法很難達到這種條件；通訊模塊System Coupling（SC）可以實現Fluent與Transient Structural（固體求解器）間的數據傳遞：Fluent先依據設定好的計算模型、求解算法與邊界條件等對氣膜模型進行仿真求解，求得的氣膜壓力分布與數值大小通過SC傳遞給Transient Structural作為結構分析的外載荷，待Transient Structural將箔片的變形結果求出后，再次通過SC將結構的位移數據傳遞給Fluent，作為流體域新的幾何邊界，如此反復迭代，直到耦合參與者均收斂。

首先編寫6DOF動態程序，6DOF的含義為轉子的6個自由度，即三個方向的旋轉自由度和平移自由度；由于試驗中采用了止推軸承，抵消了軸向的平移運動和x，y方向的旋轉運動，因此6DOF程序中應限制z方向的平移自由度和x，y方向的旋轉自由度，程序中還應包含轉子的質量和轉動慣量。同時還需編寫DEFINE_EXECUTE_AT_END宏命令，此宏命令會對求解過程中每個時間步后的軸徑分力進行提取備用；其次利用動網格方法將仿真求解的軸徑軌跡數據保存到本地備用；然后利用流體求解器中的監視器，對求解過程中每個時間步后的軸徑分速度進行提取備用；最后依據時間的先后順序，把軸徑沿x，y方向的氣膜力、速度和位移依次保存到EXCEL中備用。

由于時間步長非常短，連續5次輸出的數據中轉子運行狀態和氣膜力變化非常小，因此計算力、速度和位移的變化量時從第6組開始，然后由下式計算軸承的動態剛度、阻尼系數：

式中 ΔFx，ΔFy分別表示軸承在x，y方向的氣膜力增量；Kxx，Kyy，Kxy，Kyx和Cxx，Cyy，Cxy，Cyx分別表示軸承的動態剛度系數和阻尼系數；Δx，Δy和Δx˙，Δy˙分別表示軸承在x，y方向的位移增量和速度增量。

據此編寫MATLAB程序，利用滾動迭代法進行數據處理，并對所求結果取平均值得到箔片氣體軸承的動態特性系數。

3 動態特性仿真結果分析

3.1 轉速對軸承動態特性系數的影響

圖5為箔片氣體軸承在不同轉速下的動態剛度系數、阻尼系數的變化規律圖。由于動態特性系數的正負只代表與假定方向是否相同，與數值大小無關，所以本文中對動態特性系數的描述均是相對于其絕對值而言。由圖5（a）可以看出：主剛度系數Kxx，Kyy變化趨勢一致，均隨轉速的提高而增大，且主剛度系數要遠大于交叉剛度系數；而交叉剛度系數Kxy，Kyx隨轉速的變化則相對較為平緩。可以看出：主剛度系數主要起承載作用，隨著轉速的增加，軸承的承載能力得到提升。

由圖5（b）可以看出：隨轉速的提高，主阻尼系數Cxx，Cyy顯著減小，交叉阻尼系數Cxy，Cyx則逐漸減小并趨于0，可見動態阻尼系數均隨著轉速的提高而減小。可以看出：主阻尼系數主要起抑制渦動的作用，隨著轉速的增加，軸承穩定性下降，更易發生軸承失穩現象。這與文獻［12］的結論和規律基本一致。

3.2 長度比ls/lo對軸承動態特性系數的影響

圖6為軸承在轉速為40000 r/min時，不同長度比ls/lo下的動態剛度系數、阻尼系數的變化規律圖。由圖6（a）可以看出：主剛度系數Kxx，Kyy與交叉剛度系數Kxy，Kyx的變化趨勢一致，均隨著ls/lo

的增加而增大，但ls/lo對主剛度系數的影響要大于交叉剛度系數。由于箔片氣體軸承的特殊結構致使其動態特性由氣、固特性兩部分共同決定［1］。當軸徑的轉速和偏心率一定時，軸承的動態特性主要受箔片材料和結構的影響，當增加波箔片長度比ls/lo時，波箔片的結構剛度會增大，從而提高了波箔結構抵抗擾動的能力，軸承的動態剛度系數也隨著波箔片長度比ls/lo的增加而增大。

由圖6（b）可以看出：主阻尼系數Cxx，Cyy的變化趨勢一致，均隨著波箔片長度比ls/lo的增加而增大，而交叉阻尼系數Cxy，Cyx的變化并不明顯。軸承的動態阻尼特性一部分來自于箔片結構，而箔片的結構阻尼主要來自于平箔片與波箔片之間、波箔片與軸承座之間的庫侖摩擦力，波箔片長度比ls/lo

的改變并未對庫侖摩擦力造成影響，所以動態阻尼系數變化相對并不明顯。綜上所述，依據波箔片長度比ls/lo對軸承動態特性系數的影響規律，在波箔片的幾何形狀設計中應盡可能使波箔片長度比ls/lo在1～1.5之間，這樣既可以保證軸承具有較高的剛度，同時又能獲得較大的阻尼。

3.3 波箔片高度hb對軸承動態特性系數的影響

圖7為軸承在轉速為40000 r/min時，不同波箔片高度hb下的動態剛度系數、阻尼系數的變化規律圖。由圖7（a）與（b）可以看出：波箔片高度hb與軸承動態剛度系數、阻尼系數呈反比關系，從原理上分析，波箔片高度增加時，箔片就很容易發生變形，因此剛度系數、阻尼系數將會減小。在制造加工時，波箔片過低會給箔片成型帶來困難，所以在波箔片的幾何形狀設計中應盡可能使波箔片高度hb

的取值在0.5～0.6 mm之間。

3.4 波箔片厚度tb對軸承動態特性系數的影響

圖8為軸承在轉速為40000 r/min時，不同波箔片厚度tb下的動態剛度系數、阻尼系數的變化規律圖。由圖8（a）可以看出：主剛度系數Kxx，Kyy與交叉剛度系數Kxy，Kyx的變化趨勢一致，均隨著波箔片厚度tb的增加而增大，但其對主剛度系數的影響要大于交叉剛度系數。當波箔片厚度超過0.16 mm時，剛度系數增大趨勢變小。

由圖8（b）可以看出：主阻尼系數Cxx，Cyy的變化趨勢一致，均隨著波箔片厚度tb的增加而增大，而交叉阻尼系數Cxy，Cyx的變化并不明顯。當波箔片厚度超過0.16 mm時，阻尼系數變化緩慢。可見波箔片厚度tb對軸承動態剛度系數、阻尼系數的影響也是通過改變波箔片的結構剛度實現的。綜上所述，依據波箔片厚度tb對軸承動態特性系數的影響規律，在波箔片的幾何形狀設計中應盡可能使波箔片厚度tb的取值接近0.16 mm。

4 試驗分析與驗證

4.1 試驗原理

箔片動壓氣體軸承測試試驗臺由高速電主軸系統1（1.1～1.7）、被測軸承2、加載裝置3、數據采集與處理系統4（4.1～4.6）四部分組成，圖9所示為箔片動壓氣體軸承試驗測試分析系統原理圖。試驗機安裝在光學隔振平臺上，高速電主軸驅動系統為箔片動壓氣體軸承試驗臺的動力裝置，由高速電主軸、電主軸支座和刀具短軸組成，其最高轉速可達100000 r/min。被測軸承安裝在刀具短軸中間位置，被測軸承上端與加載裝置相連接，由加載裝置來實現載荷大小的改變。數據采集系統主要由電渦流位移傳感器、光電傳感器和溫度傳感器等采集軸承x，y方向位移、轉軸轉速并存入數據庫，數據處理系統用于實時監控與處理試驗數據以得到所需的結果。圖10為測試試驗臺實物圖。

4.2 試驗數據處理

氣體軸承轉子運行時，軸承偏離穩態平衡位置，并在其周圍進行變位運動時氣膜力的變化情況，氣膜力增量的位移擾動和速度擾動關系表達式指軸承在運行時，軸承偏離穩態平衡位置，并在其周圍進行變位運動時氣膜力的變化情況。結合氣體軸承動力學方程，計算軸承的剛度系數和阻尼系數。

當軸承軸心在靜平衡位置做微小運動時，氣膜力的增量由如下線性關系表示：

式中 m表示轉子的質量。

聯立式（2）和（3）可得到軸承的剛度、阻尼模型方程：

求解軸承剛度、阻尼系數，模型方程（4）中的每一個方程必須由4個不同的等式構成，且其系數矩陣的行列式不為0。試驗過程中，數據采集系統采樣頻率非常高，近似認為4次采樣過程中軸承只進行微小運動，軸承運行狀態和氣膜力變化微小，可忽略不計，因此將數據庫中第1～4次采樣預處理數據代入下式所示模型矩陣方程中進行計算，其中，n表示第n次采集：

可得第4次采樣時軸承的剛度、阻尼系數；按照上述方法編寫軸承剛度系數、阻尼MATLAB求解程序，使用滾動迭代法對試驗數據進行處理和計算，求得軸承的氣剛度系數、阻尼系數。

驅動電主軸進行試驗測試，數據處理系統分別提取30000，40000，50000，60000，70000，80000 r/min時軸承的振動位移、氣膜力和轉速等試驗數據，根據方程（5）編寫軸承剛度系數、阻尼MATLAB求解程序，對試驗數據進行處理計算，求解軸承的剛度系數、阻尼系數，提取所需軸承的瞬時剛度系數、阻尼系數并取平均值，得到軸承剛度系數、阻尼系數與轉速的關系如圖11所示。

4.3 試驗數據分析

從圖11中的試驗數據中可以看出：隨著轉速的增加，軸承的主剛度系數Kxx，Kyy均增大，交叉剛度系數Kxy，Kyx變化較緩；主阻尼系數Cxx，Cyy逐漸減小，交叉阻尼系數Cxy，Cyx則趨近于0。

仿真結果與試驗結果的變化趨勢是相似的，雖然在數值上存在一定差異，但誤差大小在25%以內。這是由于試驗臺具有軸承加工誤差、外界環境干擾等不可控因素；而仿真計算時具有較多理想化因素。從整體上看，仿真結果與試驗結果的一致性較好，驗證了基于ANSYS Workbench軟件仿真計算方法的正確性與可參考性。

5 結 論

（1）根據箔片氣體軸承的結構特點和潤滑機理，基于ANSYS Workbench平臺，通過將流固耦合與箔片氣體軸承的流體潤滑理論相結合，構建了適用于箔片氣體軸承的流固耦合模型，為相似結構特點的模型提供了一種思路。

（2）轉速對軸承動態特性的影響較大，主剛度系數主要起承載作用，隨著轉速的增加，軸承的承載能力得到提升；主阻尼系數主要起抑制渦動的作用，隨著轉速的增加，軸承的穩定性下降，更易發生軸承失穩現象。

（3）波箔片結構參數的改變會引起軸承結構剛度的變化，改變氣膜與箔片間的耦合關系。對于箔片氣體軸承，當波箔片的長度比在1～1.5之間，厚度為0.16 mm時，軸承具有較大的剛度和阻尼，可以提高軸承的承載能力，改善動態特性，提高軸承運行穩定性。

參考文獻

1

王云飛. 氣體潤滑理論與氣體軸承設計［M］. 北京： 機械工業出版社， 1997： 213-215.

WANG Yunfei. Gas Lubrication Theory and Gas Bearing Design［M］. Beijing： China Machine Press， 1997： 213-215.

2

賈晨輝， 張海江， 邱明， 等. 球面螺旋槽動靜壓氣體軸承動態特性試驗［J］. 航空動力學報， 2019， 34（8）： 1847-1856.

JIA Chenhui， ZHANG Haijiang， QIU Ming， et al. Test on dynamic characteristics of spherical spiral groove aerodynamic-aerostatic gas bearings［J］. Journal of Aerospace Power， 2019， 34（8）： 1847-1856.

3

李昊， 耿海鵬， 孫巖樺， 等. 箔片安裝位置對氣體軸承性能影響的仿真分析［J］. 西安交通大學學報， 2019， 53（3）： 81-87.

LI Hao， GENG Haipeng， SUN Yanhua， et al. Simulation analysis for influence of installation position on gas bearing performances［J］. Journal of Xi’an Jiaotong University， 2019， 53（3）： 81-87.

4

梁波， 羅欣洋， 寧家強， 等. 基于三維有限元波箔片模型的氣體箔片軸承承載性能研究［J］. 節能技術， 2021， 39（4）： 332-338.

LIANG Bo， LUO Xinyang， NING Jiaqiang， et al. Study on bearing performance of gas foil bearing based on three-dimensional finite element bump foil model［J］. Energy Conservation Technology， 2021， 39（4）： 332-338.

5

Zhao Xuewei， Xiao Shuhong. A finite element model for static performance analysis of gas foil bearings based on frictional contacts［J］. Tribology Transactions， 2021， 64（2）： 275-286.

6

馮凱， 胡小強， 趙雪源， 等. 三瓣式氣體箔片徑向軸承的靜動態特性［J］. 中國機械工程， 2017， 28（15）： 1826-1835.

FENG Kai， HU Xiaoqiang， ZHAO Xueyuan， et al. Static and dynamic performances of a three-pad gas foil journal bearing［J］. China Mechanical Engineering， 2017， 28（15）： 1826-1835.

7

Xu Haojie， Yang Jiaopeng， Gao Lei， et al. The influences of bump foil structure parameters on the static and dynamic characteristics of bump-type gas foil bearings［J］. Proceedings of the Institution of Mechanical Engineers， Part J： Journal of Engineering Tribology， 2020， 234（10）： 1642-1657.

8

吳垚， 楊利花， 徐騰飛， 等. 氣體動壓徑向軸承超薄氣膜潤滑動特性分析［J］. 振動工程學報， 2019， 32（5）： 908-917.

WU Yao， YANG Lihua， XU Tengfei， et al. Analysis on dynamic characteristics of the ultra-thin gas film lubrication in cylindrical gas journaI bearing［J］. Journal of Vibration Engineering， 2019， 32（5）： 908-917.

9

Jia Chenhui， Cui Zhiwu， Qiu Ming， et al. Flow field calculation and dynamic characteristic analysis of spherical hybrid gas bearings based on passive grid［J］. Journal of the Brazilian Society of Mechanical Sciences and Engineering， 2019， 41（1）： 1-18.

10

廖靜. 基于流固耦合的水潤滑橡膠軸承潤滑特性研究［D］. 重慶： 重慶大學， 2014.

LIAO Jing. The lubrication performance research of water-lubricated bearing based on fluid-structure interaction［D］. Chongqing： Chongqing University， 2014.

11

Wang Bin， Ding Qian， Yang Tianzhi. Soft rotor and gas bearing system： two-way coupled fluid-structure interaction［J］. Journal of Sound and Vibration， 2019， 445： 29-43.

12

焦映厚， 徐磊， 陳照波， 等. 波箔型徑向氣體軸承動態特性分析［C］∥ 第10屆全國轉子動力學學術討論會論文集. 昆明， 2012： 155-160.

JIAO Yinghou， XU Lei， CHEN Zhaobo， et al. Dynamic characteristics analysis of wave foil radial gas bearings［C］∥ Proceedings of the 10th National Symposium on Rotor Dynamics. Kunming， 2012： 155-160.