基于6DOF模型及動網(wǎng)格的動靜壓軸承剛度阻尼數(shù)值計算*

王攀 劉保國 馮偉 趙耿

(河南工業(yè)大學(xué) 機(jī)電工程學(xué)院, 鄭州 450001)

引言

液體動靜壓軸承具有回轉(zhuǎn)精度高、動態(tài)剛性和阻尼減振性能好、使用壽命長等優(yōu)點(diǎn),在超高速精密磨削電主軸領(lǐng)域有著廣泛的應(yīng)用前景[1].電主軸系統(tǒng)的高速穩(wěn)定性與軸承動力學(xué)特性緊密相關(guān),因此對于軸承動態(tài)剛度和阻尼的精確計算成為研究的重點(diǎn).

目前油膜軸承的數(shù)值計算模型主要有兩種,一種是運(yùn)用流體潤滑理論并通過條件假設(shè)建立經(jīng)典的雷諾方程,如Rowe等[2]基于經(jīng)典雷諾方程的有限差分法,用有限擾動法得到了非線性化的動靜壓軸承動力特性系數(shù).賀玉嶺等[3]采用類似的方法,建立了油膜在靜平衡位置做微小擾動的計算模型,進(jìn)而得到動靜壓軸承的剛度和阻尼系數(shù). 然而,該方法不能很好描述擴(kuò)散效應(yīng)、擠壓效應(yīng)、慣性效應(yīng)、粘溫效應(yīng)等[4-7]對油膜軸承所造成的影響,計算誤差較大.另一種是通過Navier-Stokes方程建立軸承內(nèi)部三維流場模型及邊界條件,并以CFD軟件為求解器進(jìn)行求解.該方法在涉及復(fù)雜的流動幾何形狀時更為有效,因此近年來得到了更多的應(yīng)用.Guo等[8]用CFD程序?qū)σ后w動靜壓軸承的靜態(tài)和動態(tài)特性進(jìn)行了計算,通過與標(biāo)準(zhǔn)潤滑理論數(shù)值計算結(jié)果對比證明了該方法的有效性.為了反映油膜網(wǎng)格在外界干擾下的扭曲和變形,熊萬里[9]等提出基于動網(wǎng)格模型的剛度阻尼計算方法,使油膜力計算結(jié)果與高速時的工況更為吻合.采用動網(wǎng)格法求解軸承剛度和阻尼時,首先需要準(zhǔn)確計算軸頸在外載荷作用時的靜平衡位置,以便對平衡位置施加位移擾動和速度擾動.本文基于6DOF模型計算軸頸平衡位置,采用動網(wǎng)格更新方法實(shí)現(xiàn)軸頸在該位置的擾動,通過求解Navier-Stokes方程得到軸承剛度和阻尼.

1 油膜剛度阻尼計算方法

采用6DOF模型及動網(wǎng)格的計算方法是在CFD軟件FLUENT基礎(chǔ)上,通過加載6DOF自定義程序,計算外載荷作用下軸頸的非線性軸心軌跡,從而得到軸頸在外載荷作用下的靜平衡位置.通過嵌入UDF程序DEFINE-CG-MOTION以動網(wǎng)格更新方法來實(shí)現(xiàn)對處于靜平衡位置的軸頸施加擾動；為有效抑制網(wǎng)格在施加擾動過程中的扭曲變形,通過嵌入UDF程序DEFINE-PROFILE來定義軸頸(油膜內(nèi)壁)為旋轉(zhuǎn)運(yùn)動.最后通過求解Navier-Stokes方程得到軸頸擾動前后位置變化的瞬態(tài)油膜力,結(jié)合數(shù)值計算的差分方法,得到動靜壓軸承油膜剛度、阻尼動力特性系數(shù).

1.1 6DOF計算模型

ANSYS FLUENT中的6 DOF模型,可以利用作用于物體的力和力矩來計算物體重心的平移和旋轉(zhuǎn)運(yùn)動.其在慣性坐標(biāo)系下求解重心平移運(yùn)動的控制方程為：

(1)

物體的旋轉(zhuǎn)運(yùn)動在體坐標(biāo)系下的控制方程為：

(2)

(3)

R代表變換矩陣,

(4)

式中,Cλ=cos(λ),Sλ=sin(λ),角φ,θ,ψ分別為繞x,y,z軸旋轉(zhuǎn)的歐拉角.平移速度和旋轉(zhuǎn)角速度可由式(1)和(2)分別迭代計算得出,在6 DOF模型中網(wǎng)格位置的更新由平移速度和旋轉(zhuǎn)角速度實(shí)現(xiàn).

1.2 加載UDF宏的動網(wǎng)格更新方法

動網(wǎng)格方法可用來模擬流體域的形狀由于邊界移動而隨時間產(chǎn)生的變化,油膜網(wǎng)格運(yùn)動需要用光順模型來調(diào)整變形邊界區(qū)域的網(wǎng)格,本文選用基于彈簧的光順模型和基于擴(kuò)散的光順模型.

(1)基于彈簧的光順,任何兩個網(wǎng)格節(jié)點(diǎn)之間的邊界可以理想化為相互連接的彈簧,且邊界節(jié)點(diǎn)處的位移將產(chǎn)生成比例的彈簧力,網(wǎng)格節(jié)點(diǎn)上的作用力可寫為：

(5)

式中,Δxi和Δxj為相鄰節(jié)點(diǎn)之間的位移,ni為相鄰節(jié)點(diǎn)的數(shù)目,kij為相鄰節(jié)點(diǎn)之間的彈簧剛度系數(shù),其相鄰節(jié)點(diǎn)彈簧剛度系數(shù)的定義為:

(6)

式中,kfac為彈簧常數(shù)因子需要自己設(shè)定,取值范圍在0與1之間.

(2)基于擴(kuò)散光順,網(wǎng)格運(yùn)動由擴(kuò)散方程控制：

(7)

1.3 剛度阻尼差分計算模型

當(dāng)軸頸受到擾動在靜平衡位置偏移一微小距離時,油膜反向承載力的增量對該微小擾動距離的比值稱為剛度.在分析計算時,常常把微小偏移分解為沿坐標(biāo)軸x和y方向的分量,油膜反力的增量也分解成沿x和y方向的分量,分別對其進(jìn)行差分計算,得四個剛度系數(shù)：

(8)

式中,ΔFdij為擾動位移引起的油膜力變化,Δx、Δy為擾動位移.用i,j表示x,y中的某一個.

當(dāng)軸頸在靜平衡位置獲得一擾動速度時,油膜阻尼力的增量與該擾動速度的比值稱為阻尼.把油膜阻尼力增量和軸頸擾動速度的增量分別分解到x軸和y軸上,可以得到四個阻尼系數(shù)：

(9)

2 油膜軸承剛度阻尼數(shù)值計算

2.1 油膜模型建立及網(wǎng)格劃分

以深淺腔液體動靜壓軸承為研究對象,軸承結(jié)構(gòu)如圖1所示,結(jié)構(gòu)參數(shù)如表1所示.

圖1 深淺腔動靜壓軸承結(jié)構(gòu)圖Fig.1 Structure diagram of Deep/shallow pocket hybrid bearing

表1 軸承基本結(jié)構(gòu)參數(shù)Table 1 Basic structural parameters of the bearing

在Solidworks里建立油膜實(shí)體模型,如圖2所示.將模型導(dǎo)入ANSYS ICEM CFD中,采用O型Block分塊的結(jié)構(gòu)化網(wǎng)格劃分方式對油膜模型進(jìn)行網(wǎng)格劃分.考慮油膜厚度方向油膜剪切應(yīng)力梯度大,沿厚度方向的網(wǎng)格等分為10層,網(wǎng)格數(shù)為100萬個,劃分后的網(wǎng)格模型如圖3所示,進(jìn)油孔部位及油膜厚度部位的局部放大如圖4和圖5所示.

圖2 三維油膜實(shí)體Fig.2 Three-dimensional oil film entity

2.2 邊界條件的確定

將網(wǎng)格文件導(dǎo)入FLUENT進(jìn)行油膜流場數(shù)值仿真,邊界條件取進(jìn)油口為Pressure-inlet,供油壓力為Ps；出油口為Pressure-outlet,出油口是軸承與軸頸兩表面之間的微小間隙,出油口壓力等于外界大氣壓力,其值設(shè)為0；油膜內(nèi)壁面設(shè)置為旋轉(zhuǎn)壁面,轉(zhuǎn)速為500rad/s,其它面設(shè)置為靜止壁面,油膜計算參數(shù)表2所示.

圖3 油膜網(wǎng)格劃分模型Fig.3 Model of oil film mesh

圖4 油孔過渡區(qū)O型劃分Fig.4 O-Shape division of oil-hole transition zone

圖5 油膜厚度方向局部放大Fig.5 Partial amplification of oil film thickness direction

表2 油膜計算參數(shù)Table 2 Oil film Calculation parameters

2.3 動網(wǎng)格UDF宏程序

動靜壓軸承油膜內(nèi)壁在主軸運(yùn)轉(zhuǎn)過程中受到外加載荷的作用而產(chǎn)生偏移,此過程屬于被動型動網(wǎng)格計算問題,可以編寫6DOF自定義程序來解決.首先需要明確模型中6個方向的自由度,動靜壓軸承油膜在實(shí)際運(yùn)動過程中,需要限制其在軸向的竄動,以及繞三個坐標(biāo)軸方向的轉(zhuǎn)動.定義軸承所受外加載荷為Fx=0,Fy≠0,轉(zhuǎn)軸質(zhì)量為m,軸頸轉(zhuǎn)向?yàn)槟鏁r針,描述油膜內(nèi)壁受外加載荷的運(yùn)動示意圖如圖6所示.油膜內(nèi)壁在外加載荷Fy的作用下做剛體運(yùn)動,油膜內(nèi)壁受力后向右下角移動.

圖6 油膜內(nèi)壁在外載荷作用下的運(yùn)動示意圖Fig.6 Motion diagram of oil film inner wall under external load

通過編寫UDF程序DEFINE-CG-MOTION來指定特定區(qū)域的網(wǎng)格運(yùn)動,其方法是在每個時間步長內(nèi)為動態(tài)網(wǎng)格區(qū)域提供線速度和角速度,FLUENT利用這些速度更新動態(tài)區(qū)域上的節(jié)點(diǎn)位置.

2.4 動靜壓油膜軸承動力特性計算流程

在6DOF程序里定義油膜內(nèi)壁受力為Fx=0N,Fy=1000N.將6DOF程序加載到FLUENT求解器并作用于油膜內(nèi)壁.在外載荷作用下油膜內(nèi)壁網(wǎng)格產(chǎn)生偏移并穩(wěn)定于一位置,此時油膜由偏心所產(chǎn)生的油膜反力為Fx=-0.529N,Fy=998.7N.軸頸在平衡位置一般應(yīng)滿足條件|Fx/Fy|<0.001,經(jīng)驗(yàn)算滿足條件.

外加載荷力作用下軸頸的軸心軌跡變化曲線如圖7所示,可以看出其軸心軌跡是收斂型的,軸心由初始位置原點(diǎn)出發(fā),經(jīng)過一定幅度的渦動,最終平衡穩(wěn)定于一點(diǎn),此時油膜反力與軸頸重力和外加載荷保持平衡.

軸頸x方向和y方向位移隨時間的變化如圖8所示,軸頸在初始階段受到較大幅度的震蕩,隨著不斷的迭代計算軸頸的偏移不斷減小,在計算時間t=0.05s左右,軸頸趨于穩(wěn)定并在靜平衡位置附近做微小幅度的渦動,其幅值在千分之一微米級,說明此系統(tǒng)處于穩(wěn)定狀態(tài).

圖7 外加載荷力作用下軸心軌跡變化曲線Fig.7 The change curve of axis trajectory under the action of external load force

圖8 x方向和y方向位移隨時間t的變化Fig.8 The change of x direction and y direction displacement with time t

軸頸x方向和y方向的油膜力隨時間的變化曲線如圖9所示,當(dāng)軸頸達(dá)到穩(wěn)態(tài)的時候,油膜反力抵消了外載荷Fy的作用,軸頸處于穩(wěn)態(tài)狀態(tài).

圖9 x方向和y方向受力隨時間t的變化Fig.9 The change of force in x direction and y direction with time t

表3 油膜剛度阻尼系數(shù)Table 3 Oil film stiffness and damping coefficient

2.5 不同轉(zhuǎn)速下動靜壓軸承剛度阻尼系數(shù)

采用上述動靜壓軸承剛度阻尼計算方法,對不同轉(zhuǎn)速的動靜壓軸承剛度阻尼進(jìn)行計算,得到不同轉(zhuǎn)速下動靜壓軸承油膜的4個剛度系數(shù)和4個阻尼系數(shù),如圖10和圖11所示.從圖中可以看出,隨著軸頸轉(zhuǎn)速的增加,軸承直接剛度系數(shù)Kxx、Kyy和交叉剛度系數(shù)Kxy、Kyx絕對值逐漸增大.直接阻尼系數(shù)Cxx、Cyy和交叉阻尼系數(shù)Cyx、Cxy絕對值也逐漸增大.

圖10 各轉(zhuǎn)速下油膜剛度系數(shù)Fig.10 The stiffness coefficient of oil film at each rotational speed

3 結(jié)論

本文采用6DOF模型及動網(wǎng)格的計算方法對具有典型結(jié)構(gòu)的液體動靜壓軸承的剛度阻尼系數(shù)進(jìn)行數(shù)值計算,通過研究分析得出以下結(jié)論:

1)采用6DOF模型可以反映油膜厚度在外載荷作用下的運(yùn)動變化過程,能夠更為準(zhǔn)確地計算軸頸靜平衡位置,使計算過程更符合工程實(shí)際；

圖11 各轉(zhuǎn)速下油膜阻尼系數(shù)Fig.11 The damping coefficient of oil film at each rotational speed

2)利用UDF宏程序以動網(wǎng)格方法實(shí)現(xiàn)對軸頸靜平衡位置的擾動,通過計算軸頸位移擾動、速度擾動前后的瞬態(tài)油膜力,利用差分法求得軸承剛度和阻尼系數(shù)；

3)分析了不同轉(zhuǎn)速下軸承剛度和阻尼的變化規(guī)律,結(jié)果表明隨著主軸轉(zhuǎn)速的增加軸承剛度系數(shù)和阻尼系數(shù)絕對值不斷增大.