基于動態子結構的軸承動力學及摩擦學耦合分析

楊欣，顧根香，孫思聰，周建明，李賽力

1 上海船用柴油機研究所，上海 200090

2 武漢理工大學 能源與動力工程學院，湖北 武漢 430063

0 引 言

滑動軸承被廣泛應用于艦船等各類機械系統中，作為主要的運動副和激勵源，它是功率損耗和振動的主要來源之一，動力學及摩擦學特性對其性能及可靠性影響很大[1-2]。因此，對滑動軸承摩擦學和動力學特性進行試驗與計算分析具有重要意義。滑動軸承摩擦學和動力學性能研究的3個最重要因素是運動部件的動力學特性、接觸面間的摩擦學特性以及其試驗驗證，但基于有限元法的建模往往會產生較多的自由度，而接觸面的摩擦特性也會引入非線性從而導致計算資源消耗的進一步增加。此外，針對發動機內部的軸承的試驗條件比較苛刻[3-4]。

常規的有限元離散化法由于會產生過多的自由度而導致求解十分困難，為了在保持計算精度的同時考慮油膜的非線性，近年來發展了動態子結構法，即在建立復雜結構系統力學模型時，從量級上大幅縮減整體結構的自由度，從而保證計算精度[5-7]。針對軸頸和軸瓦間的油膜作用，基于Hertz接觸模型提供了一種簡化的非線性計算方法[8]。Greenwood等[9]和Patir等[10]分別對彈性粗糙表面的粗糙接觸機理和粗糙表面間的油膜動壓潤滑機理進行研究，得到了Greenwood/Tripp粗糙接觸模型和基于平均流量系數的彈性流體動力潤滑（EHD）模型，奠定了彈性流體動力潤滑理論基礎。Bukovnik等[11]等比較了不同類型滑動軸承仿真模型和油液潤滑條件的影響。盧伯聰等[12]對主軸承熱彈性流體動力潤滑（TEHD）進行研究，并在TEHD求解中考慮了軸瓦熱效應的影響，從理論上更接近軸承實際的工作狀態，但相關研究主要集中在理論計算分析層面，缺乏試驗數據的支撐。

本文擬采用仿真與試驗相結合的方法，首先基于子結構法建立各零部件的體單元，并通過EHD模型考慮連接副的摩擦學特性，進而建立考慮滑動軸承油膜作用的多柔性體主軸承磨損試驗臺動力學模型；接著，進行動力學和摩擦學的仿真計算并與試驗結果進行對比分析，驗證計算方法的正確性，然后在此基礎上揭示軸承轉子系統的動力學與摩擦學耦合機理；最后，對比基于非線性彈簧模型和EHD模型的計算精度，進而得到主軸承熱彈性流體潤滑計算的精確建模方法。

1 基本理論

1.1 動態子結構縮減理論

采用有限元方法求解往往會產生規模巨大的離散化模型，加大了求解的難度。子結構法是通過將結構由整體劃分為簡單的子結構，然后基于特定的算法（靜態縮減、古典動態子結構法、約束模態綜合法或自由模態綜合法等）將原本的大自由度系統縮減為保留主自由度及包含質量、剛度、恢復矩陣的子結構模型。靜態縮減往往會使固有頻率偏高，且其假設低頻段慣性力對副自由度的影響比靜力效應小，因此在多體系統動力學中，采用動態縮減進行求解的精度更高，其運動方程為[12]：

式中：M，K分別為模型的質量和剛度矩陣；x，f分別為位移和力矢量；下標m，s分別為主自由度和副自由度的坐標。

由式（1）中的第2個方程式，可以得到

其解為：

取其第i階特征值λj2的值Pj進行計算，式（3）變為

式中：I為單位矩陣；Ts為坐標轉換矩陣。將式（4）代入式(1)，變換后縮減模型的運動方程為

1.2 Greenwood/Tripp粗糙接觸模型

隨著供油條件變差，當油膜厚度減小至一定程度時，兩粗糙表面間部分微凸體直接接觸，從而產生較大的摩擦力。針對這種現象，Greenwood等[9]提出了一種考慮表面粗糙形貌名義平面的彈性接觸理論，得到了接觸變形與表面形貌的相互關系。對于微凸體高度呈標準正態分布的表面，根據Greenwood/Tripp粗糙接觸理論，微凸體接觸壓力pa(h)為

式中：K為彈性因子；E'為等效彈性模量，Pa；Hs為無量綱凸峰間隙；h為潤滑油名義油膜厚度，m；為平均凸峰高度，m；σ為兩接觸表面綜合粗糙度，m；ε1，ε2分別為兩接觸面材料的泊松比；E1，E2為兩接觸面材料的彈性模量，Pa。

1.3 粗糙接觸表面平均雷諾方程

聯立流體連續性方程和運動方程（Navier-Stokes方程），假設滑油密度、粘度在油膜厚度方向為常值，可得到簡化后的一般形式為雷諾方程。為了考慮流體膜厚度達到與固體表面粗糙度接近的數量級時表面微觀形貌對潤滑性能的影響，以一般形式的雷諾方程為基礎，Patir等[10]引入了壓力流量系數、剪切流量系數概念，得到平均雷諾方程為

式中：x，z分別為軸向和徑向坐標軸；?x，?z分別為x，z向壓力流量系數；ρ為流體密度，kg/m3；η為機油粘度，mPa·s；h為油膜平均厚度，m；p為油膜壓力，Pa；θ為機油填充率；U為軸頸表面線速度，m/s；?s為剪切流量系數；t為時間，s。

1.4 EHD數學模型的計算流程

基于動態子結構彈性流體動力潤滑（EHD）的計算需聯立能量方程、平均雷諾方程、彈性位移方程及載荷平衡方程，并進行求解[8]。其中，能量方程、彈性位移方程采用有限元法求解；平均雷諾方程采用有限差分法求解，且差分網格由軸瓦內表面的縮減模型主自由度節點決定。在計算過程中，采用平均雷諾方程求解得到油膜壓力（包含粗糙接觸產生的接觸壓力）后代入彈性位移方程，從而得到各個節點的彈性位移和油膜厚度并迭代至收斂。其計算流程如圖1 所示。

圖1 EHD仿真計算流程圖Fig.1 Flow chart of EHD simulation

2 軸承磨損試驗臺EHD建模

基于動態子結構縮減理論，利用AVL_Excite Power Unit模塊，以軸承磨損試驗臺為研究對象進行軸承熱彈性流體潤滑特性分析。軸承磨損試驗臺主要包含驅動電機、傳動軸系、扭矩傳感器、被測單元以及液壓和供油系統。其中被測單元包含2個支撐軸承座（每個軸承座含2個支撐軸承）、1個被測軸及連桿；驅動電機通過傳動軸為被測單元的被測軸提供轉速邊界條件；液壓系統通過液壓缸將載荷施加于連桿；供油系統為系統各軸承提供穩定的供油邊界條件。軸承磨損試驗臺的三維結構如圖2所示，其基本參數如表1所示。

表1 軸承磨損試驗臺主要參數Table1 Main parameters of bearing wear test bed

圖2 軸承磨損試驗臺的三維結構Fig.2 Three-dimensional structure of bearing wear test bed

2.1 多體動力學建模

在實際工作中，除被測單元外，其他部件無宏觀運動，建模時，只需將被測部分的固定件和運動件視為柔性體，液壓載荷以外載荷形式輸入，右側電機的作用通過右側節點的恒定轉速進行考慮即可。試驗臺包含兩側的支撐座共4個滑動軸承以及被測軸承，因此，所建的動力學模型包含2個軸承座和被測偏心軸以及連桿體單元，其中在偏心軸與各軸承座之間建立2個EHD連接副，在偏心軸與連桿之間建立1個EHD連接副，得到的多體系統動力學模型如圖3所示。

圖3 軸承磨損試驗臺動力學模型Fig.3 Dynamic model of bearing wear test bed

2.2 縮減模型的建立

偏心軸和連桿是軸承磨損試驗臺多體系統動力學模型的主要運動部件，其動力學特性的準確性將直接影響計算結果的精度。為了提高計算效率，采用自由界面子結構法（Craig-Chang法）建立了縮減模型，如表2所示。為了保證縮減模型的動力學特性與有限元模型一致，分別對有限元模型和縮減模型進行了模態計算，其各階模態頻率的對比如表3和表4所示。經對比分析可知，偏心軸、連桿縮減體的各階模態頻率與有限元模型的誤差不超過2%，說明縮減模型與有限元模型的動力學特性基本一致。

表2 主要零部件有限元模型與縮減模型Table2 Finite element model and condensed model of main parts

表3 偏心軸有限元模型與縮減模型模態對比Table3 Comparison of modal frequency between finite element models and condensed model of eccentric shaft

表4 連桿有限元模型與縮減模型模態對比Table4 Comparison of modal frequency between finite element model and condensed model of connecting rod

3 試驗驗證及結果分析

3.1 測試系統

為了驗證計算模型的正確性，在軸承磨損試驗臺布置了相應的測點進行試驗數據采集。首先，通過扭矩儀間接測量軸承的摩擦功耗，如圖2所示。圖4（a）所示為主要測點布置圖。在連桿兩側伸出的結構上安裝夾角為90°的位移傳感器，用于測量軸心軌跡，其原理如圖4（b）所示，軸承磨損試驗臺的測量參數如表5所示。

表5 軸承磨損試驗臺主要測試參數Table5 Main test parameters of bearing wear test bed

圖4 軸承磨損試驗臺動力學模型Fig.4 Dynamic model of bearing wear test bed

3.2 計算結果驗證與分析

3.2.1 結果驗證

模型的驗證從2個方面進行：1）對比360°曲柄轉角內軸承磨損試驗臺摩擦力扭矩的計算值與試驗值，驗證其摩擦學特性；2）對比計算與試驗所得的軸心軌跡，驗證軸承磨損試驗臺模型轉子的動力學特性，如圖5所示。

圖5 軸承磨損試驗臺計算結果驗證Fig.5 Verification of calculation results of bearing wear test bed

經對比發現，試驗測得的扭矩值大于計算值，這是因為試驗臺在安裝過程中存在一定的軸系不對中情況，產生了部分阻力，且自由端儲存的滑油也會對軸產生拖拽力矩，不過在摩擦力力矩方面計算值與試驗值的誤差在10%以內；計算所得與試驗測得的軸心軌跡輪廓及位置基本一致，驗證了該計算模型在摩擦學及動力學上的正確性。

3.2.2 結果分析

運用本文第2節的建模方法，建立包含5個軸承的軸承磨損試驗臺EHD多體系統動力學模型，按X方向，從左到右,各軸承號分別為1#～5#。在驗證模型精度的基礎上，對各軸承的潤滑特性進行計算分析。表6對比列出了5個軸承在3 000 r/min工況下的最大油膜壓力pmax、最小油膜厚度hmin、最小膜厚比λ、一周期內平均粗糙接觸比例γ及摩擦功耗Wm。

表6 各軸承EHD潤滑計算結果Table6 EHD calculation results of diffrent bearings of the test bed

由表6可知，1#～5#軸承的潤滑狀態參數關于3#軸承（被測軸承）對稱，被測軸承承受了最大的載荷，且最大油膜壓力較其他4個支撐軸承高出一個量級，最小油膜厚度達到了綜合粗糙度量級，使得膜厚比偏小，出現了局部粗糙接觸。根據Stribeck曲線[13]，兩側支撐軸承最小膜厚比大于10，處于完全液動潤滑狀態，而被測軸承則處于混合摩擦狀態，軸頸與軸瓦間的粗糙接觸比例和摩擦功耗呈現相同的規律。被測軸承最小油膜厚度出現的時刻與活塞上的止點相對應，此時，倒置的連桿處于最高位置，下軸瓦與偏心軸表面接觸，液壓系統對連桿施加的載荷處于峰值。由于支撐軸承的數量為被測軸承的4倍，單個軸承的軸瓦承受的載荷遠小于被測軸承下瓦，此時，被測軸承下瓦的油膜壓力遠大于支撐軸承的油膜壓力，且軸頸與軸瓦間出現了部分粗糙接觸，產生了粗糙接觸壓力，進而使軸瓦進入混合潤滑狀態，粗糙接觸比例增高，摩擦功耗增高。

分別提取5個軸承在被測軸承產生最小油膜厚度時刻的油膜總壓和粗糙接觸壓力，其分布圖分別如圖6和圖7所示。由圖中不難發現，最外側固定支撐軸承的油膜壓力的分布梯度要小于中間支撐軸瓦的油膜壓力分布梯度。另由圖7所示粗糙接觸壓力分布彩圖，發現支撐軸瓦無粗糙接觸壓力，說明即便是處于最惡劣的工況，支撐軸瓦仍可保持熱動潤滑狀態，而被測軸瓦下瓦的兩側則出現了局部粗糙接觸。

圖6 各軸承油膜總壓分布彩圖Fig.6 Oil film total pressure distribution of different bearing

圖7 各軸承油膜粗糙接觸壓力分布彩圖Fig.7 Oil film asperity contact pressure distribution of different bearings

由上述分析可知，被測的3#軸承其潤滑條件最差。為了解其潤滑特性，分析不同間隙被測軸承的動力學及摩擦學特性，對比了不同間隙（1～1 000 μm）下，一個周期內的最大油膜壓力、最小油膜厚度和平均摩擦功耗。

由圖8（a）和圖8（b）可以看出，隨著半徑間隙的增大，軸承最大油膜壓力隨之升高，但最小油膜厚度在一定范圍內是先增大后減小，當半徑間隙達20 μm時，最小油膜厚度達到峰值，摩擦功耗最低。當半徑間隙非常小時，軸頸與軸瓦間出現微凸體接觸，軸承處于邊界潤滑狀態，不易形成足夠厚的油膜，而隨著間隙的增大，最小油膜厚度逐漸增大，油膜支撐力增強，但當間隙過大時，在端泄的作用下，油膜支撐能力將下降，最小油膜厚度也隨之下降。同理，摩擦功耗在半徑間隙非常小的時候因邊界潤滑摩擦系數較大，隨著間隙的增大逐漸處于液動潤滑狀態，摩擦力下降；而當間隙過大時，液動潤滑剪切應力會升高，同時由于油膜支撐力下降，局部出現的粗糙接觸也會導致摩擦功耗上升。

圖8 潤滑特性隨間隙變化的趨勢Fig.8 Variation of lubrication characteristics with clearance

間隙為10和1 000 μm這2種工況下的軸承總摩擦功耗及軸承載荷如圖9所示。由圖可知，兩種工況下的軸承載荷比較接近，但在載荷極值處，即液壓系統載荷峰值處，后者較前者力的幅度增大了幾百牛，總摩擦耗功也同樣從10 μm時的0 W上升到了幾百瓦，且摩擦功耗峰值點與載荷極值點對應，說明間隙過大產生的沖擊導致局部潤滑狀態突變，進而使摩擦功耗增加。

圖9 一周期內摩擦功耗和載荷Fig.9 Friction power consumption and load in one cycle

3.3 不同建模方式計算結果對比

由上節分析可知，軸承動力學與摩擦學之間存在較強的耦合關系。除上節敘述的通過求解平均雷諾方程進行軸承動力學分析之外，軸承油膜的建模還有其他常用方法，比如簡化的非線性彈簧單元，即將軸頸和軸瓦間的高壓油膜簡化為軸頸與軸瓦間的若干非線性彈簧阻尼單元，求解時，通過激活處于壓縮狀態的彈簧來迭代計算得到相應的彈簧力，進而得到軸承載荷。為了對比基于EHD模型和非線性彈簧模型建模方法的差別，在相同動力學模型的基礎上修改被測軸承連接副的單元類型為非線性彈簧單元并進行計算。將非線性彈簧單元與不同間隙EHD連接副算得的被測試驗臺軸承動力學、摩擦學結果進行對比，其時域圖如圖9（b）所示。由圖中不難發現，當軸承間隙較小時，采用非線性彈簧模型與EHD模型計算得到的軸承載荷基本一致，但隨著軸承間隙的增大，油膜因端泄的作用其支撐力下降；軸頸與軸瓦沖擊的增大在載荷時域圖中表現為明顯的波動。非線性彈簧不存在實際工況中可能出現的沖擊現象，采用計算得到的結果在間隙較大的情況下依然是較順滑的載荷曲線，這與實際工況差異較大。對采用2種建模方法計算得到的軸承載荷進行了傅里葉變換，結果如圖10所示。由圖可知，這2類模型在低頻段重合較好，但采用彈簧單元計算得到的高頻分量幅值較EHD單元要低，故對于軸承動力學的求解，使用EHD模型可以提高振動響應的求解精度。

圖10 2種建模方式下載荷頻域力的對比Fig.10 Comparison of load force between two modeling methods in frequency domain

4 結 論

本文以軸承磨損試驗臺為研究對象，基于子結構法建立了多體系統動力學模型，然后采用仿真計算與試驗相結合的方法研究了滑動軸承摩擦學與動力學的耦合關系，主要得到如下結論：

1） 采用子結構法建立的軸承磨損試驗臺多體系統動力學模型，通過對動力學和摩擦學分別進行的試驗驗證，證實了建模方法的正確性。

2） 在試驗工況下，軸承磨損試驗臺的4個支撐軸承處于完全液動潤滑狀態，被測軸承處于混合潤滑狀態，這種惡劣的潤滑條件主要是因上止點的載荷較高且受力集中所致。

3） 隨著間隙的增大，被測軸承的潤滑狀態先改善并在20 μm間隙時達到最佳，隨后惡化；隨著間隙的增大，最大壓力隨之增大，最小油膜厚度先增加后減小，摩擦功耗則是先減小后增大，因此在進行軸瓦設計時，需確定合適的間隙以同時滿足潤滑性能需求。

4） 當間隙較小時，采用非線性彈簧可以得到較高精度的軸承載荷，但當軸承間隙較大時，在外載的作用下軸頸與軸瓦之間存在沖擊，采用EHD模型的求解精度更高。