李 強(qiáng),許偉偉,王振波,金有海,鄭水英

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

滑動(dòng)軸承動(dòng)力特性的數(shù)值計(jì)算方法

李 強(qiáng)1,許偉偉2,王振波1,金有海1,鄭水英3

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

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

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

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

1 滑動(dòng)軸承CFD計(jì)算模型

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

1.1 控制方程

滑動(dòng)軸承中潤(rùn)滑流場(chǎng)可以用質(zhì)量、動(dòng)量守恒方程和空穴邊界方程來描述。

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

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

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

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

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

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

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

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

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

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

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

2 滑動(dòng)軸承動(dòng)力特性的計(jì)算

2.1 動(dòng)網(wǎng)格方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 計(jì)算初值的影響

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

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

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

3.2 計(jì)算結(jié)果

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

3.2.1 剛度系數(shù)

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

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

3.2.2 阻尼系數(shù)

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

渦動(dòng)軌跡如圖7所示。

圖7 軸頸圓形渦動(dòng)軌跡Fig.7 Circle whirling orbit of journal

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

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

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

3.2.3 結(jié)果對(duì)比

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

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

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

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

4 結(jié) 論

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

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

(3)CFD瞬態(tài)計(jì)算位移擾動(dòng)下的潤(rùn)滑流場(chǎng)時(shí),計(jì)算初值擾動(dòng)比較大,會(huì)影響瞬態(tài)計(jì)算結(jié)果的準(zhǔn)確性。

[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ā)展計(jì)劃項(xiàng)目(2009AA04Z413);國家自然科學(xué)基金項(xiàng)目(51275452)

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

李強(qiáng),許偉偉,王振波,等.滑動(dòng)軸承動(dòng)力特性的數(shù)值計(jì)算方法[J].中國石油大學(xué)學(xué)報(bào):自然科學(xué)版, 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.