滑動軸承動力特性的數(shù)值計算方法

李 強,許偉偉,王振波,金有海,鄭水英

(1.中國石油大學化學工程學院,山東青島 266580;2.中國石油大學儲運與建筑工程學院,山東青島 266580; 3.浙江大學化工機械研究所,浙江杭州 310027)

滑動軸承動力特性的數(shù)值計算方法

李 強1,許偉偉2,王振波1,金有海1,鄭水英3

(1.中國石油大學化學工程學院,山東青島 266580;2.中國石油大學儲運與建筑工程學院,山東青島 266580; 3.浙江大學化工機械研究所,浙江杭州 310027)

將計算流體動力學與簡諧激勵法應用于滑動軸承動力特性系數(shù)的求解,通過采用全新的變流域動網(wǎng)格技術(shù)提出一種基于瞬態(tài)流場計算的滑動軸承動特性的計算方法。利用所提出的方法計算典型滑動軸承的剛度、阻尼系數(shù),并與已有的經(jīng)典計算結(jié)果進行比較。結(jié)果表明:兩種方法的計算結(jié)果基本吻合,提出的新方法不僅考慮了計算初值的影響,而且計入了實際存在的油膜破裂現(xiàn)象,滿足精度要求,有效可行。

滑動軸承;計算流體力學(CFD);動特性系數(shù);動網(wǎng)格;數(shù)值計算

滑動軸承廣泛地應用于大型、高速旋轉(zhuǎn)機械,其油膜力既是轉(zhuǎn)子－軸承系統(tǒng)阻尼的主要來源,也是導致機組穩(wěn)定性下降的重要原因,滑動軸承動態(tài)特性對轉(zhuǎn)子振動有重要影響。能精確描述滑動軸承動力性能的數(shù)學模型是非線性的,在非線性模型還不完善的前提下,在線性假定的基礎上采用4個剛度系數(shù)和4個阻尼系數(shù)表示小擾動下油膜的動態(tài)特性比較準確,可以滿足通常的工程實際需要。為此國內(nèi)外學者從理論和實驗兩方面做了大量工作[1-5]。通過計算流體力學方法(CFD)直接求解N-S方程來分析滑動軸承的3D潤滑流場,該方法便于求解復雜區(qū)域上的問題,對于事先未知的自由邊界或求解區(qū)域內(nèi)部不同介質(zhì)的交界面比較容易處理,因此基于CFD技術(shù)的3D流場計算方法廣泛地應用于滑動軸承性能的研究[6-9]。考慮到現(xiàn)有動網(wǎng)格方法在計算滑動軸承瞬態(tài)流場時網(wǎng)格會出現(xiàn)較大的畸變,筆者利用CFD-FLUENT軟件對液體動壓滑動軸承進行建模,采用全新的變流域動網(wǎng)格技術(shù)對滑動軸承內(nèi)瞬態(tài)流場進行計算分析,提出一種滑動軸承動特性系數(shù)的計算方法,并將計算結(jié)果和已有的數(shù)值計算結(jié)果進行對比,驗證提出方法的有效性和可行性。

1 滑動軸承CFD計算模型

實際工作時滑動軸承在軸承開擴區(qū)有油和氣混合存在,因此為了更準確地描述滑動軸承動力特性,本文建立的CFD計算模型可以考慮兩相流對滑動軸承性能的影響。

1.1 控制方程

滑動軸承中潤滑流場可以用質(zhì)量、動量守恒方程和空穴邊界方程來描述。

在邊界運動的任意控制體V內(nèi)積分形式的非定常不可壓黏性流的連續(xù)性方程和動量方程的通式為

式中,ρm為氣液兩相的混合密度,kg/m3;φ為通量; dV為控制體V的邊界;vm為氣液兩相的速度矢量;vs為動網(wǎng)格的網(wǎng)格變形速度矢量;Γ為擴散系數(shù);Sφ為通量φ的源項。

為了充分地考慮實際運行中滑動軸承內(nèi)潤滑油的分離和重新合并的情況,本文中選用的空穴模型[10]是基于質(zhì)量守恒邊界條件的,

式中,vv為氣相速度矢量;f為氣相質(zhì)量分數(shù);γ為有效交換系數(shù);Re和Rc分別為空穴的生成和凝聚率。

1.2 建模及網(wǎng)格劃分

以圓柱滑動軸承為研究對象,圖1給出了滑動軸承的3D網(wǎng)格結(jié)構(gòu)。軸承水平中分面兩側(cè)各有1個進油槽,潤滑油由兩側(cè)進油口進入軸承間隙內(nèi),由軸向兩端流出。參照文獻[11],軸承直徑D=50 mm,軸承寬度B=25 mm,軸頸半徑間隙c=0.05 mm,油槽包角α=30°,油槽軸向長度b=3.75 mm,軸頸偏心率ε=0.5,轉(zhuǎn)速Ω=1 000 rad/s,潤滑油的動力黏度為μ=1.25×10－2N·s/m2,氣態(tài)潤滑油參數(shù)取空氣參數(shù),進油壓力p=103 kPa,流體流動狀態(tài)為層流。

利用前處理軟件Gambit對滑動軸承進行網(wǎng)格劃分,由于滑動軸承潤滑流場在空間3個方向的尺寸大小相差較大,須采用結(jié)構(gòu)化網(wǎng)格來提高數(shù)值計算精度和節(jié)省計算時間。對不同網(wǎng)格密度下CFD計算模型的計算結(jié)果進行比較分析后發(fā)現(xiàn),油膜間隙處不同的網(wǎng)格數(shù)導致計算結(jié)果差別較大;而周向和軸向網(wǎng)格密度對計算結(jié)果影響不大,因此這兩個方向的網(wǎng)格尺寸主要從長寬比角度考慮(最佳長寬比為1∶1)。選用軸承間隙徑向5層網(wǎng)格,軸向和周向網(wǎng)格密度為0.2。

圖1 滑動軸承網(wǎng)格模型Fig.1 Grid structure of journal bearing model

1.3 FLUENT計算參數(shù)設置

對于滑動軸承的潤滑流場,選擇3D耦合求解器,隱式格式求解;流場計算的邊界條件要用到進口、出口和壁面3種邊界條件。滑動軸承的進口為壓力邊界條件,給定進口總壓力;滑動軸承的出口同樣為壓力邊界條件,壓力根據(jù)實際情況設定為環(huán)境大氣壓;軸瓦表面為固定無滑移邊界,近壁面應用標準壁面函數(shù),軸頸表面繞偏心位置旋轉(zhuǎn)。

采用有限體積法離散控制方程、連續(xù)方程,動量方程采用一階迎風格式,壓力差值格式采用PRESTO格式,壓力速度耦合采用SIMPLEC算法。計算時所有方程的殘差都小于1×10－4,并在計算過程中對進、出口流體質(zhì)量流量進行監(jiān)控,數(shù)值基本相等時計算收斂。

2 滑動軸承動力特性的計算

2.1 動網(wǎng)格方法

由于軸承徑向間隙與其他方向的結(jié)構(gòu)尺寸相差比較大,利用FLUENT軟件自帶的動網(wǎng)格方法進行網(wǎng)格更新將產(chǎn)生較大的網(wǎng)格畸變,影響瞬態(tài)計算的準確性,因此在進行伴隨軸頸移動的網(wǎng)格更新時采用自行開發(fā)的變流域動網(wǎng)格技術(shù)。首先將軸承間隙流場用結(jié)構(gòu)化網(wǎng)格進行劃分,當軸頸移動時,軸承間隙流場最內(nèi)圈上的所有節(jié)點都隨軸頸一起移動,最外圈和油槽上的節(jié)點保持不動,中間的節(jié)點按照結(jié)構(gòu)化網(wǎng)格劃分原理進行移動,因此可以通過編程計算出每個節(jié)點坐標的變化值,然后通過FLUENT的UDF接口將各個節(jié)點移至新的位置。圖2給出了軸頸移動到一定位置后軸承網(wǎng)格變化情況(圖中油膜間隙被大大地夸大)。從圖2可以看出,通過該方法移動后的網(wǎng)格即使在很大的軸頸偏心下依然能保證良好的光滑性和規(guī)整性,而且基本不會出現(xiàn)負體積。

軸頸中心的位置隨著軸頸的移動不斷變化,因此軸頸表面的旋轉(zhuǎn)速度不能按常規(guī)的繞固定點旋轉(zhuǎn)的方法設定,本文中通過FLUENT－UDF程序捕捉到每一步的軸頸中心,然后將旋轉(zhuǎn)速度分解為x、y兩個方向的速度施加到軸頸表面。

圖2 軸頸移動后網(wǎng)格節(jié)點截面圖Fig.2 Front view of grids in moved journal bearing

2.2 剛度阻尼系數(shù)計算方法

軸頸在靜平衡位置受到外載荷作用時,對軸頸油膜產(chǎn)生位移和速度擾動,油膜力將發(fā)生變化,力的變化和擾動之間的關(guān)系一般是非線性的,為簡化分析,當擾動較小時,常將這種關(guān)系線性化,如圖3所示。油膜力可近似由下列線性關(guān)系表達:

式中,Fx、Fy為油膜力在x、y方向的分量;Fx0、Fy0為靜態(tài)平衡位置時,油膜力在x、y方向的分量;kij、cij分別為軸頸中心在靜平衡位置(x0,y0)時油膜的剛度(N/m)和阻尼系數(shù),N·s/m。

由式(3)可以得到油膜力的增量,即動態(tài)油膜力,寫成矩陣的形式為

式中,[K]為剛度系數(shù)矩陣;[C]為阻尼系數(shù)矩陣。

圖3 動力特性系數(shù)示意圖Fig.3 Schematic diagram of dynamic coefficients

動力特性系數(shù)決定了小擾動下轉(zhuǎn)子－軸承系統(tǒng)的動特性。提出了一種基于CFD動網(wǎng)格技術(shù)的動特性系數(shù)計算新方法,首先利用軸頸在x、y兩個方向的平移求出剛度系數(shù),然后在此基礎上利用簡諧激振原理計算求得阻尼系數(shù)。

圖4 動力特性系數(shù)求解流程圖Fig.4 Schematic diagram of solution of dynamic coefficients

圖4給出了利用新方法求解剛度和阻尼系數(shù)的過程。首先利用位移增量原理求出滑動軸承的剛度系數(shù),使軸頸從平衡位置水平位移Δx(這時Δy= 0),計算時須在軸頸上加力ΔFx和ΔFy。這時可根據(jù)定義求得剛度系數(shù)kxx=ΔFx/Δx,kyx=ΔFy/Δx。同理通過使軸頸從平衡位置豎直位移Δy(這時Δx =0),可得到kxy=ΔFx/Δy,kyy=ΔFy/Δy。

同樣以滑動軸承的靜平衡位置為基準,分別在x、y方向施加頻率為ω,相位差90°的已知簡諧位移:

通過自編FLUENT-UDF程序積分求解得到這時作用在軸頸上的油膜力:

假定轉(zhuǎn)子的質(zhì)量忽略不計,則將式(5)、(6)代入牛頓第二定律可以很容易得到阻尼系數(shù):

除了利用軸頸平移先求得剛度系數(shù)外,還可以對軸頸分別施加兩次線性獨立的簡諧位移激勵,從而直接求解出所需的4個剛度系數(shù)和阻尼系數(shù)。

3 數(shù)值計算結(jié)果

3.1 計算初值的影響

鑒于CFD求解方法能更好地考慮軸承內(nèi)部3D瞬態(tài)流場的影響和表征流體流動形態(tài),目前利用CFD技術(shù)求解滑動軸承動特性的方法主要有兩種:一是利用相對坐標系的方法將瞬態(tài)流場近似為穩(wěn)態(tài)流場求解,該方法僅局限于渦動軌跡為圓的圓柱軸承,若渦動軌跡或軸承結(jié)構(gòu)稍作變化則無法使用;第二種較常用的方法是利用小擾動法直接瞬態(tài)求解,但該方法沒有考慮計算初始值的影響,結(jié)果的準確性有待商榷。

非穩(wěn)態(tài)計算屬于微分方程的初值計算問題,不同的初值在計算初期得到的計算結(jié)果并不一樣。圖5給出了軸頸做指定的直線運動時油膜力的變化情況。從圖5可以看出,數(shù)值計算初期,結(jié)果振蕩比較大,隨著計算的進行,數(shù)值結(jié)果逐漸趨于穩(wěn)定。因此,傳統(tǒng)的不考慮初值影響而只根據(jù)軸頸移動時的計算初始值求解得到的動特性系數(shù)是不準確的。

圖5 位移擾動下油膜力的變化過程Fig.5 Oil film versus displacement step under displacement perturbation

3.2 計算結(jié)果

為了驗證上述計算滑動軸承動力特性系數(shù)方法的適用性和準確性,對文獻[11]中提供的典型滑動軸承模型進行計算,并與該文獻提供的計算結(jié)果進行對比。其中VT-FAST,DyRoBes-BePerf和VT-EXPRESS都是目前應用比較廣泛的基于經(jīng)典Reynolds方程求解的轉(zhuǎn)子動力學計算軟件。

3.2.1 剛度系數(shù)

為了保證剛度系數(shù)求解的準確性,軸頸的平移距離要足夠小,分別取位移步長(0.1、0.2 μm)使軸頸從靜平衡位置沿x正方向和y負方向移動,得到的油膜力和擬合曲線如圖6所示。從圖6可知,隨著軸頸在x和y方向的移動,所需施加在軸頸上的力與位移呈明顯的線性比例關(guān)系,該比值即為剛度系數(shù)(N/m):

圖6 剛度系數(shù)的擬合求解Fig.6 Fitting of stiffness coefficients

3.2.2 阻尼系數(shù)

不同于剛度系數(shù)的求解,阻尼系數(shù)是與軸頸擾動速度有關(guān)的量,因此本文中利用簡諧激勵原理進行求解。首先令軸頸以下列形式繞靜平衡位置做圓軌跡渦動:

渦動軌跡如圖7所示。

圖7 軸頸圓形渦動軌跡Fig.7 Circle whirling orbit of journal

然后根據(jù)計算值擬合轉(zhuǎn)子表面的油膜力,可得擬合公式為

將上述擬合值代入式(7)即可得到4個阻尼系數(shù)為

油膜力及其擬合值如圖8所示。為了消除初始值的影響,通常需要將初期的計算結(jié)果舍棄。

3.2.3 結(jié)果對比

表1、2中的第4組數(shù)據(jù)分別為根據(jù)圖6、8擬合計算出的剛度系數(shù)和阻尼系數(shù)。從表中可以看出,針對簡單滑動軸承結(jié)構(gòu),本文計算的剛度和阻尼值與傳統(tǒng)計算方法得到的值基本一致,雖然交叉系數(shù)的符號不一樣,但經(jīng)對比發(fā)現(xiàn),符號與文獻[12]中的經(jīng)典數(shù)據(jù)完全一致,因此本文提出的基于CFD技術(shù)的滑動軸承動特性系數(shù)求解方法是可行的,且求解精度滿足要求。該方法的優(yōu)點是針對多油楔、帶溝槽等復雜結(jié)構(gòu)滑動軸承的動特性系數(shù)求解不需要重新推導復雜公式,并且可以很容易地考慮油膜本身質(zhì)量的影響,尤其針對高速輕型轉(zhuǎn)子軸承系統(tǒng),慣性項的影響是至關(guān)重要的。

圖8 x、y方向的油膜力變化及其擬合值Fig.8 Oil film force in x-axis and y-axis direction during circle whirling

表1 圓柱滑動軸承剛度系數(shù)Table 1 Stiffness coefficients of journal bearing

表2 圓柱滑動軸承阻尼系數(shù)Table 2 Damping coefficients of column journal bearing

4 結(jié) 論

(1)提出的新的變流域動網(wǎng)格方法理論上適合大部分結(jié)構(gòu)的滑動軸承,實時地反映軸承潤滑流場的動態(tài)變化。

(2)新方法的計算結(jié)果與已有的經(jīng)典數(shù)據(jù)基本吻合,驗證了該計算方法的可行性和有效性。

(3)CFD瞬態(tài)計算位移擾動下的潤滑流場時,計算初值擾動比較大,會影響瞬態(tài)計算結(jié)果的準確性。

[1] TIEU A K,QIU Z L.Identification of sixteen dynamic coefficients of two journal bearings from experimental unbalance response[J].Wear,1994,177(1):63-69.

[2] ZHENG T,HASEBE N.Calculation of equilibrium position and dynamic coefficients of a journal bearing using free boundary theory[J].ASME Journal of Tribology, 2000,122:616-621.

[3] PERRINATO B,FLACK R D.Test results for a highly preloaded three-lobe journal bearing-effect of load orientation on static and dynamic characteristics[J].Journal of the Society of Tribology and Lubrication Engineers,2001, 57(9):23-30.

[4] WANG J K,KHONSARI M M.A new derivation for journal bearing stiffness and damping coefficients in polar coordinates[J].Journal of Sound and Vibration,2006, 290:500-507.

[5] BRAD A M,SAMUEL A H.Identifying bearing rotor-dynamic coefficients using an extended Kalman filter[J].Tribology Transactions,2009,52:671-679.

[6] GANDJALIKHAN N S A.Inertia effect on the thermohydrodynamic characteristics of journal bearings[J].Journal of Engineering Manufacture,2005,219(6):459-467.

[7] SHEIKH N Y,GANDJALIKHAN N S A.Three-dimensional numerical analysis of hydrodynamic characteristics of axial groove journal bearings running with ferrofluids under magnetic field[J].Journal of Engineering Tribology,2010,224(7):609-619.

[8] PANDAY K M,CHOUDHURY P L,KUMAR N P.Numerical unsteady analysis of thin film lubricated journal bearing[J].International Journal of Engineering and Technology,2012,4(2):185-191.

[9] DELIGANT M,PODEVIN P,DESCOMBES G.CFD model for turbocharger journal bearing performances[J].Applied Thermal Engineering,2011,31:811-819.

[10] SINGHAL A K,ATHAVALE M M,LI H Y,et al.Mathematical basis and validation of the full cavitation model[J].Journal of Fluids Engineering,2002,124 (3):617-624.

[11] GUO Z,HIRANO T,KIRK R G.Application of CFD analysis for rotating machinery:part 1-hydrodynamic,hydrostatic bearings and squeeze film damper[J].ASME Journal of Engineering for Gas Turbines and Tower, 2005,127:445-451.

[12] LUND J W.Rotor-bearing dynamics[R].Lecture Notes:Technical University of Denmark,1979.

(編輯 沈玉英)

Numerical calculation method of dynamic characteristics of journal bearing

LI Qiang1,XU Weiwei2,WANG Zhenbo1,JIN Youhai1,ZHENG Shuiying3
(1.College of Chemical Engineering in China University of Petroleum,Qingdao 266580,China; 2.College of Transport&Storage and Civil Engineering in China University of Petroleum,Qingdao 266580,China; 3.Institute of Chemical Machinery,Zhejiang University,Hangzhou 310027,China)

The computational fluid dynamics(CFD)and harmonic excitation method were applied to the numerical calculation of dynamic characteristics of journal bearing.By employing a new mesh movement approach based on structured grid,a new approach for calculating the dynamic characteristics of journal bearing was proposed based on the transient flow calculation.The stiffness and damping coefficients of a typical bearing were calculated by applying the new approach.The results obtained from the method were compared with previous classic computation results.The results show that the computation results of two methods are consistent.The effects of the computational initial value and the oil film fracture phenomenon are considered in this method,which is suitable for most of the journal bearing structures.The numerical method has good accuracy,and the method is valid.

journal bearing;computational fluid dynamic(CFD);dynamic coefficient;dynamic mesh;numerical calculation

TH 133

A

1673-5005(2014)05-0165-07

10.3969/j.issn.1673-5005.2014.05.024

2013－12－22

國家“863”高技術(shù)研究發(fā)展計劃項目(2009AA04Z413);國家自然科學基金項目(51275452)

李強(1984－),男,講師,博士,從事轉(zhuǎn)子動力學、振動測試和故障診斷研究。E－mail:liqiangsydx@163.com。

李強,許偉偉,王振波,等.滑動軸承動力特性的數(shù)值計算方法[J].中國石油大學學報:自然科學版, 2014,38(5):165-171.

LI Qiang,XU Weiwei,WANG Zhenbo,et al.Numerical calculation method of dynamic characteristics of journal bearing [J].Journal of China University of Petroleum(Edition of Natural Science),2014,38(5):165-171.