基于3D瞬態流場計算的滑動軸承非線性油膜力分析

2018-11-01 01:21:14許偉偉

振動與沖擊 2018年20期

李強，張碩，許偉偉

(1. 中國石油大學(華東) 化學工程學院，山東青島 266580; 2. 中國石油大學(華東) 儲運與建筑工程學院，山東青島 266580)

液體動壓徑向滑動軸承以其承載能力大、功耗小、運轉精度高等優點廣泛應用于汽輪機組、工業壓縮機、核電機組及各類大型機床，是機械工業中使用廣泛、要求嚴格的關鍵基礎部件之一。

滑動軸承油膜力既是轉子-軸承系統阻尼的主要來源，也可能是導致機組穩定性下降的重要原因，是否能準確獲得非線性油膜力直接決定了轉子-軸承系統非線性動力學分析的計算精度和效率。目前普遍采用的非線性油膜力計算方法主要有4類，一是解析或半解析方法，如將滑動軸承簡化為無限長、無限短軸承油膜力模型[1-2]；二是采用有限元或有限差分法直接求解Reynolds方程[3-4]；三是基于Poincare變換的數據庫方法；四是基于CFD的流場數值計算方法，其中前3類是基于Reynolds方程的方法。近年來，隨著大型商用計算流體力學軟件的出現，人們開始通過直接求解Navier-Stokes方程來研究滑動軸承的動力學行為。Chen等[5]利用CFD軟件對滑動軸承、階梯軸承、阻尼器等進行了穩態計算，討論了在傳統的求解Reynolds方程的方法中被忽略的慣性項、黏性項的變化對上述結構性能的影響；Lin等[6]利用CFD技術對轉子-軸承系統的潤滑性能進行了分析，詳細討論了軸承的熱效應和空穴效應；Biswas等[7]利用FLUENT軟件對氣體潤滑三油楔軸承進行了瞬態分析，并比較了不同時刻的壓力和溫度值；Panday等[8]利用FLUENT軟件研究了滑動軸承的瞬態動力學行為，并計算了不同長徑比下滑動軸承的承載力、摩擦力等；Li等[9]提出了基于FLUENT的滑動軸承流體模型和基于FORTRAN的多圓盤撓性轉子動力學模型的耦合算法，并分析了軸頸傾斜對轉子-軸承性能的影響。采用CFD建模計算，不僅可以考慮軸承內部的徑向流場變化、流體慣性項的影響，對于事先未知的自由邊界等問題更容易處理，理論上適用于任意軸承結構型式，因此基于CFD技術的3D流場計算方法越來越廣泛地應用于滑動軸承性能的研究。

目前，滑動軸承3D流場計算最大的問題是，滑動軸承油膜流場在三個空間尺度上的尺寸大小非常不均勻，如果采用CFD軟件中已有的動網格技術在軸頸移動幾步后就會出現比較大的網格畸變，使瞬態流場無法繼續計算[10-11]。為了準確計算滑動軸承非線性油膜力，通過研究，本文提出一種理論上適用于任何固定瓦滑動軸承瞬態流場計算的變流域動網格技術，保證了瞬態流場計算的順利進行。在此基礎上，本文首先計算并對比了小擾動下的線性和非線性油膜力，發現數值計算結果與理論分析一致，驗證了數值計算方法的正確性；其次，通過對大擾動工況下滑動軸承油膜力的瞬態計算發現，大擾動工況下的轉子-軸承系統非線性動力學分析有必要計入滑動軸承油膜非線性的影響；最后，在考慮了油膜力的非線性因素下，闡明了油膜力的切向和徑向分量與滑動軸承穩定性的關系，以及瞬態潤滑流場下油膜力分量的變化過程。

1 控制方程

FLUENT中所求解的連續性方程和動量方程對于可壓縮和不可壓縮流體的流動均適用。其中的動網格模型可以用來模擬由于邊界運動引起的流域形狀隨時間變化的瞬態流動。

連續性方程為

(1)

動量方程為

(2)

(3)

式中：μ是流體黏度;I是單位應力張量; 式(3)右邊第二項表示體積膨脹效應。

2 計算模型與方法

2.1 計算模型

本文主要對比討論不同幅值的擾動載荷下油膜力的變化情況，而與具體的轉子系統無關，因此計算模型簡化為滑動軸承支撐下的單跨轉子進行分析，計算模型如圖1所示，應當指出，為了方便說明，軸承間隙被放大，其中轉子的質量為M=17.527 kg。

圖1 轉子計算模型Fig.1 Schematic diagram of rotor

滑動軸承軸頸中心的運動方程為

(4)

式中：M為軸承承受轉子的質量;e為偏心距;ω為轉子角速度;Fx、Fy為油膜力在x,y方向的分力，是通過對油膜壓力積分得到的，計算公式為

(5)

式中：R為轉子半徑;p表示壓力;θ為壓力與豎直方向夾角。

本文在進行非線性油膜力計算時，采用的滑動軸承模型如圖2所示。

圖2 滑動軸承計算模型Fig.2 Schematic diagram of journal bearing

為了與式(5)計算得到的非線性油膜力對比，下面給出線性油膜力計算公式

(6)

式中：Fx0、Fy0分別為靜平衡位置時，油膜力在x,y方向的分量; [K]、[C]分別為剛度系數矩陣和阻尼系數矩陣。

圖3 軸頸中心渦動時油膜力分量關系圖Fig.3 Oil film force components at oil whirling

滑動軸承的動特性是與軸承內的流體動態力密切相關的，如圖3所示。將油膜動壓產生的力沿軸心軌跡分解，分為沿著軸頸渦動方向的切向油膜力Fτ和垂直軸頸渦動方向的徑向油膜力Fr，計算公式為

(7)

式中：Ob為滑動軸承的中心；O′為軸頸的渦動中心； Δx、Δy分別為一個時間步內軸頸中心沿x和y方向的位移值。文中令切向油膜力的方向與軸頸渦動方向相同時為正值，徑向油膜力以指向渦動軌跡曲率中心的方向為負值方向。

2.2 數值計算方法與邊界條件

滑動軸承潤滑流場為小間隙結構，徑向間隙很小，尤其是在軸頸大偏心下，選用高階精度的離散格式容易引起解的振蕩而達不到收斂精度，而且應用本文的動網格方法可生成與流場適應的網格模型從而減輕低階離散格式帶來的假擴散現象。因此本文計算中，在保證計算精度的同時，為了提高計算速度和加快收斂，數值計算方法采用 “基于壓力求解”，壓力速度耦合方法，在穩態計算中為提高計算效率采用SIMPLEC算法，在非穩態計算中為了提高計算穩定性采用PISO算法，連續方程、動量方程采用一階迎風格式，壓力差值格式采用Linear差值，控制方程求解過程中考慮重力的影響，重力加速度設為豎直向下9.8 m/s2，忽略溫度的影響[12-13]。

通過軸承潤滑流場的尺寸和操作條件，可判定軸承流場處于層流狀態；操作壓力設為0 Pa，計算域的進口在兩邊進油口位置，邊界類型為 “壓力入口”，計算域的出口在軸承軸向的兩端，邊界類型為“壓力出口”，出口位置的壓力為大氣壓；軸頸面設為旋轉面，穩態計算時直接指定旋轉速度，動態計算時通過UDF指定旋轉速度；固壁邊界條件設置為無滑移條件，近壁面應用標準壁面函數。

2.3 轉子-軸承系統耦合計算方法

轉子-軸承系統的流固耦合計算采用弱耦合方法，計算流程圖如圖4所示。首先，基于FLUENT環境進行滑動軸承的三維瞬態流場計算，將計算得到的油膜壓力通過式(5)獲得滑動軸承瞬態油膜力并寫入數據文件。其次，將得到的瞬態油膜力和轉子參數代入式(4)可以得到當前時間歩內軸頸中心在x，y方向的加速度，進而求得速度和位移，并寫入數據文件。最后，FLUENT通過自定義程序(UDF)調用上述軸頸位置參數完成軸承流域動網格更新和軸頸表面旋轉速度的設置，然后進行下一個時間歩的計算，直至計算結束。

圖4 瞬態轉子-軸承系統計算流程圖Fig.4 Transient calculation process of rotor-journal bearing system

3 動網格方法

在現有動網格技術的基礎上，針對滑動軸承油膜間隙小的結構特點開發了適用于瞬態潤滑流場計算的動網格方法，該方法中，網格節點的數量不變，節點間拓撲關系保持不變。

圖5為軸頸擾動時，潤滑流場中網格位置變化示意圖。首先，選用六面體單元對滑動軸承潤滑流場劃分網格，保證潤滑流場中每個節點的位置坐標都可以通過計算得到。當軸頸渦動時，流場最內圈上的所有節點都隨著軸頸一起移動，最外圈及油槽上的節點保持不動，中間的節點按一定算法移動，如下式所示。

(8)

圖5 軸承間隙周向網格移動示意圖Fig.5 Mesh updating at radial clearance

該算法編程后可以通過UDF接口載入FLUENT中，并將節點移至新的位置，每完成一個時間步長的計算，就更新一次流域網格。通過大量的數值計算發現，通過該動網格方法移動后的網格即使在很大的軸頸偏心下依然能保證良好的光滑性和規整性，而且不容易出現負體積。

不同動網格方法更新結果如表1所示，由更新后的網格可以看出，FLUENT提供的三種動網格模型對滑動軸承小間隙模型的適用性較差，特別是大擾動情況下，難以保證計算穩定性、計算精度和效率。

表1 不同動網格更新方法比較Tab.1 Comparison of different dynamic mesh methods

4 滑動軸承的油膜力分析

當軸頸在靜平衡位置受到擾動時，軸頸受到的油膜力和擾動之間的關系一般是非線性的。而當軸頸中心在靜平衡位置附近作小振幅運動時，可以近似認為這種關系是線性的，本節對比分析了不同擾動幅值下的瞬態油膜力。

4.1 小擾動下油膜力分析

當轉速取小于一階臨界轉速的ω=800 rad/s，不平衡量取e=1 μm和e=2 μm時，計算得到的軸心軌跡如圖6所示。從圖中可知，當不平衡載荷較小時，軸心在靜平衡位置附近做軌跡為橢圓的小幅度渦動，渦動的頻率為轉頻，這與線性分析結果基本一致。

根據傳統的潤滑理論，小擾動情況下，隨著不平衡量的增大，軸心軌跡和油膜力是以靜平衡位置為中心呈比例增大。對應于圖6軸心軌跡的線性油膜力和非線性油膜力的對比關系如圖7所示。可以看出，當不平衡量e=1 μm時，線性油膜力與非線性油膜力一致性比較好，這時可以將系統進行線性化處理；當不平衡量增加到e=2 μm時，線性和非線性油膜力隨時間的變化規律基本一致，但由于渦動中心地上移，使非線性油膜力和線性油膜力在位置上出現了偏移。

圖6 軸心軌跡隨不平衡量的變化Fig.6 Journal orbits for different unbalance

圖7 小擾動下線性油膜力與非線性油膜力的對比Fig.7 Comparison of linear oil film force and nonlinear oil film force under small disturbance

4.2 大擾動下油膜力分析

圖8和表2分別給出了不同擾動下滑動軸承的軸心軌跡和渦動中心的變化過程。從圖中可知，當不平衡量較小時，軸心在靜平衡位置附近作軌跡為橢圓的小幅度渦動，隨著不平衡量的增加，軸頸振幅增大，軸心軌跡不再收斂于橢圓形狀，且并沒有一個固定的平衡位置，這時當不平衡量再變化時，軌跡的中心在變化，偏離了原來的靜平衡位置，從軸心軌跡來看，即使轉子的轉速和軸承的靜載荷都不變化，只要動載荷發生變化，軸頸的渦動中心也會隨之改變，動載荷越大，渦動中心越靠近坐標原點位置。

圖8 油膜間隙對橢圓軸承軸心軌跡的影響Fig.8 Influence of oil film clearance on the journal orbits of elliptical bearing

傳統的潤滑理論認為，對于確定的轉子-軸承系統，軸頸中心的靜平衡位置決定于軸頸轉速，即軸承的動力特性系數僅僅是軸頸轉速的函數，因此關于滑動軸承動力特性系數的表格都是按照Sommerfeld數排列的，只考慮軸承參數、靜載荷和軸頸轉速的影響，而并沒有考慮動載荷的影響。由上面對不同擾動載荷下軸心軌跡的計算可知，該結論在大擾動下并不成立。

表2 滑動軸承的渦動中心隨不平衡量的變化Tab.2 The whirling centers of the journal bearing with different unbalance

圖9給出了大擾動下根據瞬態流場計算得出的非線性油膜力和線性油膜力的對比關系，從圖中可以看出，根據靜平衡位置處的動特性系數得出的線性油膜力與非線性油膜力差別明顯，這時不只是兩種油膜力曲線的位置不同，兩種油膜力曲線的形狀也產生了明顯變化，此時不能再用線性模型來求解瞬態油膜力。因此在進行轉子-軸承系統的非線性動力學分析時有必要計入滑動軸承油膜非線性的影響。

圖9 大擾動下線性油膜力與非線性油膜力的對比Fig.9 Comparison of linear oil film force and nonlinear oil film force under large disturbance

4.3 瞬態油膜力的徑向和切向分量

圖10為轉速ω=800 rad/s，不平衡量e=2 μm時軸心軌跡所對應的瞬態切向和徑向油膜力，從圖中可以看出，與油膜力分量Fx和Fy相比，瞬態油膜力的切向分量和徑向分量更能反映軸承穩定性的好壞，當徑向油膜力始終為負值時表明油膜的支撐剛度比較大，為軸頸提供足夠的支撐而防止軸頸渦動軌跡變大；當切向油膜力始終為負值時表明系統阻尼足夠大，使軸頸穩定渦動而不失穩。因此可以很容易從圖10中切向油膜力和徑向油膜力的大小判斷出，當不平衡量e=2 μm時軸頸是以工頻穩定渦動的。

圖10 滑動軸承瞬態油膜力分量Fig.10 Transient oil film force components of journal bearing

圖11、12分別給出了不平衡量e=18 μm時不同的軸心渦動軌跡點對應的徑向、切向油膜力和瞬態油膜壓力場。當軸頸中心渦動到A點和B點附近時，受軸承右邊油槽和進油口的影響，收斂楔的長度被大大縮短，降低了油膜的支撐剛度，表現為徑向油膜力很小；當渦動中心繼續向C、D點移動，雖然此時軸心在第二象限，但這時圓柱軸承的上軸瓦并沒有形成明顯的油膜收斂楔而產生大的油膜壓力，反而隨著軸頸中心的下移，逐漸在下軸瓦形成較大的油膜壓力，這時由于整體油膜壓力分布受油槽影響比較少，因此油膜支撐剛度比較理想，但渦動中心的位置比較靠上，形成的壓力峰值比較小，導致切向油膜力比較小；隨著軸頸渦動中心繼續向E點移動，軸承左邊油槽和進油口對油膜支撐剛度的影響慢慢變大，表現為徑向油膜力有變大的趨勢；但當軸心下移到F點附近時，由于這時軸心在x，y兩個方向的偏心都比較大，形成了很大的油膜壓力，而且F點距左油槽的距離明顯比B點距右油槽的距離遠，從壓力云圖上看，壓力分布并沒有明顯的被油槽分割開，因此大的油膜壓力使得這時切向油膜力和徑向油膜力都比較大。

圖11 e=18 μm時滑動軸承的軸心軌跡和瞬態油膜力分量Fig.11 Journal orbit and transient oil film force of journal bearing at e=18 μm

圖12 e=18 μm時不同時刻軸頸油膜壓力云圖Fig.12 The pressure of oil film at different time at e=18 μm

為了研究半頻渦動對切向油膜力的影響，令軸頸以相同的轉速和不同的渦動頻率作相同的渦動軌跡，得到的切向油膜力如圖13所示。從圖中可以看出，當渦動頻率為轉速的一半時，即軸頸進行半頻渦動時，切向油膜力分量會變小，導致系統阻尼減小，容易引起系統失穩。因此當系統出現低頻渦動時，系統的阻尼會大大降低而影響系統的穩定性。