滑動軸承壓力分布及動特性系數

史冬巖，張成，任龍龍，張亮，彭梁

(哈爾濱工程大學 機電工程學院，黑龍江 哈爾濱 150001)

滑動軸承在旋轉機械中應用廣泛，其工作時轉子與軸承之間形成壓力油膜.油膜壓力作用在軸瓦的合金層上，循環交變的應力是導致軸瓦變形失效的主要原因.對油膜壓力的計算和軸瓦合金層應力分布的研究是對滑動軸承進行設計和失效分析的重要理論依據.油膜不僅起著承受載荷、減輕摩擦、消除磨損等作用，從動力學觀點看來，油膜的動特性對整個轉子系統的動力特性有很大影響［1-3］.它也是轉子—支承—基礎這個系統中的一個環節［4］.

本論文以流體動壓潤滑滑動軸承為研究對象，以二維流動的雷諾方程作為研究的基礎.對雷諾方程的求解是滑動軸承研究的關鍵問題之一，早期對雷諾方程的求解是基于解析解法，但解析解法無法獲得較精確的解，隨著計算機技術的發展.數值計算方法對雷諾方程求解己經成為一個主流趨勢.有限差分法編程簡單，求解方便，本文采用有限差分法求解雷諾方程［5-6］，在求解得出油膜壓力分布的基礎上研究軸瓦合金層的應力分布情況，并對滑動軸承的動態特性進行研究，求解出滑動軸承的動特性系數，為有限元分析中滑動軸承的簡化提供幫助.

1 差分法求解壓力分布

從層流運動的油膜中取一個微小的單元體作為研究對象，可以導出雷諾方程的一般形式［7］:

式中:x為周向坐標，z為軸向坐標，p為油膜壓力，h為油膜厚度，U為軸頸速度，μ為油膜粘度.為了將方程寫成最緊湊的形式，將式(1)無量綱化，可得雷諾方程的無量綱形式［8］為

式中:φ為偏位角，l為軸承寬度，λ=2z/l，H=1+εcosφ，H為無量綱油膜厚度，ε為偏心率.P為無量綱油膜壓力，由式(2)可以得出，無量綱油膜壓力P的分布取決于偏心率ε和寬徑比d/l.

本文選擇采用有限差分法求解雷諾方程.將軸瓦的油膜區域劃分為網格，如圖1所示，用各個節點上的壓力值構成各階差商，近似取代節點上的壓力值.所得的一組離散的壓力數值，也就近似表達了油膜中的壓力分布.先把整個油膜區域離散成長方形的網格，將網格節點按所在的列數和行數順序編號，沿φ方向的列數用i編號，沿λ方向的列數用j編號，每個節點位置用(i，j)二維編號表示［9］.

圖1 網格劃分及差商示意Fig.1 Schematic diagram of mesh generation and difference quotient

對節點(i，j)上的一階導數，可用其相鄰節點上的P值構成的中差商來表達，為了提高計算精度，采用半步長上的值構成的中差商表示一階導數，對于(i，j)上的二階導數，可先用相鄰步長上的一階導數的中差商表示，然后將式中的一階導數用相鄰節點值的中差商表示，則式(2)可表示為

式(3)適用于全部內節點 i=2，3，…，m、j=2，3，…，n，共有 (m-1)(n-1)個方程，可構成一個方程組，根據給定邊界條件可解出各內節點Pi，j值.

引入雷諾邊界條件最常用的有效且簡便的做法是:在網格區域每行上均由起始邊向終止邊方向逐點計算，如果算出某點壓力為負，即取為零.此點位置即可作為該行上油膜自然破裂邊的近似位置.該點以后各點壓力均取為零，不按式(3)計算.每次迭代均如此處理，則破裂邊近似位置會逐漸逼近應有的自然破裂邊界，整個壓力分布也就逼近計入雷諾邊界條件的壓力分布.

軸承參數取自文獻［10］的一個計算示例，該示例中數據如表1所示.計入表1數據，取m=60，n=40得出一組無量綱油膜壓力分布，其分布情況如圖2所示.

表1 軸承參數Table 1 Bearing parameter

圖2 油膜壓力分布Fig.2 Pressure distribution of oil film

從壓力分布圖2中可以看出，360°包角有限寬徑向軸承的無量綱油膜壓力的分布為近似拋物面分布.無量綱油膜壓力在某一段逐漸增大到最大壓力值，之后急劇下降，在φ＞180°的某一區域，壓力降為零，壓力變為零的點就是油膜的自然破裂點.在油膜壓力增大到最大的過程中，油膜壓力變化平滑，在超過峰值后，油膜壓力變化幾乎突變為零，這與實際情況中楔形油膜的變化情況一致.

2 合金層應力應變分布求解

根據表1數據建立滑動軸承網格模型，共生成單元數129 840個、節點數105 652個，外層鋼材料采用五面體單元進行離散，內層的合金層采用六面體單元進行離散.鋼被層和合金層材料屬性如表2所示.

表2 軸承材料參數Table 2 Bearing material parameter

網格模型導入到ANSYS中，在前面編程求得的油膜壓力分布式大量的無規則的、離散的值，選擇用APDL命令進行加載.油膜壓力加載情況如圖3所示.求解得出的合金層應力應變分布如圖4～6所示.

圖3 油膜壓力加載Fig.3 Oil film pressure loading

圖4 滑動軸承軸瓦應力分布Fig.4 Stress distribution of sliding bearing bush

圖5 滑動軸承軸瓦剪應力分布Fig.5 Shear stress distribution of sliding bearing bush

圖6 滑動軸承軸瓦應變分布Fig.6 Strain distribution of sliding bearing bush

由圖4合金層應力分布情況，可以看出，滑動軸承在承受油膜壓力所產生的應力的分布與油膜壓力的分布基本相同.在滑動軸承的寬度方向，應力從外截面到中截面逐漸增加，合金層軸向應力的峰值位于中截面油膜厚度最小處［11］，峰值為20.702 MPa.

由圖5滑動軸承軸瓦剪應力分布情況，可以看出，剪應力存在的區域為壓力峰值周圍和壓力梯度大的區域.同時，油膜壓力最大值在滑動軸承的中截面處并具有最大壓力梯度，剪應力的最大值在滑動軸承的中截面處合金層與鋼被結合處，因為該處壓力梯度最大，最大剪應力為8.375 MPa.在軸承表面，剪應力最大值在油膜壓力分布峰值處，最大值為5.322 MPa.

由圖6滑動軸承的應變分布情況，可以看出，滑動軸承在油膜壓力作用下的徑向變形與油膜壓力分布十分相似，在一定區域內，變形隨著油膜壓力的增大而逐漸變大，當油膜壓力達到最大值時，變形也達到最大.隨著油膜壓力急劇降低而使變形變小，應變隨著油膜厚度的增大而逐漸減小.

3 滑動軸承動力特性系數計算

油膜不僅起著承受載荷、減輕摩擦、消除磨損等作用，從動力學觀點看來，它也是轉子—支承—基礎這個系統中的一個環節.油膜通常起著非線性的彈簧和阻尼作用.由于這類問題大多只涉及到較小振幅或無限小的振幅，所以常可將油膜近似看成具有線性化了的彈簧常數和阻尼特性，通常即稱這些線性化了的動力特性為油膜剛度和阻尼，下面介紹油膜剛度和阻尼系數［7］的求解.

3.1 等效力學模型

定義油膜剛度系數為單位位移所引起的油膜力增量，即

定義油膜阻尼系數為單位速度所引起的油膜力增量，即

式中，各系數的第1個下標代表力的方向，第2個下標代表位移或速度的方向.油膜剛度系數和阻尼系數統稱為油膜動力特性系數或動力系數.動態分析時滑動軸承的力學模型可表示成如圖7所示情況［12］.

圖7 滑動軸承等效力學模型Fig.7 Equivalent mechanical model of sliding bearing

3.2 擾動壓力方程求解

徑向滑動軸承的非定常運動雷諾方程為

式中:h為油膜厚度，p為油膜壓力，μ為潤滑油膜動力粘度，φ為軸頸軸向位置，z為軸向位置，r為軸頸半徑，Ω是軸頸角速度，t為時間.

設Δx為橫向擾動、Δy為垂向擾動、h0為靜平衡油膜厚度 ，將(6)式按 h=h0+ Δxsin φ -Δycos φ展開，并忽略高階小量，可以獲得各項擾動壓力的微分方程為［13］

式中，P'分別表示Px、Py、P˙x、P˙y、P0為靜平衡時的油膜壓力.第1個方程為Δx對應項，第2個方程為Δy對應項，第3個方程為Δ˙x對應項，第4個方程為Δ˙y對應項.各擾動方程與無擾動下雷諾方程形式一樣，只是右邊項不同，故求解方法和求解無擾動下雷諾方程一樣，在求解前需把擾動方程無量綱化.令

式中:ψ =c/r.

將式(8)中各無量綱因子帶入式(7)即可得出擾動壓力微分方程的無量綱形式如下:

式中:Pi分別表示 P1、P2、P3、P4的無量綱形式.

式(9)表示的是無量綱擾動壓力與靜壓力分布、油膜厚度、周向位置之間的關系.由式(9)中第3個、第4個方程可以看出動態速度引起的油膜壓力增量P3、P4與靜平衡位置時油膜壓力分布情況和油膜厚度分布情況無關.

3.3 無量綱動特性系數的計算

求解擾動壓力即解方程(9)，求解方法和求解無擾動下雷諾方程一樣.對于這些擾動壓力，邊界條件是:在完整油膜區的全部邊界上，這些擾動壓力均為零.在計算時，先按求解出軸心靜平衡位置相應的壓力分布，以及由雷諾邊界條件確定的破裂邊位置，然后按方程(9)計算各擾動壓力.當計算出各擾動壓力 P1、P2、P3、P4后，經再次積分求出油膜力的各項增量，即無量綱剛度系數如式(10)所示，無量綱阻尼系數如式(11)所示:

求解得出無量綱油膜動特性系數如表3所示.由表3可以看出，滑動軸承油膜無量綱剛度和阻尼系數隨著偏心率的增大而增大，在偏心率小于0.5時，隨偏心率的變化不大;當偏心率大于0.5后，隨偏心率的增大而迅速增大;當偏心率超過0.8時，該趨勢變的更加明顯.另外可以看出，交叉阻尼系數近似相等.這一點也可以有力的證明前面所采用的算法是可行的，因為由滑動軸承動態特性的經典理論知道滑動軸承油膜交叉阻尼是相等的，在本文中由于采用數值計算方法，二者近似相等.

表3 無量綱油膜動特性系數(L/D=0.5)Table 3 Dynamic characteristics of oil film with dimensionless(L/D=0.5)

4 結論

本文對滑動軸承的油膜壓力和軸瓦合金層應力分布進行了研究，在求解出油膜壓力分布的基礎上對滑動軸承的動態特性進行了研究，求解出滑動軸承的動特性系數，主要得出了以下結論:

1)由油膜壓力分布圖可得出滑動軸承的壓力分布為近似拋物面分布.無量綱油膜壓力在某一段逐漸增大到最大壓力值，之后急劇下降，在φ＞180°的某一區域，壓力降為零.

2)滑動軸承合金層在承受油膜壓力所產生的應力分布與油膜壓力的分布基本一致，應力峰值為20.702 MPa.剪應力存在的區域為壓力峰值周圍和壓力梯度大的區域，最大剪應力為8.375 MPa.

3)滑動軸承油膜無量綱剛度系數和阻尼系數隨著偏心率的增大而增大，交叉阻尼系數近似相等.

［1］姚熊亮，孫士麗，陳玉.高頻動載軸承內油膜壓力特性［J］.機械工程學報，2010，46(17):93-99.

YAO Xiongliang，SUN Shili，CHEN Yu.The pressure behavior of oil film in bearing subjected to high-frequency dynamic load［J］.Journal of Mechanical Engineering，2010，46(17):93-99.

［2］鄧玫，孫軍，符永紅，等.計及軸受載變形的粗糙表面軸承熱彈性流體動力潤滑分析［J］.機械工程學報，2010，46(15):95-101.

DENG Mei，SUN Jun，FU Yonghong，et al.Thermoelastohydrodynamic lubrication analysis of bearing considering shaft deformation and surface roughness［J］.Journal of Mechanical Engineering，2010，46(15):95-101.

［3］姚熊亮，張成，孫士麗.考慮可壓縮性及慣性力的油膜力研究［J］.中國艦船研究，2010，5(6)33-40.

YAO Xiongliang，ZHANG Cheng，SUN Shili.Analysis of oil film force considering compressibility and inertial force［J］.Chinese Journal of Ship Research，2010，5(6):33-40.

［4］全永昕.工程摩擦學［M］.杭州:浙江大學出版社，1994:222-250.

QUAN Yongxin.Engineering tribology［M］.Hangzhou:Zhejiang University Press，1994:222-250.

［5］王兆伍，楊家富，徐尚賢.流體潤滑軸承靜動態特性的有限分析法［J］.南京林業大學學報，1996，20(4):32-42.

WANG Zhaowu，YANG Jiafu，XU Shangxian.The finite analytic method on static and dynamic performance of journal bearings［J］.Journal of Nanjing Forestry University，1996，20(4):32-42.

［6］高創寬，齊秀梅.雷諾方程數值解中的幾個問題［J］.太原重型機械學院學報，1993，18(6):30-41.

GAO Chuangkuan，QI Xiumei.Several problems in the numerial solution to Reynolds'equation［J］.Journal of Taiyuan University of Science and Technology，1993，18(6):30-41.

［7］張直明.滑動軸承的流體動力潤滑理論［M］.北京:高等教育出版社，1986:1-5，34-46，65-72.

ZHANG Zhiming.Hydrodynamic lubrication theory of sliding bearing［M］.Beijing:Higher Education Press，1986:1-5，34-46，65-72.

［8］黃民毅.液體動壓滑動軸承設計的數值計算［J］.四川工業學院學報，1998，17(3):37-41.

HUANG Minyi.Numerical calculation in the design of sliding bearings with liquid dynamic lubrication［J］.Sichuan U-niversity of Science and Technology，1998，17(3):37-41.

［9］王寧.基于MATLAB的滑動軸承壓力分布的數值計算［D］.大連:大連理工大學，2006:1-3.

WANG Ning.Numerical Calculation to the pressure distribution of journal bearing based on the Matlab［D］.Dalian:Dalian University of Technology，2006:1-3.

［10］機械設計手冊委員會.機械設計手冊(第三卷)［M］.北京:機械工業出版社，2004:21-25.

Standard Handbook of Machine Design Committee.Standard handbook of machine design-Ⅲ［M］.Beijing:China Machine Press，2004:21-25.

［11］唐倩，方志勇，朱才朝，等.滑動軸承油膜壓力及合金層應力分布［J］.中南大學學報，2008，39(4):776-779.

TANG Qian，FANG Zhiyong，ZHU Caichao，et al.Oil film pressure and stress distribution in alloy layer of journal bearing［J］.J Cent South Univ，2008，39(4):776-779.

［12］曹樹謙，丁千，陳予恕，等.具有滑動軸承的穩態轉子系統有限元建模分析［J］.汽輪機技術，1999，41(6):347-350.

CAO Shuqian，DING Qian，CHEN Yushu，et al.Analysis on modeling steady rotor system with sliding bearings by using FEM［J］.Turbine Technology，1999，41(6):347-350.

［13］鐘一鄂.轉子動力學［M］.北京:清華大學出版社，1987:41-50.

ZHONG Yi’e.Rotor dynamics［M］.Beijing:Tsinghua U-niversity Press，1987:41-50.

［14］SUN Jun，GUI Changlin.Hydrodynamic lubrication analysis of journal bearing considering misalignment caused by shaft deformation［J］.Tribology International，2004，37:841-848.

［15］孫軍，王震華，桂長林.計入曲軸受載變形的粗糙表面曲軸軸承彈性流體動力潤滑分析［J］.機械工程學報，2009，45(1):135-140.

SUN Jun，WANG Zhenhua，GUI Changlin.Elastohydrodynamic lubrication analysis of crankshaft bearing considering crankshaft deformation under load and roughness surface［J］.Journal of Mechanical Engineering，2009，45(1):135-140.

［16］SUN Jun，GUI Changlin.Effect of lubrication status of bearing on crankshaft strength［J］.Journal of Tribology，Transactions of the ASME，2007，129:887-894.

［17］柳江林，孫軍，桂長林，等.軸頸傾斜軸承的熱流體動力潤滑分析［J］.潤滑與密封，2007(9):60-63，32.

LIU Jianglin，SUN Jun，GUI Changlin，et al.Thermohydrodynamic lubrication analyses of misaligned journal bearing［J］.Lubrication Engineering，2007(9):60-63，32.