王攀 劉保國(guó) 馮偉 趙耿

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

引言

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

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

1 油膜剛度阻尼計(jì)算方法

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

1.1 6DOF計(jì)算模型

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

(1)

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

(2)

(3)

R代表變換矩陣,

(4)

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

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

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

(1)基于彈簧的光順,任何兩個(gè)網(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)動(dòng)由擴(kuò)散方程控制：

(7)

1.3 剛度阻尼差分計(jì)算模型

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

(8)

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

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

(9)

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

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

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

圖1 深淺腔動(dòng)靜壓軸承結(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)格劃分方式對(duì)油膜模型進(jìn)行網(wǎng)格劃分.考慮油膜厚度方向油膜剪切應(yīng)力梯度大,沿厚度方向的網(wǎng)格等分為10層,網(wǎng)格數(shù)為100萬(wàn)個(gè),劃分后的網(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)行油膜流場(chǎng)數(shù)值仿真,邊界條件取進(jìn)油口為Pressure-inlet,供油壓力為Ps；出油口為Pressure-outlet,出油口是軸承與軸頸兩表面之間的微小間隙,出油口壓力等于外界大氣壓力,其值設(shè)為0；油膜內(nèi)壁面設(shè)置為旋轉(zhuǎn)壁面,轉(zhuǎn)速為500rad/s,其它面設(shè)置為靜止壁面,油膜計(jì)算參數(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 油膜計(jì)算參數(shù)Table 2 Oil film Calculation parameters

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

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

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

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

2.4 動(dòng)靜壓油膜軸承動(dòng)力特性計(jì)算流程

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

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

軸頸x方向和y方向位移隨時(shí)間的變化如圖8所示,軸頸在初始階段受到較大幅度的震蕩,隨著不斷的迭代計(jì)算軸頸的偏移不斷減小,在計(jì)算時(shí)間t=0.05s左右,軸頸趨于穩(wěn)定并在靜平衡位置附近做微小幅度的渦動(dòng),其幅值在千分之一微米級(jí),說明此系統(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方向位移隨時(shí)間t的變化Fig.8 The change of x direction and y direction displacement with time t

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

圖9 x方向和y方向受力隨時(shí)間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)速下動(dòng)靜壓軸承剛度阻尼系數(shù)

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

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

3 結(jié)論

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

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

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

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

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