一種滑動軸承油膜性能計算的動網格方法

吳 超 尹雪梅 李蒙蒙 李奕君 王 文

1.鄭州輕工業大學機電工程學院，鄭州，4500022.鄭州輕工業大學能源與動力工程學院，鄭州，4500023.上海大學機電工程與自動化學院，上海，200072

0 引言

計算油膜軸承性能的關鍵問題就是求解動態流體潤滑方程，得到油膜的壓力分布[1-4]。通常情況下，采用動態雷諾方程計算軸承的性能時，難以精確反映轉速所引起的油流周向慣性效應、動態擠壓效應和靜壓效應之間的線性耦合關系及其對油膜三維壓力場、溫度場和速度場的影響，因此，有必要直接求解Navier-Stokes方程，精確研究軸承參數對油膜性能的影響[5]。文獻[6-8]采用計算流體動力學(computational fluid dynamics，CFD)開展了靜壓軸承承載特性的研究，證實了在表征復雜求解域流體流動形態方面Navier-Stokes方程可彌補雷諾方程的不足。文獻[9-10]利用Fluent軟件，采用靜網格方法，研究了氣穴現象對滑動軸承性能的影響，但靜網格方法不能計算變載荷下軸心的軌跡。文獻[11-12]將軸頸的旋轉動邊界轉換為靜邊界，并利用Smoothing動網格模型，成功計算了油膜軸承的動特性系數。文獻[13-15]自定義動網格更新程序，提出了基于彈性變形的動網格調整法，能用來求解瞬態軸心軌跡，并利用弱耦合算法研究了軸承剛度隨方向的變化以及渦動中心與載荷、不平衡量的關系。文獻[16]針對求解復雜轉子-軸承系統非線性動力學特性的問題，基于Smoothing動網格技術，提出了一種計算流體力學和計算轉子動力學的流固耦合新方法，計算結果表明，該方法能夠得到精確的軸心軌跡，并能準確分析復雜轉子-軸承系統非線性動力學特性。文獻[17]利用非定常動網格技術建立了考慮軸頸渦動頻率與渦動軌跡的滑動軸承動力特性求解模型，研究了不同的軌跡下軸頸渦動頻率和偏心率對滑動軸承動力特性的影響。以上在油膜軸承性能計算中運用的動網格方法是基于Smoothing模型或其變形模型基礎上進行的，能夠方便地求解出油膜軸承的特性參數，分析軸承-轉子系統的動力學性能。

Smoothing動網格模型[18]的特點是不改變網格節點間的拓撲關系，只進行網格形狀的改變，能夠計算變載荷下軸承的軸心軌跡，但Smoothing模型網格位移的更新方法容易造成網格畸變，當偏心率過大時，垂直于軸頸的網格線就會出現嚴重傾斜，甚至出現負網格的情況，導致計算發散。鑒于Fluent軟件在表征流體流動形態方面的優勢及其計算油膜軸承動態性能方面的不足，本文提出了一種基于Fluent方法的滑動軸承油膜性能計算的動網格更新方法。在它的每一時間步內，所有網格節點更新方式是依據軸頸中心坐標確定的，使徑向網格線始終垂直于軸頸表面，保證了移動前后膜厚方向的網格線沿圓周方向上等均分布，網格不發生扭曲變形，避免了網格節點更新產生的累計誤差。通過與典型算例和實驗結果對比分析來驗證該方法的有效性，并考察了進油壓力和載荷對轉子靜平衡位置的影響。

1 動網格更新算法原理

1.1 已有動網格方法

動網格模型可以用來模擬由于邊界運動而引起的流域隨時間的變化。邊界運動可以是主動的規定動作(例如，用戶可以定義剛體中心的速度和角速度)，也可以是被動的待求解的動作(例如，六自由度模型中，已知運動狀態的剛體邊界，受合力后確定下一時刻的速度和加速度)。根據每一個時間步的運動邊界的位置，Fluent軟件自動對動網格更新處理。對于動網格模型，用戶需要確定初始網格和運動區域。3種常用的動網格更新算法[18]分別是Smoothing、Layering和Remeshing。Smoothing方法是不改變網格節點間的拓撲關系，只進行網格形狀的改變，但當網格出現畸變時計算容易發散。Layering方法是隨著動邊界的移動，在邊界處發生網格的增加或合并。Remeshing方法是將控制區內的所有網格重新劃分，一般適用于非結構化網格。理論上，考慮到油膜幾何尺寸的特點，應當選擇Smoothing動網格更新方式，但在軸承性能計算過程中也出現了網格畸變，顯然Smoothing方法也不能很好地解決軸承油膜的瞬態計算問題。

文獻[13]提出的基于Smoothing的網格更新模型示意圖見圖1。網格節點位移按公式Δxi=niΔx/n和 Δyi=niΔy/n進行計算(n為膜厚方向網格的層數，ni為網格節點所在的層數，Δx、Δy為軸頸中心在該時間步的位移)。利用該方法，垂直于軸頸表面的網格線就會出現傾斜，導致計算結果誤差較大，特別是多步計算后的位置累計誤差。

圖1 文獻[13]網格更新方式Fig.1 Mesh updating method in references

1.2 自定義動網格更新算法原理

為解決上述難點，本文提出了一種新的網格更新算法，在網格更新前后，膜厚方向的網格線始終指向軸頸的中心，從而避免了出現網格傾斜。本文網格結構更新算法的示意圖見圖2，網格節點位置計算原理見圖3。

圖2 本文提出的新的網格更新方式Fig.2 New mesh updating method in this paper

圖3 網格節點更新原理圖Fig.3 Schematic diagram of mesh node updating

為使網格不發生扭曲變形，就需要指向軸頸中心的網格線方向始終不改變，即始終垂直于軸頸表面，同時要保證在移動前后膜厚方向的網格線沿圓周方向上等均分布(膜厚網格線的夾角相等)，即更新前后網格線平行，也就是圖3中的向量O1N1與O2N2共線。具體方法如下。

記軸承的中心為O，移動前后的軸頸的中心分別為O1和O2，油膜網格在膜厚方向上分成k層網格，假設N1為t1時刻網格上任意一個待更新節點，節點N1以及與其相對應的更新后的節點N2均在第ki層網格上，根據此刻軸頸所受的油膜力、軸頸坐標、軸頸中心位移速度、軸頸質量可以得到下一時刻t2的軸頸中心坐標，根據時刻t2的軸頸中心坐標和該節點相對軸心的單位向量(與向量O1N1平行或者共線)，作射線分別交軸頸圓和軸承圓于點A2、B2。又知待求點N2所在的層數為ki，即已知LA2N2=(ki/k)LA2B2，LOB2為軸承半徑，LO2A2為軸頸半徑。N2坐標求解算法步驟如下。

(1)求出LOO2；

(2)已知向量O1N1的單位法向量e，由于向量O1N1與向量O2N2共線，故可得到向量O2N2的單位法向量e；

(3)求解向量O2O與向量O2B2的夾角余弦cosα；

(4)利用三角形關系式求解LO2B2：

(5)求出LA2B2=LO2B2-LO2A2；

已知軸頸中心移動后的坐標O2，根據上述算法可求出節點位置N1移動到新位置的坐標N2。

2 軸心靜平衡位置的求解方法

2.1 數學模型

根據牛頓第二定律，運動方程[13]為

(1)

式中，m為轉子質量的一半；S為轉子中心距原點的位移，G=mg；F為轉子受的油膜合力。

已知時間步Δt，上一時間步對應的速度v0，根據牛頓第二定律求解每一個時間點軸頸中心的相對位移和速度計算公式如下：

(2)

v=v0+(G+F)Δt/m

(3)

求解油膜力F的分量Fx和Fy的公式如下：

(4)

(5)

式中，Fx、Fy分別為x軸和y軸方向的油膜力；A為某一網格單元的面積；p為根據N-S方程計算出的軸承油膜各點的壓力。

2.2 求解軸心靜平衡位置流程及宏函數的選用

基于Fluent軟件，利用該方法計算油膜軸承所支撐軸頸靜平衡位置的流程如圖4所示，首先設定軸頸初始速度和初始位移，通過udf宏DEFINE_GRID_MOTION計算出油膜力Fx和Fy，再根據式(2)和式(3)計算出軸心坐標和速度。同時按照新提出的算法原理進行網格更新，在該時間步內迭代完后進入下一個時間步。下一個時間步用到上一時間步的軸心坐標和速度(讀取儲存于文件2和文件3最后一次更新的數值)。一直循環下去，直到軸頸中心趨近于一定點，即可停止計算，該定點即穩態下徑向油膜軸承所支撐軸頸中心的靜平衡位置。

圖4 軸心靜平衡位置計算流程圖Fig.4 Computational flow chart about static balance position of the rotor

一般來說，動網格節點的位移由坐標確定，但滑動軸承油膜是一種長徑比很大的網格，用坐標判斷節點位置容易出現錯誤，可采用節點全局編號判定節點的位置，強制精確控制每一個節點的位移變化。利用DEFINE_GRID_MOTION宏對網格節點位置進行強制性精確定義，可有效降低網格畸變的可能性，由于軸頸表面各點的速度不一樣，可通過編寫DEFINE_PROFILE宏程序控制軸頸表面各點的速度。按照以上步驟，可以利用動網格方法實現對滑動軸承油膜狀態的仿真。

3 算法的正確性和有效性驗證

本文提出的網格更新方法與文獻[13]所使用的網格更新方法的主要區別就是更新后油膜厚度方向的網格線是否垂直于軸頸表面。要驗證所提動網格算法的正確性和有效性，需要進行兩方面的驗證：用于靜態油膜力計算的精度驗證；用于軸頸中心靜平衡位置求解的累積誤差驗證。

3.1 算例1：用于靜態油膜力計算的精度驗證(單次油膜力計算)

利用所提出的網格結構和文獻[13-15]的網格結構，分別計算軸頸轉動角速度為500 rad/s、偏心率為0.5和0.9的圓柱軸承的油膜力，來對比這兩種網格在靜態油膜力計算(用于靜網格計算時，單次網格計算)的精度。算例選取的軸承幾何參數和潤滑油物性參數見表1。軸承供油方式采用兩側雙向進油并且設置軸向油槽。

由于油膜厚度尺寸很小，故油膜厚度方向的網格劃分對計算結果影響最大，需要先進行網格獨立性驗證。網格獨立性計算時選擇多項流混合模型，氣穴模型選用Singhal全空化模型，連續性方程、動量方程和氣穴方程均選用一階迎風格式進行求解，速度耦合格式選擇SIMPILE格式。軸

表1 軸承和潤滑油物性參數

承油膜網格初始劃分的網格數(膜厚方向×周向×軸向)為5×600×80，通過改變膜厚方向的網格層數，并以凈流量小于進口和出口流量兩者之間的最小值的1%作為判斷收斂的標準。膜厚方向上設置不同的網格層數，以承載力作為觀測量，網格獨立性驗證結果見表2。結果表明，本文所用的網格密度計算結果的最大偏差均小于4%。

表2 網格獨立性驗證結果

綜合考慮計算效率和計算精度，把油膜膜厚方向網格劃分為6層(6×600×80)，見圖5。算例中壓力-速度耦合格式分別選SIMPILE、SIMPILEC、PISO方式，其他設置同上。采用瞬態計算方法，利用文獻中的網格方法和本文的網格方法計算得到的承載力結果和圓周方向的壓力分布結果(取軸承寬度中間L/2處，偏心率為0.9)分別見表3和圖6。

圖5 油膜網格Fig.5 Oil film mesh

表3 兩種靜網格結構計算的承載力對比

(b) 本文方法的計算結果

由表3和圖6可以看出，基于兩種網格方法得到的軸承承載力和瓦塊圓周方向的壓力計算結果相同，均與文獻[19]結果相接近，表明兩種網格(油膜厚度方向的網格線垂直于或者不垂直于軸頸表面)用于軸承油膜穩態油膜力計算時，均能滿足精度要求(<9%)。由于在計算過程中網格不更新(靜網格法)，無論油膜厚度方向的網格線是否傾斜，累計計算誤差都不存在，故兩種網格的計算結果基本無差別。這就證明了Smoothing網格方法和本文提出的網格方法來計算軸承穩態油膜力(靜網格法)的精度是滿足要求的。

3.2 算例2：用于靜平衡位置求解的累計誤差驗證(多次油膜力計算)

根據文獻[13]的算例，選取軸頸角速度為500 rad/s，轉子質量為17.527 kg(171.76 N)，軸承、潤滑油的參數和油膜網格劃分同算例1，對軸頸中心軌跡的收斂過程及靜平衡位置進行了仿真，結果見圖7，計算得到其偏心率為0.221。根據文獻[19](未給出進油壓力，其他條件完全一樣)，在相同軸承參數下，據此載荷計算出軸承偏心率為0.267，與本文計算結果相比，相對誤差是17.2%。文獻[13]和本文在計算時考慮了氣穴的影響，而文獻[19]沒有考慮氣穴的影響。

(a) 二維圖

(b) 三維圖

文獻[13]的靜平衡位置所對應的偏心率為0.046。將文獻[13]結果和本文的計算結果分別與文獻[18]相比較，文獻[13]的計算結果差別更大，偏心率明顯不在同一數量級，表明文獻[13]的網格經多次更新后，出現了較大的累計誤差。這與文獻[13]在網格更新時沿膜厚方向的網格線與軸頸表面不垂直有關，也與其網格節點坐標更新的計算方式與軸頸表面的位置有關，從而導致多次迭代后求解軸頸中心靜平衡位置時出現了較大的誤差。

為驗證所編制程序的穩定性，從3個不同起始位置開始分別計算相同條件下的軸心靜平衡位置，軸心軌跡收斂曲線如圖8所示。可以看出，3條曲線最終收斂于一點，表明計算的起始位置對軸頸的靜平衡位置并無影響，驗證了所編制的網格更新程序的穩定性。

圖8 不同起點位置的軸心收斂軌跡Fig.8 Axial convergence trajectory of different starting positions

3.3 算例3：進油壓力對軸心靜平衡位置的影響

為了分析進油壓力對軸承所支撐轉子軸心靜平衡位置的影響，在給定氣化壓力29 185 Pa、轉子質量17.527 kg、角速度500 rad/s和 000 rad/s的條件下，計算了供油壓力在0.1～0.5 MPa之間時轉子軸心的靜平衡位置。軸承的偏心率和偏位角隨進油壓力的變化關系如圖9所示。算例的其他計算參數見表1。

(a) 偏心率

(b) 偏位角

由圖9可以看出，隨著進油壓力的增大，偏心率減小(<7%)，偏位角增大(<9%)。這是因為在載荷和轉速不變的前提下，隨著進油壓力的增大，軸承油膜的最大壓力增大，承載能力增強，偏心率減小。同時由于進油口在水平方向上，故隨著進油壓力增大，偏位角增大。計算結果表明進油壓力對軸承的承載性能有影響，但影響率小于10%。結合算例2的結果，進一步證明了所提出的動網格更新方法用于軸心軌跡計算時累計誤差相對較小。

通過以上算例驗證了所提出的動網格方法在油膜性能計算的有效性、可行性和穩定性。

4 結論

(1)本文提出了一種用于滑動軸承油膜性能計算的動網格更新方法，驗證了其正確性。本文方法能夠保證在網格移動前后膜厚方向的網格線沿圓周方向上等均分布，并使指向軸頸中心的網格線方向始終垂直于軸頸表面，網格不發生傾斜。

(2) 本文方法應用于求解滑動軸承所支撐轉子的靜平衡位置時，能夠減少網格計算的累計誤差，提高計算的精度；同時發現利用Fluent軟件動網格方法計算油膜軸承流場時，膜厚方向網格傾斜會對計算結果產生較大的累計誤差。

(3)在相同的條件下，隨著進油壓力增大，軸承偏心率減小，偏位角增大，但進油壓力對軸承的承載性能的影響率一般小于10%。