吳 超 尹雪梅 李蒙蒙 李奕君 王 文
1.鄭州輕工業(yè)大學(xué)機(jī)電工程學(xué)院,鄭州,4500022.鄭州輕工業(yè)大學(xué)能源與動(dòng)力工程學(xué)院,鄭州,4500023.上海大學(xué)機(jī)電工程與自動(dòng)化學(xué)院,上海,200072
計(jì)算油膜軸承性能的關(guān)鍵問(wèn)題就是求解動(dòng)態(tài)流體潤(rùn)滑方程,得到油膜的壓力分布[1-4]。通常情況下,采用動(dòng)態(tài)雷諾方程計(jì)算軸承的性能時(shí),難以精確反映轉(zhuǎn)速所引起的油流周向慣性效應(yīng)、動(dòng)態(tài)擠壓效應(yīng)和靜壓效應(yīng)之間的線(xiàn)性耦合關(guān)系及其對(duì)油膜三維壓力場(chǎng)、溫度場(chǎng)和速度場(chǎng)的影響,因此,有必要直接求解Navier-Stokes方程,精確研究軸承參數(shù)對(duì)油膜性能的影響[5]。文獻(xiàn)[6-8]采用計(jì)算流體動(dòng)力學(xué)(computational fluid dynamics,CFD)開(kāi)展了靜壓軸承承載特性的研究,證實(shí)了在表征復(fù)雜求解域流體流動(dòng)形態(tài)方面Navier-Stokes方程可彌補(bǔ)雷諾方程的不足。文獻(xiàn)[9-10]利用Fluent軟件,采用靜網(wǎng)格方法,研究了氣穴現(xiàn)象對(duì)滑動(dòng)軸承性能的影響,但靜網(wǎng)格方法不能計(jì)算變載荷下軸心的軌跡。文獻(xiàn)[11-12]將軸頸的旋轉(zhuǎn)動(dòng)邊界轉(zhuǎn)換為靜邊界,并利用Smoothing動(dòng)網(wǎng)格模型,成功計(jì)算了油膜軸承的動(dòng)特性系數(shù)。文獻(xiàn)[13-15]自定義動(dòng)網(wǎng)格更新程序,提出了基于彈性變形的動(dòng)網(wǎng)格調(diào)整法,能用來(lái)求解瞬態(tài)軸心軌跡,并利用弱耦合算法研究了軸承剛度隨方向的變化以及渦動(dòng)中心與載荷、不平衡量的關(guān)系。文獻(xiàn)[16]針對(duì)求解復(fù)雜轉(zhuǎn)子-軸承系統(tǒng)非線(xiàn)性動(dòng)力學(xué)特性的問(wèn)題,基于Smoothing動(dòng)網(wǎng)格技術(shù),提出了一種計(jì)算流體力學(xué)和計(jì)算轉(zhuǎn)子動(dòng)力學(xué)的流固耦合新方法,計(jì)算結(jié)果表明,該方法能夠得到精確的軸心軌跡,并能準(zhǔn)確分析復(fù)雜轉(zhuǎn)子-軸承系統(tǒng)非線(xiàn)性動(dòng)力學(xué)特性。文獻(xiàn)[17]利用非定常動(dòng)網(wǎng)格技術(shù)建立了考慮軸頸渦動(dòng)頻率與渦動(dòng)軌跡的滑動(dòng)軸承動(dòng)力特性求解模型,研究了不同的軌跡下軸頸渦動(dòng)頻率和偏心率對(duì)滑動(dòng)軸承動(dòng)力特性的影響。以上在油膜軸承性能計(jì)算中運(yùn)用的動(dòng)網(wǎng)格方法是基于Smoothing模型或其變形模型基礎(chǔ)上進(jìn)行的,能夠方便地求解出油膜軸承的特性參數(shù),分析軸承-轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)性能。
Smoothing動(dòng)網(wǎng)格模型[18]的特點(diǎn)是不改變網(wǎng)格節(jié)點(diǎn)間的拓?fù)潢P(guān)系,只進(jìn)行網(wǎng)格形狀的改變,能夠計(jì)算變載荷下軸承的軸心軌跡,但Smoothing模型網(wǎng)格位移的更新方法容易造成網(wǎng)格畸變,當(dāng)偏心率過(guò)大時(shí),垂直于軸頸的網(wǎng)格線(xiàn)就會(huì)出現(xiàn)嚴(yán)重傾斜,甚至出現(xiàn)負(fù)網(wǎng)格的情況,導(dǎo)致計(jì)算發(fā)散。鑒于Fluent軟件在表征流體流動(dòng)形態(tài)方面的優(yōu)勢(shì)及其計(jì)算油膜軸承動(dòng)態(tài)性能方面的不足,本文提出了一種基于Fluent方法的滑動(dòng)軸承油膜性能計(jì)算的動(dòng)網(wǎng)格更新方法。在它的每一時(shí)間步內(nèi),所有網(wǎng)格節(jié)點(diǎn)更新方式是依據(jù)軸頸中心坐標(biāo)確定的,使徑向網(wǎng)格線(xiàn)始終垂直于軸頸表面,保證了移動(dòng)前后膜厚方向的網(wǎng)格線(xiàn)沿圓周方向上等均分布,網(wǎng)格不發(fā)生扭曲變形,避免了網(wǎng)格節(jié)點(diǎn)更新產(chǎn)生的累計(jì)誤差。通過(guò)與典型算例和實(shí)驗(yàn)結(jié)果對(duì)比分析來(lái)驗(yàn)證該方法的有效性,并考察了進(jìn)油壓力和載荷對(duì)轉(zhuǎn)子靜平衡位置的影響。
動(dòng)網(wǎng)格模型可以用來(lái)模擬由于邊界運(yùn)動(dòng)而引起的流域隨時(shí)間的變化。邊界運(yùn)動(dòng)可以是主動(dòng)的規(guī)定動(dòng)作(例如,用戶(hù)可以定義剛體中心的速度和角速度),也可以是被動(dòng)的待求解的動(dòng)作(例如,六自由度模型中,已知運(yùn)動(dòng)狀態(tài)的剛體邊界,受合力后確定下一時(shí)刻的速度和加速度)。根據(jù)每一個(gè)時(shí)間步的運(yùn)動(dòng)邊界的位置,F(xiàn)luent軟件自動(dòng)對(duì)動(dòng)網(wǎng)格更新處理。對(duì)于動(dòng)網(wǎng)格模型,用戶(hù)需要確定初始網(wǎng)格和運(yùn)動(dòng)區(qū)域。3種常用的動(dòng)網(wǎng)格更新算法[18]分別是Smoothing、Layering和Remeshing。Smoothing方法是不改變網(wǎng)格節(jié)點(diǎn)間的拓?fù)潢P(guān)系,只進(jìn)行網(wǎng)格形狀的改變,但當(dāng)網(wǎng)格出現(xiàn)畸變時(shí)計(jì)算容易發(fā)散。Layering方法是隨著動(dòng)邊界的移動(dòng),在邊界處發(fā)生網(wǎng)格的增加或合并。Remeshing方法是將控制區(qū)內(nèi)的所有網(wǎng)格重新劃分,一般適用于非結(jié)構(gòu)化網(wǎng)格。理論上,考慮到油膜幾何尺寸的特點(diǎn),應(yīng)當(dāng)選擇Smoothing動(dòng)網(wǎng)格更新方式,但在軸承性能計(jì)算過(guò)程中也出現(xiàn)了網(wǎng)格畸變,顯然Smoothing方法也不能很好地解決軸承油膜的瞬態(tài)計(jì)算問(wèn)題。
文獻(xiàn)[13]提出的基于Smoothing的網(wǎng)格更新模型示意圖見(jiàn)圖1。網(wǎng)格節(jié)點(diǎn)位移按公式Δxi=niΔx/n和 Δyi=niΔy/n進(jìn)行計(jì)算(n為膜厚方向網(wǎng)格的層數(shù),ni為網(wǎng)格節(jié)點(diǎn)所在的層數(shù),Δx、Δy為軸頸中心在該時(shí)間步的位移)。利用該方法,垂直于軸頸表面的網(wǎng)格線(xiàn)就會(huì)出現(xiàn)傾斜,導(dǎo)致計(jì)算結(jié)果誤差較大,特別是多步計(jì)算后的位置累計(jì)誤差。
圖1 文獻(xiàn)[13]網(wǎng)格更新方式Fig.1 Mesh updating method in references
為解決上述難點(diǎn),本文提出了一種新的網(wǎng)格更新算法,在網(wǎng)格更新前后,膜厚方向的網(wǎng)格線(xiàn)始終指向軸頸的中心,從而避免了出現(xiàn)網(wǎng)格傾斜。本文網(wǎng)格結(jié)構(gòu)更新算法的示意圖見(jiàn)圖2,網(wǎng)格節(jié)點(diǎn)位置計(jì)算原理見(jiàn)圖3。
圖2 本文提出的新的網(wǎng)格更新方式Fig.2 New mesh updating method in this paper
圖3 網(wǎng)格節(jié)點(diǎn)更新原理圖Fig.3 Schematic diagram of mesh node updating
為使網(wǎng)格不發(fā)生扭曲變形,就需要指向軸頸中心的網(wǎng)格線(xiàn)方向始終不改變,即始終垂直于軸頸表面,同時(shí)要保證在移動(dòng)前后膜厚方向的網(wǎng)格線(xiàn)沿圓周方向上等均分布(膜厚網(wǎng)格線(xiàn)的夾角相等),即更新前后網(wǎng)格線(xiàn)平行,也就是圖3中的向量O1N1與O2N2共線(xiàn)。具體方法如下。
記軸承的中心為O,移動(dòng)前后的軸頸的中心分別為O1和O2,油膜網(wǎng)格在膜厚方向上分成k層網(wǎng)格,假設(shè)N1為t1時(shí)刻網(wǎng)格上任意一個(gè)待更新節(jié)點(diǎn),節(jié)點(diǎn)N1以及與其相對(duì)應(yīng)的更新后的節(jié)點(diǎn)N2均在第ki層網(wǎng)格上,根據(jù)此刻軸頸所受的油膜力、軸頸坐標(biāo)、軸頸中心位移速度、軸頸質(zhì)量可以得到下一時(shí)刻t2的軸頸中心坐標(biāo),根據(jù)時(shí)刻t2的軸頸中心坐標(biāo)和該節(jié)點(diǎn)相對(duì)軸心的單位向量(與向量O1N1平行或者共線(xiàn)),作射線(xiàn)分別交軸頸圓和軸承圓于點(diǎn)A2、B2。又知待求點(diǎn)N2所在的層數(shù)為ki,即已知LA2N2=(ki/k)LA2B2,LOB2為軸承半徑,LO2A2為軸頸半徑。N2坐標(biāo)求解算法步驟如下。
(1)求出LOO2;
(2)已知向量O1N1的單位法向量e,由于向量O1N1與向量O2N2共線(xiàn),故可得到向量O2N2的單位法向量e;
(3)求解向量O2O與向量O2B2的夾角余弦cosα;
(4)利用三角形關(guān)系式求解LO2B2:
(5)求出LA2B2=LO2B2-LO2A2;
已知軸頸中心移動(dòng)后的坐標(biāo)O2,根據(jù)上述算法可求出節(jié)點(diǎn)位置N1移動(dòng)到新位置的坐標(biāo)N2。
根據(jù)牛頓第二定律,運(yùn)動(dòng)方程[13]為
(1)
式中,m為轉(zhuǎn)子質(zhì)量的一半;S為轉(zhuǎn)子中心距原點(diǎn)的位移,G=mg;F為轉(zhuǎn)子受的油膜合力。
已知時(shí)間步Δt,上一時(shí)間步對(duì)應(yīng)的速度v0,根據(jù)牛頓第二定律求解每一個(gè)時(shí)間點(diǎn)軸頸中心的相對(duì)位移和速度計(jì)算公式如下:
(2)
v=v0+(G+F)Δt/m
(3)
求解油膜力F的分量Fx和Fy的公式如下:
(4)
(5)
式中,F(xiàn)x、Fy分別為x軸和y軸方向的油膜力;A為某一網(wǎng)格單元的面積;p為根據(jù)N-S方程計(jì)算出的軸承油膜各點(diǎn)的壓力。
基于Fluent軟件,利用該方法計(jì)算油膜軸承所支撐軸頸靜平衡位置的流程如圖4所示,首先設(shè)定軸頸初始速度和初始位移,通過(guò)udf宏DEFINE_GRID_MOTION計(jì)算出油膜力Fx和Fy,再根據(jù)式(2)和式(3)計(jì)算出軸心坐標(biāo)和速度。同時(shí)按照新提出的算法原理進(jìn)行網(wǎng)格更新,在該時(shí)間步內(nèi)迭代完后進(jìn)入下一個(gè)時(shí)間步。下一個(gè)時(shí)間步用到上一時(shí)間步的軸心坐標(biāo)和速度(讀取儲(chǔ)存于文件2和文件3最后一次更新的數(shù)值)。一直循環(huán)下去,直到軸頸中心趨近于一定點(diǎn),即可停止計(jì)算,該定點(diǎn)即穩(wěn)態(tài)下徑向油膜軸承所支撐軸頸中心的靜平衡位置。
圖4 軸心靜平衡位置計(jì)算流程圖Fig.4 Computational flow chart about static balance position of the rotor
一般來(lái)說(shuō),動(dòng)網(wǎng)格節(jié)點(diǎn)的位移由坐標(biāo)確定,但滑動(dòng)軸承油膜是一種長(zhǎng)徑比很大的網(wǎng)格,用坐標(biāo)判斷節(jié)點(diǎn)位置容易出現(xiàn)錯(cuò)誤,可采用節(jié)點(diǎn)全局編號(hào)判定節(jié)點(diǎn)的位置,強(qiáng)制精確控制每一個(gè)節(jié)點(diǎn)的位移變化。利用DEFINE_GRID_MOTION宏對(duì)網(wǎng)格節(jié)點(diǎn)位置進(jìn)行強(qiáng)制性精確定義,可有效降低網(wǎng)格畸變的可能性,由于軸頸表面各點(diǎn)的速度不一樣,可通過(guò)編寫(xiě)DEFINE_PROFILE宏程序控制軸頸表面各點(diǎn)的速度。按照以上步驟,可以利用動(dòng)網(wǎng)格方法實(shí)現(xiàn)對(duì)滑動(dòng)軸承油膜狀態(tài)的仿真。
本文提出的網(wǎng)格更新方法與文獻(xiàn)[13]所使用的網(wǎng)格更新方法的主要區(qū)別就是更新后油膜厚度方向的網(wǎng)格線(xiàn)是否垂直于軸頸表面。要驗(yàn)證所提動(dòng)網(wǎng)格算法的正確性和有效性,需要進(jìn)行兩方面的驗(yàn)證:用于靜態(tài)油膜力計(jì)算的精度驗(yàn)證;用于軸頸中心靜平衡位置求解的累積誤差驗(yàn)證。
利用所提出的網(wǎng)格結(jié)構(gòu)和文獻(xiàn)[13-15]的網(wǎng)格結(jié)構(gòu),分別計(jì)算軸頸轉(zhuǎn)動(dòng)角速度為500 rad/s、偏心率為0.5和0.9的圓柱軸承的油膜力,來(lái)對(duì)比這兩種網(wǎng)格在靜態(tài)油膜力計(jì)算(用于靜網(wǎng)格計(jì)算時(shí),單次網(wǎng)格計(jì)算)的精度。算例選取的軸承幾何參數(shù)和潤(rùn)滑油物性參數(shù)見(jiàn)表1。軸承供油方式采用兩側(cè)雙向進(jìn)油并且設(shè)置軸向油槽。
由于油膜厚度尺寸很小,故油膜厚度方向的網(wǎng)格劃分對(duì)計(jì)算結(jié)果影響最大,需要先進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證。網(wǎng)格獨(dú)立性計(jì)算時(shí)選擇多項(xiàng)流混合模型,氣穴模型選用Singhal全空化模型,連續(xù)性方程、動(dòng)量方程和氣穴方程均選用一階迎風(fēng)格式進(jìn)行求解,速度耦合格式選擇SIMPILE格式。軸
表1 軸承和潤(rùn)滑油物性參數(shù)
承油膜網(wǎng)格初始劃分的網(wǎng)格數(shù)(膜厚方向×周向×軸向)為5×600×80,通過(guò)改變膜厚方向的網(wǎng)格層數(shù),并以?xún)袅髁啃∮谶M(jìn)口和出口流量?jī)烧咧g的最小值的1%作為判斷收斂的標(biāo)準(zhǔn)。膜厚方向上設(shè)置不同的網(wǎng)格層數(shù),以承載力作為觀測(cè)量,網(wǎng)格獨(dú)立性驗(yàn)證結(jié)果見(jiàn)表2。結(jié)果表明,本文所用的網(wǎng)格密度計(jì)算結(jié)果的最大偏差均小于4%。
表2 網(wǎng)格獨(dú)立性驗(yàn)證結(jié)果
綜合考慮計(jì)算效率和計(jì)算精度,把油膜膜厚方向網(wǎng)格劃分為6層(6×600×80),見(jiàn)圖5。算例中壓力-速度耦合格式分別選SIMPILE、SIMPILEC、PISO方式,其他設(shè)置同上。采用瞬態(tài)計(jì)算方法,利用文獻(xiàn)中的網(wǎng)格方法和本文的網(wǎng)格方法計(jì)算得到的承載力結(jié)果和圓周方向的壓力分布結(jié)果(取軸承寬度中間L/2處,偏心率為0.9)分別見(jiàn)表3和圖6。
圖5 油膜網(wǎng)格Fig.5 Oil film mesh
表3 兩種靜網(wǎng)格結(jié)構(gòu)計(jì)算的承載力對(duì)比
(b) 本文方法的計(jì)算結(jié)果
由表3和圖6可以看出,基于兩種網(wǎng)格方法得到的軸承承載力和瓦塊圓周方向的壓力計(jì)算結(jié)果相同,均與文獻(xiàn)[19]結(jié)果相接近,表明兩種網(wǎng)格(油膜厚度方向的網(wǎng)格線(xiàn)垂直于或者不垂直于軸頸表面)用于軸承油膜穩(wěn)態(tài)油膜力計(jì)算時(shí),均能滿(mǎn)足精度要求(<9%)。由于在計(jì)算過(guò)程中網(wǎng)格不更新(靜網(wǎng)格法),無(wú)論油膜厚度方向的網(wǎng)格線(xiàn)是否傾斜,累計(jì)計(jì)算誤差都不存在,故兩種網(wǎng)格的計(jì)算結(jié)果基本無(wú)差別。這就證明了Smoothing網(wǎng)格方法和本文提出的網(wǎng)格方法來(lái)計(jì)算軸承穩(wěn)態(tài)油膜力(靜網(wǎng)格法)的精度是滿(mǎn)足要求的。
根據(jù)文獻(xiàn)[13]的算例,選取軸頸角速度為500 rad/s,轉(zhuǎn)子質(zhì)量為17.527 kg(171.76 N),軸承、潤(rùn)滑油的參數(shù)和油膜網(wǎng)格劃分同算例1,對(duì)軸頸中心軌跡的收斂過(guò)程及靜平衡位置進(jìn)行了仿真,結(jié)果見(jiàn)圖7,計(jì)算得到其偏心率為0.221。根據(jù)文獻(xiàn)[19](未給出進(jìn)油壓力,其他條件完全一樣),在相同軸承參數(shù)下,據(jù)此載荷計(jì)算出軸承偏心率為0.267,與本文計(jì)算結(jié)果相比,相對(duì)誤差是17.2%。文獻(xiàn)[13]和本文在計(jì)算時(shí)考慮了氣穴的影響,而文獻(xiàn)[19]沒(méi)有考慮氣穴的影響。
(a) 二維圖
(b) 三維圖
文獻(xiàn)[13]的靜平衡位置所對(duì)應(yīng)的偏心率為0.046。將文獻(xiàn)[13]結(jié)果和本文的計(jì)算結(jié)果分別與文獻(xiàn)[18]相比較,文獻(xiàn)[13]的計(jì)算結(jié)果差別更大,偏心率明顯不在同一數(shù)量級(jí),表明文獻(xiàn)[13]的網(wǎng)格經(jīng)多次更新后,出現(xiàn)了較大的累計(jì)誤差。這與文獻(xiàn)[13]在網(wǎng)格更新時(shí)沿膜厚方向的網(wǎng)格線(xiàn)與軸頸表面不垂直有關(guān),也與其網(wǎng)格節(jié)點(diǎn)坐標(biāo)更新的計(jì)算方式與軸頸表面的位置有關(guān),從而導(dǎo)致多次迭代后求解軸頸中心靜平衡位置時(shí)出現(xiàn)了較大的誤差。
為驗(yàn)證所編制程序的穩(wěn)定性,從3個(gè)不同起始位置開(kāi)始分別計(jì)算相同條件下的軸心靜平衡位置,軸心軌跡收斂曲線(xiàn)如圖8所示??梢钥闯?,3條曲線(xiàn)最終收斂于一點(diǎn),表明計(jì)算的起始位置對(duì)軸頸的靜平衡位置并無(wú)影響,驗(yàn)證了所編制的網(wǎng)格更新程序的穩(wěn)定性。
圖8 不同起點(diǎn)位置的軸心收斂軌跡Fig.8 Axial convergence trajectory of different starting positions
為了分析進(jìn)油壓力對(duì)軸承所支撐轉(zhuǎn)子軸心靜平衡位置的影響,在給定氣化壓力29 185 Pa、轉(zhuǎn)子質(zhì)量17.527 kg、角速度500 rad/s和 000 rad/s的條件下,計(jì)算了供油壓力在0.1~0.5 MPa之間時(shí)轉(zhuǎn)子軸心的靜平衡位置。軸承的偏心率和偏位角隨進(jìn)油壓力的變化關(guān)系如圖9所示。算例的其他計(jì)算參數(shù)見(jiàn)表1。
(a) 偏心率
(b) 偏位角
由圖9可以看出,隨著進(jìn)油壓力的增大,偏心率減小(<7%),偏位角增大(<9%)。這是因?yàn)樵谳d荷和轉(zhuǎn)速不變的前提下,隨著進(jìn)油壓力的增大,軸承油膜的最大壓力增大,承載能力增強(qiáng),偏心率減小。同時(shí)由于進(jìn)油口在水平方向上,故隨著進(jìn)油壓力增大,偏位角增大。計(jì)算結(jié)果表明進(jìn)油壓力對(duì)軸承的承載性能有影響,但影響率小于10%。結(jié)合算例2的結(jié)果,進(jìn)一步證明了所提出的動(dòng)網(wǎng)格更新方法用于軸心軌跡計(jì)算時(shí)累計(jì)誤差相對(duì)較小。
通過(guò)以上算例驗(yàn)證了所提出的動(dòng)網(wǎng)格方法在油膜性能計(jì)算的有效性、可行性和穩(wěn)定性。
(1)本文提出了一種用于滑動(dòng)軸承油膜性能計(jì)算的動(dòng)網(wǎng)格更新方法,驗(yàn)證了其正確性。本文方法能夠保證在網(wǎng)格移動(dòng)前后膜厚方向的網(wǎng)格線(xiàn)沿圓周方向上等均分布,并使指向軸頸中心的網(wǎng)格線(xiàn)方向始終垂直于軸頸表面,網(wǎng)格不發(fā)生傾斜。
(2) 本文方法應(yīng)用于求解滑動(dòng)軸承所支撐轉(zhuǎn)子的靜平衡位置時(shí),能夠減少網(wǎng)格計(jì)算的累計(jì)誤差,提高計(jì)算的精度;同時(shí)發(fā)現(xiàn)利用Fluent軟件動(dòng)網(wǎng)格方法計(jì)算油膜軸承流場(chǎng)時(shí),膜厚方向網(wǎng)格傾斜會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生較大的累計(jì)誤差。
(3)在相同的條件下,隨著進(jìn)油壓力增大,軸承偏心率減小,偏位角增大,但進(jìn)油壓力對(duì)軸承的承載性能的影響率一般小于10%。