徑向滑動軸承潤滑油膜流動-傳熱過程仿真

孫 正，黃鈺期，俞小莉

(浙江大學 能源工程學院，杭州 310027)

0 引 言

徑向滑動軸承廣泛應用于各類機械傳動機構。隨著動力裝置功率密度的不斷提升，軸承的轉速越來越高、承受的載荷越來越大、同時由于結構緊湊性限制等因素，軸承的冷卻條件也相對惡劣，熱效應對軸承性能的影響變得越來越顯著[1-3]。失效現象時有發生。

通過試驗與數值仿真相結合的手段，研究摩擦副在各種工況下的潤滑性能，為軸承設計提供依據，是解決和避免上述問題的基礎。在傳統的軸承潤滑數值計算中，多采用基于Reynolds方程的方法[4-9]。Reynolds方程的推導需對質量、動量守恒方程引入簡化假設，因此丟棄了大量油膜厚度方向的計算信息，無法精確分析油膜內部的流動傳熱現象，無法處理真實的空化現象。在某些情況下，這些簡化可能導致計算結果與實際情況之間存在較大誤差。

近年來，研究人員轉而嘗試采用計算流體力學(Computational fluid dynamics，CFD)方法求解完整的質量、動量守恒方程，希望能夠克服傳統方法的不足[10-16]。浙江大學Li等[12]計算了軸頸傾斜情況下，潤滑油膜的瞬態流動以及軸承系統整體的動力學特性問題。東南大學張楚等[17]建立了滑動軸承流場數值計算方法，結果表明，CFD的有效載荷計算結果與試驗更加吻合。北京理工大學王康等[18]嘗試采用CFD方法對內燃機連桿大頭軸承潤滑性能進行評估。北京科技大學Li等[19]的研究也表明，Reynolds方法無法準確求解摩擦副表面粗糙度較大的情況。Gertzos等[20]計算了非牛頓流體滑動軸承的潤滑性能。重慶大學孟凡明等[21]研究對比了不同CFD求解器對最大油膜壓力、承載力等計算結果的影響。通過分析上述文獻可以發現，相比于Reynolds方法，CFD方法可考慮的影響因素更為全面、適用范圍更為廣泛，可為動壓潤滑數值仿真計算提供更為豐富和更為精確的結果，將成為未來動壓潤滑數值仿真的趨勢。

雖然CFD方法已應用于流體動力潤滑過程的研究，并且為模擬動壓潤滑中的空化現象提供了必要的基礎，但是在目前多數研究中，空化現象的處理仍采用簡化假設的方法[1，20，22-27]，尚未引入全空化模型[28]。同時，現有研究對熱效應的影響也關注較少，鮮有探討不同空化模型對溫度分布計算結果的影響。

本文在已有的研究基礎上，引入空化模型，并考慮熱效應的影響，進一步完善已有的CFD潤滑計算方法。并將CFD方法的計算結果與文獻[29,30]中的試驗和Reynolds方程仿真結果進行對比，分析不同仿真方法以及空化模型對徑向滑動軸承流動-傳熱計算結果的影響。

1 研究對象

本文以應用廣泛的徑向滑動軸承為研究對象。軸承模型示意圖如圖1所示。參考文獻[29，30]對徑向滑動軸承進行系統的試驗，獲得了壓力分布以及溫度分布試驗數據。因此，本文將文獻[29，30]中的軸承模型選為本文的仿真對象。軸承系統具體參數如表1所示。載荷施加在軸瓦上，方向沿圖1中x軸方向向上。

圖1 軸承模型示意圖Fig.1 Schematic diagram of journal bearing model

表1 軸承參數Table 1 Journal bearing parameters

2 仿真模型建立

2.1 控制方程

通用形式的連續性方程為：

(1)

動量守恒方程為：

(2)

能量方程為：

(3)

式中：k為導熱系數；Sh為體積熱源項；E的表達式如下：

(4)

2.2 空化模型

(5)

式中：RB為氣泡半徑；σ為液體表面張力；ρl為液相密度；PB為氣泡表面壓力；P為遠場壓力。

采用Singhal模型描述氣穴生成和潰滅的過程，考慮液相和汽相的兩相連續性方程為：

(6)

式中：fν為汽相質量分數；Γ為擴散系數；Re、Rc分別為汽化和液化相變率。

根據Singhal模型：

當P≤Pv時，有：

Re=

(7)

當P>Pv時，有：

(8)

式中：Fvap取0.02；Fcond取0.01；κ為湍動能；fg為不可溶氣體質量分數；ρv為汽相密度；Psat為飽和蒸汽壓力；Pv為修正后的飽和蒸汽壓力，其計算公式為：

領導要重視和支持健美操在高校的發展，制定新的政策，保障健美操發展所需要的資金；加大力度改善場地、器材設備的狀況，從外部硬件上促進學校健美操教學的發展；重視體育教師培訓工作，要經常過問并督促檢查本校的體育教師培訓工作的開展情況。

(9)

2.3 邊界條件

仿真計算中的邊界條件設置如下：將圖1中的供油槽設置為壓力入口，潤滑油入口壓力和溫度分別如表1所示；將軸承間隙兩側端面設為壓力出口，表面壓力為0；將軸表面設置為旋轉面(Moving wall)，由于軸高速旋轉，假設軸表面溫度分度均勻，溫度值與文獻[29,30]中取值相同，如表1所示；軸瓦固定不動，外表面熱邊界條件設置為第3類邊界條件，取值見表1。

2.4 求解器設置

采用ANSYS Fluent 13.0軟件作為計算平臺，選用3D-Double Precision-Steady-Pressure Based求解器；選用Mixture多相流模型；由于本文算例的Re數為47.12～94.24，因此黏性流模型設置為層流；考慮流體黏性生熱；潤滑油黏度采用變物性設置，黏度-溫度關系采用與文獻[29，30]中相同的表達式；迭代求解方法采用SIMPLE算法，各物理量的空間離散均采用Second Order Upwind設置。求解過程中，計算開始時需適當調小能量和汽相質量分數的松弛因子，防止計算初期的發散，待計算殘差穩定后將松弛因子調回默認值。

2.5 網格無關性分析

由于潤滑油膜各個方向上都包含流動-傳熱信息。而潤滑油膜在油膜厚度方向上尺度很小，一般為10～100 μm，但油膜周向和軸向的尺度都為數十毫米。因此，仿真計算結果對網格質量非常敏感。在不遺漏各個方向上計算結果的基礎上，同時又能提高計算收斂速度、節省計算資源，合理的網格劃分變得非常重要。本文采用如表2所示的4種網格劃分方案，選取油膜內最大壓力和最高溫度值進行網格無關性分析。如圖2所示，油膜內最大壓力和最高溫度計算結果在網格劃分方案3以后已趨于穩定。因此，本文采用表2中網格劃分方案3進行后續研究，即油膜厚度方向布置18層網格；軸向400層網格；周向布置360層網格。

表2 網格無關性分析Table 2 Grid independence analysis

圖2 網格無關性分析Fig.2 Grid independence analysis

3 計算結果及討論

3.1 壓力分布計算結果

圖3給出了采用不同數值仿真方法，軸瓦中間截面周向壓力分布(逆時針方向)計算結果與試驗值的對比情況。

圖3 軸瓦中間截面周向壓力分布Fig.3 Circumferential pressure distribution on middle plane of bearing bush

圖3中，Reynolds方法計算結果和試驗數據均來自文獻[29，30]。從圖3可以看出：傳統Reynolds方法的計算結果，壓力最高點與試驗值的誤差分別為16.1%、3.6%和8.7%；CFD方法由于采用了黏-溫變物性設置，壓力最高點計算結果與試驗值誤差分別為4.0%、4.1%和1.2%。同時，CFD方法的壓力分布計算結果在總體趨勢上也更接近試驗值。

3.2 溫度分布計算結果

圖4為采用不同數值仿真方法，軸瓦中間截面周向溫度分布(逆時針方向)計算結果與試驗值的對比情況。

圖4 軸瓦中間截面周向溫度分布Fig.4 Circumferential temperature distribution on middle plane of bearing bush

圖4中Reynolds方法計算結果和試驗數據也均來自文獻[29,30]。從圖4可以看出，CFD方法獲得的結果無論是在最高溫度點還是整體趨勢上，都與試驗值更相符，而Reynolds方法獲得的最高溫度點計算結果普遍低于試驗值。這為上文壓力分布的計算結果提供了解釋。從表1中的黏-溫特性值可以發現，潤滑油黏度對溫度變化非常敏感，而潤滑油黏度的空間分布會極大地影響壓力分布的計算結果。由于求解時采用CFD方法，因此獲得了更為精確的流動、傳熱、黏性產熱計算結果，進而獲得了更精確的潤滑油三維溫度分布。與此同時，通過Fluent求解器中的變物性設置，將溫度分布與黏度分布實時耦合，進而獲得了更精確的壓力分布結果。上述分析也在表3中得以體現，CFD方法的壓力、溫度計算結果均與試驗值符合良好。

從圖4中還可以看出，采用CFD方法計算溫度分布，在壓力建立段(即圖1中幾何結構收斂段)效果較好。而在幾何結構發散段，CFD方法能夠模擬出溫度分布的驟降，但過于劇烈，后段出現較明顯的偏差；傳統Reynolds方法的溫度分布結果則顯得平緩，兩者與試驗值都有一定差異。

分析這一現象的原因可能是：后段潤滑油膜內的氣穴破裂，壓力驟降，因此下游低溫潤滑油在背壓作用下，回流進入氣穴潰滅形成的真空區域，與氣穴破裂液化后形成的少量液態潤滑油混合?；旌想A段中，由于倒流的潤滑油在質量上占主要權重且溫度偏低，因此出現了溫度的驟降。在CFD仿真中，Singhal空化模型判斷油膜內壓力一旦低于潤滑油的氣化壓力，則氣穴全部立即開始潰滅，反映在圖上就形成了一個明顯的溫度拐點。在實際過程中，氣穴的潰滅是相對漸進的過程，因此溫度分布試驗值曲線呈緩慢光滑下降。傳統方法完全不能模擬氣穴生成和潰滅過程，只是采用Reynolds邊界條件或者Half-Sommerfeld邊界條件簡單地將負壓置零，發散段傳熱采用簡化處理，因此溫度分布過于平緩，與溫度試驗結果曲線的趨勢存在較大差異。

總體上看，CFD方法獲得的溫度分布在趨勢上與試驗值更為接近，在氣穴處理以及發散段溫度計算中也更接近物理實際，與傳統方法相比更具優勢。但現有的空化模型在判斷、模擬氣穴潰滅過程這一點上仍有進一步改進的空間。

表3 仿真結果與試驗值標準差對比Table 3 Standard deviation of simulation results

3.3 兩相流計算結果

圖5為基于CFD方法引入全空化模型后，潤滑油膜內汽相百分比的計算結果云圖，以及產生氣穴的區域與油膜壓力的對應關系。圖5和圖6中軸頸旋轉方向均為逆時針方向。從圖中可看出：在潤滑油膜的發散段，出現了明顯的空化區域。傳統Reynolds方法，采用Reynolds邊界條件或者Half-Sommerfeld邊界條件簡單地將負壓置零，并不直接模擬氣穴的產生和潰滅過程，因此無法得到類似圖5的計算結果。同樣地，采用CFD方法但不引入全空化模型也無法獲得豐富的空化區域三維計算云圖。引入全空化模型的CFD方法可為進行徑向滑動軸承各項設計參數、結構參數對氣穴形成區域的影響規律等后續研究提供數值仿真手段，并可以獲得豐富的兩相流計算結果，相比于Reynolds方法和以往CFD方法的局限，體現出明顯的優勢。

圖5 計算結果云圖Fig.5 Contours of simulation results

圖6 空化模型對溫度分布的影響Fig.6 Influence of cavitation models on temperature distribution

3.4 空化模型對溫度分布計算結果影響

在以往文獻中[20,22]，雖然已采用了CFD方法求解徑向滑動軸承流動問題，但仍采用簡化方法處理空化現象，具體做法是利用Fluent UDF中的DEFINE_ON_DEMAND命令對可能產生氣穴的區域施加Half-Sommerfeld邊界條件。但在文獻中只展示了采用這種方法后，壓力分布計算結果與試驗值的對比，并未討論其對溫度分布的影響。

圖6為不同空化模型對溫度分布的影響，其中，軸承轉速為2000 r/min，載荷為4 kN。從圖6可以看出，不同空化現象處理方法對徑向滑動軸承幾何發散區域的溫度分布計算結果有明顯影響：①采用Half-Sommerfeld邊界，只對壓力計算結果做置零處理，因此發散區域是理應產生氣穴的區域，在數值計算處理中，完全由出口和間隙側面壓力邊界中的回流填充，而回流的潤滑油溫度取決于人為設定值，圖中上方的藍色區域即反映了這種情況。②采用Singhal模型處理空化，由于考慮了實際生成氣穴的物理過程，抑制了上述不合理的潤滑油回流問題。

因此，若想要獲取精確的徑向滑動軸承溫度分布計算結果，需要考慮氣穴的影響，引入完整的空化模型，僅采用簡化的空化現象處理方法所獲得的結果與實際情況存在較大差異。

4 結 論

(1)相比于傳統方法，采用CFD方法獲得的油膜最高壓力點計算結果與試驗值更接近；同時整體壓力分布也與試驗值更相符。

(2)在CFD方法中引入全空化模型，可以獲得更精確的流動、傳熱、黏性產熱計算結果，從而獲得更精確的潤滑油三維溫度分布。結合黏-溫分布的實時耦合手段，可使油膜壓力分布計算結果更精確。而僅在CFD方法中應用簡化空化現象處理方法時，溫度分布計算結果與實際存在較大差異。

(3)CFD方法在求解氣穴生成-潰滅、流動傳熱、黏性產熱等流動-傳熱問題上，采用的數值處理方法更符合物理實際。尤其是在對空化現象的模擬中，可以得到傳統方法無法獲取的豐富仿真結果，體現出了較明顯的優勢。不足之處在于，在油膜結構發散段的氣穴潰滅混合和溫度計算問題上，現有空化模型仍有進一步改進的空間。

參考文獻：

[1] Hartinger M，Dumont M L，Loannides S，et al. CFD modeling of a thermal and shear-thinning elastohydrodynamic line contact[J]. Journal of Tribology，2008，130：041503.

[2] Bruyere V，Fillon N，Morales-Espejel G E，et al. Computational fluid dynamics and full elasticity model for sliding line thermal elastohydrodynamic contacts[J]. Tribology International，2012，46(1)：3-13.

[3] 于曉東,付旭,劉丹,等. 環形腔多油墊靜壓推力軸承熱變形[J]. 吉林大學學報:工學版,2015,45(2):460-465.

Yu Xiao-dong，Fu Xu，Liu Dan，et al. Thermal deformation of annular recess multi-pad hydrostatic thrust bearing[J]. Journal of Jilin University(Engineering and Technology Edition)，2015,45(2):460-465.

[4] Dowson D. A generalized Reynolds equation for fluid-film lubrication[J]. International Journal of Mechanical Sciences，1962，4(2)：159-170.

[5] Yang Pei-ran，Wen Shi-zhu. A generalized Reynolds equation for non-Newtonian thermal elastohydrodynamic lubrication[J]. Journal of Tribology，1990，112(4):631-636.

[6] 滕憲斌,張均東,葉偉強. 基于熱彈性流體動壓潤滑的柴油機曲軸主軸承潤滑特性分析[J]. 船舶工程，2016,38(4):61-64,80.

Teng Xian-bin，Zhang Jun-dong，Ye Wei-qiang. Characteristic analysis of diesel engine main bearing lubricationbased on thermo-elastic hydrodynamic pressure[J]. Ship Engineering，2016,38(4):61-64,80.

[7] Sun Jun, Zhao Xiao-yong, Wang Hu. Lubrication analysis of crankshaft bearing considering crankshaft deformation[C]∥SAE Papers，2011-01-0613.

[8] Inui M, Kobayashi M, Oowaki K,et al. Analysis of oil film generation on the main journal bearing using a thin-film sensor and elasto-hydrodynamic lubrication (EHL) model[J]. Sae International Journal of Fuels & Lubricants，2013,6(1):119-125.

[9] 桂長林,何芝仙,李震,等. 求解內燃機曲軸-軸承軸心軌跡的動力學方法[J]. 機械工程學報，2009,45(8):163-170.

Gui Chang-lin，He Zhi-xian，Li Zhen，et al. Dynamic method to calculate for shaft center trajectory of bearing of internal combustion engine[J]. Journal of Mechanical Engineering，2009,45(8):163-170.

[10] Zhang Xiu-li，Yin Zhong-wei，Gao Geng-yuan,et al. Determination of stiffness coefficients of hydrodynamic water-lubricated plain journal bearings[J]. Tribology International，2015，85：37-47.

[11] Gao Geng-yuan，Yin Zhong-wei，Jiang Dan,et al. Numerical analysis of plain journal bearing under hydrodynamic lubrication by water[J]. Tribology International，2014，75:31-38.

[12] Li Qiang, Liu Shu-lian, Pan Xiao-hong,et al. A new method for studying the 3D transient flow of misaligned journal bearings in flexible rotor-bearing systems[J]. Journal of Zhejiang University—Science A(Applied Physics & Engineering)，2012，13(4):293-310.

[13] Shi Xi, Ni Ting. Effects of groove textures on fully lubricated sliding with cavitation[J]. Tribology International，2011，44(12):2022-2028.

[14] Papadopoulos C I, Nikolakopoulos P G, Kaiktsis L. Characterization of stiffness and damping in textured sector pad micro thrust bearings using computational fluid dynamics[J]. Journal of Engineering for Gas Turbines and Power，2012,134(11):112502.

[15] 符永宏,張宗濤,張彥虎,等. 基于CFD活塞環表面凹腔微造型動壓潤滑性能研究[J]. 潤滑與密封，2010,35(12):25-27.

Fu Yong-hong，Zhang Zong-tao，ZhangYan-hu，et al. Research on hydrodynamic lubrication of piston ring surface micro-pores based on CFD[J]. Lubrication Engineering，2010,35(12):25-27.

[16] 高慶水,楊建剛. 基于cfd方法的液體動壓滑動軸承動特性研究[J]. 潤滑與密封，2008,33(9):70-72.

Gao Qing-shui，Yang Jian-gang. Research on the dynamic characteristics of hydrodynamical journal bearing based on CFD analysis[J]. Lubrication Engineering，2008,33(9):70-72.

[17] 張楚,楊建剛,郭瑞,等. 基于兩相流理論的滑動軸承流場計算分析[J]. 中國電機工程學報，2010,30(29):80-84.

Zhang Chu，Yang Jian-gang，Guo Rui，et al. Simulation of journal bearing flow field using computational fluid dynamics two phase flow theory[J]. Proceedings of the CSEE，2010,30(29):80-84.

[18] 王康,張衛正,郭良平,等. 內燃機連桿大頭軸承潤滑CFD分析[J]. 內燃機學報，2006,24(2):173-178.

Wang Kang，Zhang Wei-zheng，Guo Liang-ping，et al. CFD simulation of engine connecting-rod big-end bearing lubrication[J]. Transactions of CSICE，2006,24(2):173-178.

[19] Li Jiang，Chen Hao-sheng. Evaluation on applicability of Reynolds equation for squared transverse roughness compared to CFD[J]. Journal of Tribology，2007，129(4)：963-967.

[20] Gertzos K P，Nikolakopoulos P G，Papadopoulos C A. CFD analysis of journal bearing hydrodynamic lubrication by Bingham lubricant[J]. Tribology International，2008，41(12)：1190-1204.

[21] 孟凡明,陳原培,楊濤. CFX和Fluent在徑向滑動軸承潤滑計算中的異同探討[J]. 重慶大學學報，2013,36(1):7-14.

Meng Fan-ming，Chen Yuan-pei，Yang Tao. Discussion on similarities and differences between CFX and Fluent software in calculation journal bearing lubrication[J]. Journal of Chongqing University，2013,36(1):7-14.

[22] Deligant M，Podevin P，Descombes G. CFD model for turbocharger journal bearing performances[J]. Applied Thermal Engineering，2011，31(5)：811-819.

[23] Lin Q,Wei Z,Wang N,et al. Analysis on the lubrication performances of journal bearing system using computational fluid dynamics and fluid-structure interaction considering thermal influence and cavitation[J]. Tribology International，2013,64(3):8-15.

[24] Liu H,Xu H,Ellison P J,et al. Application of computational fluid dynamics and fluid-structure interaction method to the lubrication study of a rotor-bearing system[J]. Tribology Letters，2010,38(3):325-336.

[25] Almqvist T,Almqvist A, Larsson R. A comparison between computational fluid dynamic andreynolds approaches for simulating transient EHL line contacts[J]. Tribology International，2004,37(1):61-69.

[26] Almqvist T, Larsson R. The Navier-Stokes approach for thermal EHL line contact solutions[J]. Tribology International，2002,35(35):163-170.

[27] Almqvist T, Larsson R. Thermal transient rough EHL line contact simulations by aid of computational fluid dynamics[J]. Tribology International，2008,41(8):683-693.

[28] Singhal A K, Athavale M M, Li H, et al. Mathematical basis and validation of the full cavitation model[J]. Journal of Fluids Engineering，2002,124(3):617-624.

[29] Ferron J，Frene J，Boncompain R. A study of the thermohydrodynamic performance of a plain journal bearing comparison between theory and experiments[J]. Journal of Lubrication Technology，1983，105(3)：422-428.

[30] Pierre I，Fillon M，Bouyer J. Thermohy drodynamic behavior of misaligned plain journal bearings: Theoretical and experimental approaches[J]. Tribology Transactions，2004，47(4)：594-604.