滾動軸承運動-接觸精確建模與二次開發

馬帥軍,閆柯,張曉紅,王明凱,朱永生,洪軍

(西安交通大學現代設計及轉子軸承系統教育部重點實驗室,710049,西安)

滾動軸承作為旋轉機械轉子系統的核心支撐零部件,其動力學性能直接影響裝備整機的運行狀態。軸承力學的發展經歷了靜力學、擬靜力學、擬動力學和動力學四個階段,其中動力學階段考慮因素最多,計算精度最高,最為完善[1]。然而,基于運動微分方程對軸承動力學特性進行描述時,需要對軸承各組件進行復雜的空間力系分析,且存在方程參數多、初值依賴性強、非線性弱、收斂性不強等問題,例如廣為所知的動力學軟件ADORE[2],在仿真過程中需要設置配合公差、摩擦拖動系數、疲勞系數等上百個參數,操作專業性極強,給工程應用帶來一定不便。隨著多體動力學軟件的不斷成熟,如ADAMS、Simpack、RecurDyn等,通過將軸承各個零件視為純剛體,忽略大變形,細化小變形,大幅度提升了動力學分析軟件的計算效率[3],成為當前滾動軸承工程應用分析中的首選方案。

國內學者基于ADAMS等軟件,圍繞滾動軸承動態特性開展了大量研究工作。關猛等利用Pro-Engineer軟件建立雙列圓錐滾子軸承的三維模型,導入到ADAMS軟件進行相關參數設置,進而開展了軸承的動力學特性仿真分析[4];姚廷強等直接利用ADAMS軟件中的CMD參數化命令建立了四點角接觸球軸承三維模型,并利用宏命令進行了相關參數設置,分析了不同工況下軸承的動態性能[5];張軍飛等對上述仿真步驟進行了集成,編寫了簡單的軸承仿真界面,包括參數化建模和參數設置功能,提升了軸承動力學分析的便利性[6]。考慮到滾動軸承組件間的點線高副接觸特征,其接觸變形直接影響著套圈與滾子之間的接觸載荷與油膜壓力分布特點,進而影響軸承磨損、發熱、疲勞等一系列服役性能預測的準確性[7],因此,在利用ADAMS軟件開展軸承動力學仿真時,解決軸承套圈與滾動體間的接觸問題顯得尤為重要。ADAMS軟件在處理接觸問題時,有兩種庫可供選擇:Default_Library和Parasoilds,兩者之間的主要區別如圖1所示。其中Dafault_Library對體的識別精度較低,適合分析宏觀接觸,但計算效率高;Parasoilds對體的識別精度較高,適合分析微觀接觸,但計算效率低,很難在工程中應用[8]。

(a)Default_Library庫

(b)Parasoilds庫

針對上述問題,李琦利用Fortran語言編寫了Gfosub作用力子程序,對ADAMS進行了二次開發,開展了軸承動力學性能研究,但仍然沒有完全脫離軟件自身的接觸設置,仿真分析精度較低[9];鄧四二等利用Fortran語言編寫了Gfosub作用力子程序,對ADAMS進行了二次開發,沒有依賴于軟件自身接觸,計算精度較高[10],但是建立的動力學模型是基于局部坐標系的,在處理形狀較為規則的保持架與滾動體時效果較好[11],而對于形狀復雜的異型保持架、修型滾子等稍顯不足。

針對ADAMS軟件自身接觸求解器的問題,本文基于軸承全局坐標系,利用C語言編寫Gfosub作用力子程序,以角接觸球軸承為例,建立了軸承動力學模型,編寫的子程序與軟件自身的數據傳遞原理如圖2所示,針對Gfosub子程序及ADAMS求解器的特點,選取GSTIFF I3求解方法來求解建立的動力學模型。通過將所建立模型的仿真結果與ADAMS自身接觸模型、經典的軸承力學模型結果對比,證明了本模型的接觸仿真計算精度較高。最后搭建保持架轉速測量實驗臺,驗證本模型的運動仿真分析精度。

1 動力學模型的建立

角接觸球軸承動力學模型主要由各個零件之間的相對位置關系以及滾動體、保持架、內外套圈之間的作用力構成。對模型作如下假設。

(1)軸承各個零件的質心與其幾何中心重合。

(2)軸承各個零件均為剛體,不考慮軸承零件的大變形。

圖2 ADAMS二次開發原理示意圖Fig.2 A schematic diagram of ADAMS secondary development

(3)外圈固定,限制6個自由度;內圈由于需要添加驅動轉速,故具有5個自由度;滾動體具有6個自由度;保持架具有3個自由度。

1.1 軸承中各個參考坐標系的建立

為了能夠準確地確定軸承各個零件之間的幾何位置關系,建立角接觸球軸承的整體坐標系與各個局部坐標系。各零件坐標系如圖3所示。

圖3 角接觸球軸承的各零件坐標系Fig.3 Coordinate system for each part of angular contact ball bearing

(1)Oxyz為軸承整體的慣性坐標系,O點與外圈的質心相重合(假設外圈固定,內圈旋轉),x軸為軸承的軸向,y軸垂直于x軸指向1號滾動體,z軸方向可根據右手定則確定。

(2)Oixiyizi為內圈的隨體坐標系,Oi點與內圈的質心重合。

(3)Oajxajyajzaj為滾動體的隨體坐標系,Oaj點與滾動體j的質心重合。

(4)Objxbjybjzbj為滾動體的方位坐標系,Obj點與滾動體j的質心重合,xbj軸時刻指向軸承的軸向,ybj軸垂直于xbj徑向外指,zbj軸由右手定則確定其方向。

(5)Ocxcyczc為保持架的中心坐標系,Oc點與保持架整體的形心重合。

(6)Ocjxcjycjzcj為保持架的兜孔坐標系,Ocj點與兜孔中心重合。

1.2 軸承各個零件之間的相互作用力

(1)

式中:Toi為內圈連體坐標系到軸承整體坐標系的旋轉變換矩陣;r2為滾動體中心與內圈中心的距離;?j為滾動體j相對于軸承內圈坐標系下的方位角。

(2)

式中:ri為內圈的溝曲率半徑;Dw為滾動體直徑。

圖4 各個零件之間的相互位置關系Fig.4 Position relationship between parts

雖然滾動體與滾道之間被彈流油膜分離,但對其接觸參數影響不大,仍可用Hertz理論進行接觸力的計算。第j個滾動體與內滾道之間的法向接觸力為

(3)

式中:Kij為第j個滾動體與內滾道間的Hertz接觸剛度,具體計算詳見文獻[12]。

同理可得第j個滾動體與外滾道之間的法向接觸力

(4)

1.2.2 滾動體與套圈之間的摩擦力 在滾動體與內圈滾道的橢圓接觸面上,由于接觸點處兩者的滑動速度不同,故存在相對運動趨勢或者相對運動,導致滾動體與內外滾道表面產生油膜拖動力。第j個滾動體與內圈滾道表面的拖動力Tix(y)j可由滾動體與滾道之間的接觸應力和拖動系數的乘積進行積分得到,即為

(5)

同理得第j個滾動體與外滾道之間的拖動力為

(6)

(7)

式中:Dcm為保持架兜孔中心所在的圓的直徑。

圖5 保持架與滾動體之間的幾何關系Fig.5 Geometric relationship between cage and roller

(8)

式中:Cp為保持架與滾動體之間的兜孔間隙。

保持架的兜孔和滾動體接觸的模型可視為Hertz接觸,接觸剛度Kcj可參照套圈與滾動體的接觸剛度計算方法,則兩者之間的接觸力為

(9)

1.2.4 保持架與引導套圈之間的相互作用力 保持架與引導套圈之間的相互作用是由于潤滑劑的流體動壓效應所產生的,而套圈引導面和保持架端面的作用力比較小,故可采用短滑動軸承模型[15]。保持架與引導套圈的相關幾何關系如圖6所示。

圖6 保持架與引導套圈之間受力示意Fig.6 Force diagram between cage and guide ring

在圖6中,平面坐標系Sc=(Oc,yc,zc)固定在保持架上,yc軸指向由于流體動壓效應產生的最小油膜厚度h0所在點。以外引導的保持架為例,某一時刻保持架所受的力與力矩為

(10)

式中:η0為一般環境下潤滑油的動力黏度;u1為潤滑油的拖動速度;R1、L分別為保持架定心表面半徑和寬度;C1為保持架的引導間隙;ε為保持架質心與中心的偏移量。

為了較為方便地建立保持架的動力學微分方程,將保持架所受的力轉換到軸承的整體坐標系中,其作用力可表示為

(11)

式中:Tc為保持架坐標系到整體坐標系的變換矩陣。

1.2.5 滾動體的潤滑黏性阻力 在軸承運行過程中,軸承腔內會存在大量的潤滑劑以及空氣組成的混合物,容易對滾動體產生相關的黏性阻力,對軸承性能產生影響。但是由于內部幾何形狀、運行環境都比較復雜,建立準確的黏性阻力計算模型比較困難。在此將滾動體受到的黏性阻力等同于圓柱體在流體中移動所受的阻力[16-17],故可按下式計算:

(12)

式中:wmj為滾動體公轉速度;Cd為圓柱體的黏性阻力系數;ρ0為潤滑油密度;dm為軸承的節圓直徑;S為圓柱體的平移方向面積。

1.3 微分方程的建立

(1)對保持架進行受力分析,如圖7所示,即可得到保持架的運動微分方程

(13)

式中:mc為保持架質量;Ic為保持架轉動慣量;yc、zc為保持架位移;wcx為保持架公轉角速度。

圖7 某時刻保持架的受力情況Fig.7 The force of the cage at a certain time

(2)對內圈進行受力分析,如圖8所示,即可得到內圈的運動微分方程

(14)

圖8 某時刻內圈的受力情況Fig.8 The force on the inner ring at a certain moment

式中:mi為內圈的質量;FX、FY、FZ、MY、MZ為施加在內圈上的載荷;Iiy、Iiz為內圈在各個方向的轉動慣量。

2 動力學模型在動力學軟件中的實現

2.1 作用力子程序的開發

ADAMS軟件不僅提供了易于操作的用戶界面以及強大的動力學求解功能,還留有特定的接口,方便用戶針對特殊需求基于軟件進行二次開發。然而軸承各個零件之間的受力較為復雜,需要以子程序的形式來定義零件之間的相互作用力。在開發軸承作用力子程序時,采用ADAMS軟件提供的用戶作用力子程序的模板進行編寫,這個子程序需要包含一個頭文件slv_c_utils.h,用于預定義在子程序中用于子程序和ADAMS軟件自身的數據傳輸函數[18-19]。在進行子程序編寫時,仿真過程中調用的作用力子程序中的函數名稱、形參個數和類型必須按照ADAMS預先規定的格式進行編寫。在子程序中通過SYSARY函數或SYSFNC函數提取求解器中的狀態變量,并且利用RESULT數組函數將計算的結果返回到求解器之中。

2.2 求解器的選擇

基于Gfosub編寫完成作用力子程序之后,需要在ADAMS軟件建立微分方程,并完成求解工作。而在ADAMS求解器對微分方程進行求解時,主要包含以下3個步驟:

(1)利用顯式求解方法對微分方程進行初步求解,得到微分方程解的預估值;

(2)利用隱式求解方法再次對微分方程進行求解,得到微分方程解的校正值;

(3)將預估值與校正值之間的誤差與求解器設定誤差進行對比,若滿足設定的誤差要求,則作為微分方程的解,若不滿足,則縮小步長并返回步驟1,繼續求解。

在求解過程中,ADAMS一共提供5種類型的預估-校正方法:GSTIFF、WSTIFF、HASTIFF、HHT、Newmark,各種方法的具體特點如表1所示。其中GSTIFF和WSTIFF的I3積分格式求解器的特點是將誤差收斂在位移上,GSTIFF、WSTIFF的I2積分格式和HASTIFF的SI1/SI2求解器的特點是將誤差收斂在速度上,HHT和Newmark求解器是將誤差收斂在加速度上[20]。基于軸承的相關基礎理論,第一種仿真計算精度最高,第三種仿真計算精度最差,在此選擇GSTIFF和WSTIFF的I3積分格式求解器。

表1 ADAMS中各種求解方法的特點

2.3 Gfosub子程序與自身接觸設置對比

為了驗證本文所提出的建模方法的優越性,與常規直接應用ADAMS軟件自身的接觸算法(Dafault_Library庫)進行求解的模型進行對比分析。以7008C軸承為例,利用兩種方法建立相應的仿真模型,設置相同的仿真工況:ni=4 000 r/min,Fa=450 N,Fr=500 N。選取模型中1號滾動體與內、外圈的接觸力計算結果進行對比,結果如圖9和圖10所示。

圖9 Gfosub子程序建立模型仿真接觸力Fig.9 The force simulated by the Gfosub model

圖10 ADAMS自身接觸建立模型仿真接觸力Fig.10 The force simulated by the model of ADAMS’s contact

由于ADAMS的Parasoilds求解器不能處理復雜形狀的接觸問題,而自身默認的求解器又不適合處理微觀接觸,誤差較大,其仿真的接觸力擴大十幾倍,造成其摩擦拖動力也相應擴大,無法準確反映輕載、高速下軸承的打滑率情況,而且接觸力波動較大;而通過對ADAMS進行二次開發,利用Gfosub子程序建立的軸承動力學模型的仿真精度較高,接觸力的最大值為141.59 N,利用相關經驗公式進行預估[21],得到理論接觸力為142.85 N,兩者誤差為0.82%,在可接受范圍內,從而證明了子程序計算的優越性與可行性。

圖11 保持架轉速測量驗證實驗臺Fig.11 A verification experiment platform of cage rotation speed

3 實驗驗證

3.1 轉速測量實驗臺的搭建

為了驗證所建立的軸承動力學模型的正確性,搭建了高速滾動軸承保持架轉速測量實驗臺,如圖11所示。整體結構由機械主軸-軸承系統和測試系統兩部分構成,其中測試系統主要包括用于發射和接收紅外信號的激光轉速傳感器、用于反射激光信號的反光條、用于采集脈沖信號的數據采集系統以及用于信號處理的計算機等。在測量之前,先將反光片貼在保持架的軸向端面上,再將激光速度傳感器固定在磁鐵支架,并與數據采集系統和計算機相連,調整傳感器的位置,確保發出的光源能夠直射到反光條上。在測量過程中,保持架每旋轉一周,激光轉速傳感器即可接收一次反射的紅外信號,進而產生一個電壓脈沖,最終經過計算機進行相關數據處理便可得到保持架實際轉速。

為了檢驗該裝置測量轉速的精度,先以已知轉速的機械軸為測量對象預先做了一組實驗。將機械軸由停機狀態下勻加速到12 000 r/min,通過機械軸的實際轉速和測量轉速的對比,發現誤差小于1%,說明該裝置測量精度滿足驗證要求。

3.2 實驗結果對比分析

基于上述二次開發方案,對7008C角接觸球軸承動態特性開展仿真與實驗對比驗證,分析了軸向預緊力為50、100、150、200 N,內圈驅動轉速在2 000、4 000、6 000 r/min工況下的保持架轉速,結果如圖12所示。

(a)轉速2 000 r/min

(b)轉速4 000 r/min

(c)轉速6 000 r/min圖12 保持架實際轉速測量與仿真對比Fig.12 A comparison of actual measurement and simulation of cage speed

觀察圖12可以發現:本文建立的角接觸球軸承動力學模型在2 000 r/min時,預緊力對保持架轉速影響較弱,且保持架轉速的仿真值與實驗值較為接近,最大相差0.236 3 rad/s,此時相對誤差為0.26%。在4 000 r/min下,預緊力對保持架轉速的影響微小,隨著預緊力的增加,保持架轉速逐漸增加,而保持架轉速的仿真值與實驗值極為接近,此時最大相差0.227 0 rad/s,相對誤差為0.125%。在6 000 r/min下,預緊力對保持架的轉速影響較小,而保持架轉速的仿真值與實驗值最大相差0.630 3 rad/s,此時相對誤差為0.23%。實驗測量值與模型仿真值的相差不超過0.5%,且誤差產生的原因可能是模型中些許參數存在相關誤差,或保持架在實際測量過程中存在測量誤差,但整體上證明建立的軸承模型運動分析的準確性。

4 結 論

針對ADAMS自身接觸器計算精度不足的問題,本文基于ADAMS軟件自身動力學分析優勢,對其接觸求解器進行了相關開發,建立了角接觸球軸承動力學模型,對軸承之間的載荷分布與接觸進行了精確建模,并利用實驗驗證了模型的仿真計算精度,具體結論如下:

(1)基于軸承全局坐標系,利用Gfosub子程序編寫軸承零件之間的作用力子程序,根據求解方法的特點應用GSTIFF I3進行求解,實現了面向ADAMS軟件二次開發的接觸與運動精確建模與計算,經過與經驗公式對比,接觸力仿真誤差不超過1%;

(2)通過搭建保持架轉速測量實驗臺,對建立的動力學模型進行相關驗證,發現在低速時,模型能夠準確地仿真計算出保持架的轉速,且仿真與實驗測量誤差不超過0.5%,證明了模型的可靠性;

(3)在高速時,建立的模型仿真得到的保持架轉速與實驗測得轉速存在著一定的差別,還有待于進一步細化研究。