基于接觸非線性的箔片結(jié)構(gòu)力學(xué)特性分析

摘要：箔片軸承中頂箔與波箔之間的非線性接觸力學(xué)問(wèn)題具有非線性和時(shí)變性，建立考慮非線性接觸的理論模型研究箔片結(jié)構(gòu)的力學(xué)行為，采用增量迭代法和牛頓-拉弗森法求解系統(tǒng)平衡方程，通過(guò)該模型進(jìn)行轉(zhuǎn)子推拉仿真獲得箔片結(jié)構(gòu)變形量與響應(yīng)力之間的關(guān)系以及箔片的能量耗散。結(jié)果表明:摩擦是箔片結(jié)構(gòu)產(chǎn)生滯回行為的主要因素，隨著摩擦因數(shù)的增大，箔片結(jié)構(gòu)剛度增大，滯回曲線面積也隨之增大；在一個(gè)加載、卸載回路中，加載過(guò)程比卸載過(guò)程具有更高的靜載力；靠近固定端的波箔片比靠近自由端的波箔片具有更高的剛度。關(guān)鍵詞：滑動(dòng)軸承;箔片軸承;氣體動(dòng)壓軸承;滯回曲線;能量耗散與傳統(tǒng)的油潤(rùn)滑軸承和滾動(dòng)軸承相比,箔片氣體軸承因結(jié)構(gòu)簡(jiǎn)單,使用壽命長(zhǎng),維護(hù)成本低以及能在高速、高溫等極端環(huán)境下運(yùn)行[1-2],成為高速渦輪機(jī)械的絕佳選擇[3]。在已開發(fā)的箔片軸承類型中,徑向箔片氣體軸承使用最廣泛,軸套內(nèi)部放置一個(gè)柔性箔結(jié)構(gòu),該結(jié)構(gòu)由一個(gè)波箔和一個(gè)頂箔組成,這種彈性支承結(jié)構(gòu)增強(qiáng)了軸承的承載能力和阻尼[4],克服了氣體黏度低帶來(lái)的性能缺陷?？焖?、準(zhǔn)確預(yù)測(cè)軸承性能對(duì)于箔片氣體軸承的開發(fā)至關(guān)重要,然而,由于摩擦使得軸承系統(tǒng)具有高度非線性,仍然很難準(zhǔn)確描述軸承的工作特性。文獻(xiàn)[5-6]用一種彈簧模型模擬箔片,將各個(gè)波箔片等效為無(wú)相互作用的彈簧,但忽略了頂箔的剛度,箔片之間的摩擦和波箔之間的相互作用,導(dǎo)致結(jié)構(gòu)剛度被低估。不過(guò)該模型因其結(jié)構(gòu)簡(jiǎn)單,易于編程仍被廣泛使用[7-9],文獻(xiàn)[10]就利用此類模型在研究中首次考慮了箔片結(jié)構(gòu)的阻尼效應(yīng),并預(yù)測(cè)了箔片軸承的剛度和阻尼系數(shù)。雖然上述模型已被廣泛采用,但由于轉(zhuǎn)軸與頂箔之間碰撞的相互作用被丟棄,因而具有局限性。為獲得更精確的模擬結(jié)果,開始使用有限元法對(duì)軸承進(jìn)行建模:文獻(xiàn)[11]使用二維梁?jiǎn)卧獙?duì)頂箔和波箔建模,并使用黏滑算法[12]考慮摩擦,建立了考慮滯回特性的模型,預(yù)測(cè)表明一個(gè)靜載荷對(duì)應(yīng)多個(gè)平衡位置;文獻(xiàn)[13]同樣采用梁?jiǎn)卧獮榛A(chǔ)的有限元結(jié)構(gòu)模型,但不同之處在于采用非線性接觸的數(shù)值程序處理各接觸區(qū)域的接觸力,并將該方法應(yīng)用于三維模型[14]。另外還可使用殼單元模擬箔片[15-18],如文獻(xiàn)[15]提出了一個(gè)考慮大位移的非線性模型,并將箔片建模為殼單元,文獻(xiàn)[18]以懲罰法和正則化摩擦定律為基礎(chǔ)對(duì)波箔進(jìn)行建模。還有使用有限元商業(yè)代碼[19-20]對(duì)箔片進(jìn)行建模的方法,其考慮了箔片之間的相互作用,但求解較為復(fù)雜,在建立非線性摩擦模型時(shí)數(shù)值收斂性不好。為簡(jiǎn)化計(jì)算量并保證計(jì)算精度,開發(fā)了近似有限元結(jié)構(gòu)的桁架系統(tǒng)。文獻(xiàn)[21]用一種結(jié)構(gòu)簡(jiǎn)單的桁架系統(tǒng)模擬波箔片,每個(gè)波箔片被建模為3個(gè)基本彈簧,彈簧的剛度由卡斯提利亞諾理論獲得,雖然結(jié)構(gòu)簡(jiǎn)單,但與有限元仿真結(jié)果表現(xiàn)出更好的一致性,還證明了波箔相互作用的重要性。文獻(xiàn)[22-23]采用上述桁架模型,在波箔與軸套的接觸部位增加了垂直自由度,并考慮了頂箔的影響和軸承系統(tǒng)的接觸/分離行為,箔片結(jié)構(gòu)的法向和切向接觸力分別采用增廣拉格朗日乘子法和懲罰法處理,分析結(jié)果與文獻(xiàn)[24-25]一致。為準(zhǔn)確描述摩擦對(duì)軸承系統(tǒng)特性的影響,本文提出一種考慮非線性接觸的靜態(tài)結(jié)構(gòu)模型,其能捕捉柔性結(jié)構(gòu)的非線性與時(shí)變性;主要討論波箔與軸套之間以及箔片之間的點(diǎn)線接觸狀態(tài);利用懲罰法離散每個(gè)接觸區(qū)域的力學(xué)向量,通過(guò)正則化技術(shù)改進(jìn)庫(kù)倫摩擦模型,采用牛頓-拉弗森法求解平衡方程;最后通過(guò)轉(zhuǎn)子推拉仿真得到箔片結(jié)構(gòu)的力學(xué)行為與能量耗散等特性。1理論模型1.1波箔建模以軸心為原點(diǎn)建立全局坐標(biāo)系OXY,則徑向箔片氣體軸承簡(jiǎn)化圖如圖1所示,為簡(jiǎn)化計(jì)算量,波箔采用文獻(xiàn)[21]提出的簡(jiǎn)化桁架結(jié)構(gòu)進(jìn)行建模。對(duì)波箔建立局部坐標(biāo)系Oxy,2個(gè)波箔片組成的箔片數(shù)學(xué)模型如圖2所示,可以擴(kuò)展到任意數(shù)量波箔片,因?yàn)榭紤]了波箔、頂箔與軸套之間的松緊接觸,在每個(gè)凸塊底部都需要增加一個(gè)額外的自由度(圖2中藍(lán)色箭頭),即垂直位移,并且需要重新求解波箔的剛度。圖1徑向箔片氣體軸承簡(jiǎn)化圖Fig.1Simplifieddiagramofradialfoilgasbearing圖22個(gè)波箔片組成的箔片數(shù)學(xué)模型Fig.2Amathematicalmodeloffoilcomposedoftwobumpfoils彈簧剛度與文獻(xiàn)[21]中一致,每個(gè)彈簧均被視為基本單元。由于彈簧剛度ki已經(jīng)量化,使用直接剛度法更容易求解單元的剛度矩陣,因此局部坐標(biāo)系中單元的剛度矩陣可表示為(1)進(jìn)一步轉(zhuǎn)化為全局坐標(biāo)系的局部剛度矩陣kb,即(2)(3)c=cos

β,s=sin

β,式中:T為坐標(biāo)轉(zhuǎn)換矩陣;β為兩坐標(biāo)系之間的夾角。1.2頂箔建模頂箔由歐拉-伯努利梁?jiǎn)卧M(jìn)行建模,每個(gè)波箔片對(duì)應(yīng)的頂箔被離散為3個(gè)單元,每個(gè)單元有2個(gè)節(jié)點(diǎn)和3個(gè)自由度,展開的頂箔單元水平放置在桁架結(jié)構(gòu)上,因此頂箔單元的全局坐標(biāo)系與局部坐標(biāo)系重合。在全局坐標(biāo)系中,頂箔單元的剛度矩陣ktop表示為頂箔的局部單元?jiǎng)偠染仃囉少み|金法獲得,即(4)式中:E為材料彈性模量,GPa;A為單元橫截面積,m2;I為橫截面積的慣性矩,m4;Le為單元長(zhǎng)度,m。箔片的整體剛度矩陣K=diag(Kb,Ktop)通過(guò)單元節(jié)點(diǎn)的編號(hào)組合單元?jiǎng)偠染仃嚝@得,其中Kb=∑kb,Ktop=∑ktop。全局結(jié)構(gòu)的廣義位移坐標(biāo)向量其中ub為波箔單元節(jié)點(diǎn)的位移坐標(biāo)向量,ut為頂箔單元節(jié)點(diǎn)的位移坐標(biāo)向量。1.3摩擦接觸建模頂箔和軸套視為主體,波箔視為從體。在外力作用下,箔片的變形增加了箔片之間的接觸面積,理論上,每個(gè)波箔與頂箔之間應(yīng)建立多個(gè)接觸點(diǎn)對(duì),由于與箔片的長(zhǎng)度相比變形的尺寸較小,因此僅選擇每個(gè)波箔上的最高節(jié)點(diǎn)作為接觸節(jié)點(diǎn)。箔片結(jié)構(gòu)及接觸點(diǎn)示意圖如圖3所示,波箔上的綠色節(jié)點(diǎn)為接觸節(jié)點(diǎn),也被視為檢測(cè)點(diǎn),被檢測(cè)節(jié)點(diǎn)位于紅色單元中,藍(lán)色節(jié)點(diǎn)是與接觸對(duì)相關(guān)的節(jié)點(diǎn)。圖3箔片結(jié)構(gòu)及接觸點(diǎn)示意圖Fig.3Diagramoffoilstructureandcontactpoints頂箔與波箔之間的接觸示意圖如圖4所示,Q為波箔的接觸節(jié)點(diǎn),P為Q在頂箔上的投影,P與Q之間的相對(duì)位移為(下角標(biāo)t表示切向方向,n表示法向方向)圖4頂箔與波箔之間的接觸示意圖Fig.4DiagramofcontactbetweentopfoilandbumpfoiluP-uQ=Ncuc,(5)Nc=[-N1(ξ)-N2(ξ)

I2×2],式中:uP為P點(diǎn)位移向量,需要通過(guò)節(jié)點(diǎn)P1和P2的位移向量插值得到;uP,uP分別為節(jié)點(diǎn)P1,P2的位移向量,為轉(zhuǎn)角向量;uQ為Q點(diǎn)的位移向量;Nc為形狀函數(shù),利用插值法得到相應(yīng)坐標(biāo);uc為形狀位移;Ni(ξ)為梁?jiǎn)卧男魏瘮?shù);ξ為P點(diǎn)對(duì)應(yīng)量綱一的坐標(biāo),ξ(P1)=-1,ξ(P2)=1;I為單位矩陣。上述向量在全局坐標(biāo)系中定義,為方便計(jì)算,將其引入到接觸條件中,在接觸點(diǎn)對(duì)的局部坐標(biāo)系中重新定義(5)式,即(6)L=-I2×2,(7)式中:uB,uA分別為P,Q在局部坐標(biāo)系中的位移向量;L為坐標(biāo)變換矩陣。將接觸面離散化后,每個(gè)接觸區(qū)域的分布接觸力轉(zhuǎn)化為離散形式,節(jié)點(diǎn)Q的法向接觸力Fn表示為(8)式中:α為懲罰數(shù),實(shí)際計(jì)算中不可能將α設(shè)置為無(wú)窮大以獲得準(zhǔn)確解,需要根據(jù)實(shí)際問(wèn)題選擇一個(gè)可接受的值,本文α取1.0×108

N/mm。利用庫(kù)倫摩擦模型將節(jié)點(diǎn)Q的切向接觸力Ft表示為(9)式中:分別為當(dāng)前加載步驟與前一個(gè)加載步驟的切向相對(duì)位移,為摩擦因數(shù)。時(shí),接觸點(diǎn)處于黏著狀態(tài),時(shí),接觸點(diǎn)在大小和方向上存在不受限制的相對(duì)滑動(dòng)。由于接觸狀態(tài)的不連續(xù),數(shù)值計(jì)算難以收斂,因此使用符號(hào)函數(shù)來(lái)近似平滑不連續(xù),即(10)式中:c1為控制正則化摩擦模型與庫(kù)倫摩擦模型之間接近度的平滑參數(shù),本文c1取100m。結(jié)合(9)—(10)式得到Q點(diǎn)切向力為(11)根據(jù)反作用力原理可得到相關(guān)節(jié)點(diǎn)對(duì)和接觸節(jié)點(diǎn)對(duì)的等效節(jié)點(diǎn)接觸力向量,即(12)FB=[Ft

Fn]T,式中:FB為B點(diǎn)的接觸力向量;qc等效節(jié)點(diǎn)接觸力向量。由于考慮了摩擦接觸,非線性行為被引入模型,假設(shè)將所有單元組合后形成的系統(tǒng)的全局非線性方程為K(u)=FL=Ku-Qc,(13)式中:K(u)為u的非線性函數(shù);FL為外力向量。牛頓-拉弗森法是求解非線性方程的基本方法之一,其從一個(gè)假設(shè)解開始,通過(guò)不斷的迭代逼近真實(shí)解,直到滿足收斂準(zhǔn)則。假設(shè)第l步的近似解已知,記作ul,則下一步的近似解ul+1近似為一階泰勒展開式,即(14)進(jìn)一步寫為(15)(16)Qc=[(qc,1)T

(qc,2)T

…(qc,N)T]T,(17)D=[(uc,1)T

(uc,2)T

…(uc,N)T]T,(18)式中:Δul為位移向量的增量;為第l步的全局切向剛度矩陣;Qc,D分別為所有接觸點(diǎn)對(duì)與相關(guān)節(jié)點(diǎn)的接觸力向量和位移向量?？梢酝ㄟ^(guò)使用相應(yīng)的選擇矩陣從全局向量中提取。任一接觸點(diǎn)對(duì)的等效節(jié)點(diǎn)接觸力qc的導(dǎo)數(shù)為(19)最后得到線性化的系統(tǒng)方程為(20)(21)通過(guò)(20)式得到增量Δul,則新的近似解為ul+1=ul+Δul,通常該解滿足不了原非線性方程,當(dāng)殘差變得足夠小時(shí),u是正確解。平衡方程的收斂準(zhǔn)則為(22)2結(jié)果與分析轉(zhuǎn)子不對(duì)中時(shí),隨著轉(zhuǎn)軸在垂直方向移動(dòng),箔片結(jié)構(gòu)在垂直和水平方向都發(fā)生變形,箔片的垂直變形決定了軸承的承載能力,而水平運(yùn)動(dòng)產(chǎn)生阻尼。箔片氣體軸承結(jié)構(gòu)參數(shù)見(jiàn)表1,利用Matlab將數(shù)學(xué)模型轉(zhuǎn)化為程序語(yǔ)言進(jìn)行仿真,轉(zhuǎn)軸的最大徑向位移設(shè)為79.5μm(即2.5C),程序模擬轉(zhuǎn)軸在正負(fù)X,Y方向上逐步壓向箔片結(jié)構(gòu),當(dāng)轉(zhuǎn)軸位移超過(guò)C時(shí),轉(zhuǎn)軸的表面和箔片結(jié)構(gòu)開始相互作用,達(dá)到最大徑向位移后,轉(zhuǎn)軸逐漸被拉回初始位置。表1箔片氣體軸承結(jié)構(gòu)參數(shù)Tab.1Structuralparametersoffoilgasbearings2.1箔片結(jié)構(gòu)力學(xué)行為不同摩擦因數(shù)下X,Y方向靜載荷與轉(zhuǎn)軸位移的關(guān)系如圖5所示:除-X方向外仿真結(jié)果與文獻(xiàn)[20]完全非線性模型的結(jié)果非常吻合,因?yàn)榉抡娉绦蜿U明了加載和卸載,所以預(yù)期的滯回曲線清晰可見(jiàn);靜載荷為轉(zhuǎn)軸擠壓頂箔時(shí)受到頂箔的反作用力,方向?yàn)檗D(zhuǎn)軸徑向移動(dòng)方向,并與接觸點(diǎn)處切線垂直,加載過(guò)程的靜載荷比卸載過(guò)程的靜載荷大,這是因?yàn)樾遁d開始時(shí)一部分力平衡了摩擦力,接觸點(diǎn)不會(huì)立即滑動(dòng);載荷與偏轉(zhuǎn)曲線在輕載荷區(qū)域不是線性的,而在重載荷區(qū)域幾乎是直線,這是因?yàn)樵谳p載荷區(qū)域只有少數(shù)波箔發(fā)生變形,當(dāng)變形波箔的數(shù)量不再增加時(shí),載荷與偏轉(zhuǎn)曲線的斜率幾乎恒定;隨著摩擦因數(shù)的增大,環(huán)路變得越來(lái)越明顯,當(dāng)沒(méi)有摩擦?xí)r,加載過(guò)程曲線與卸載過(guò)程曲線完全重合;在-X方向上,仿真預(yù)測(cè)結(jié)果明顯小于完全非線性模型預(yù)測(cè)結(jié)果,文獻(xiàn)[22]指出,雖然固定端附近的幾個(gè)波箔的剛度非常大,但自由端附近的波箔剛度非常小,因此完全非線性模型可能會(huì)高估-X方向上的結(jié)構(gòu)剛度。

圖5不同摩擦因數(shù)下靜載荷與轉(zhuǎn)軸位移的關(guān)系Fig.5Relationshipbetweenstaticloadandshaftdisplace-mentunderdifferentfrictioncoefficients摩擦因數(shù)分別為0,0.10,轉(zhuǎn)軸達(dá)到最大徑向位移(2.5C)時(shí)4個(gè)方向波箔上的靜載荷仿真結(jié)果如圖6所示:波箔從固定端開始編號(hào),摩擦因數(shù)對(duì)各波箔之間剛度的相對(duì)大小有影響;對(duì)于-X方向,靠近固定端的第1和第2個(gè)波箔比其他波箔的靜載荷大得多,自由端附近波箔的靜載荷很小。

圖6轉(zhuǎn)軸位移為2.5C時(shí)各波箔的承載力Fig.6Loadcapacityofeachbumpfoilwhenshaftdisplacementis2.5C2.2固定端對(duì)箔片剛度的影響選擇+X方向分析固定端對(duì)箔片剛度的影響,如圖7所示,隨著轉(zhuǎn)軸在+X方向上被推拉(最大徑向位移為2.5C),第10到第18個(gè)波箔與轉(zhuǎn)軸表面接觸,波箔14的垂直位移最大,波箔14左右兩側(cè)各選一個(gè)波箔(波箔12與16)進(jìn)行對(duì)比分析。圖7波箔與轉(zhuǎn)軸表面接觸Fig.7Contactsofbumpfoilswithshaft轉(zhuǎn)軸在+X方向上被推拉至徑向位移為2.5C的過(guò)程中,波箔12,14,16的滯回曲線如圖8所示:無(wú)論加載還是卸載過(guò)程,靠近固定端的波箔具有更大的剛度;在整個(gè)回路中,波箔垂直位移與靜載荷之間呈線性關(guān)系,符合實(shí)際情況。轉(zhuǎn)軸最大徑向位移分別為60,70,80μm時(shí)波箔14的滯回曲線如圖9所示:在卸載過(guò)程中曲線不重合,這是固定端施加反作用力的結(jié)果。通過(guò)上述仿真結(jié)果可以看出,固定端對(duì)箔片結(jié)構(gòu)力學(xué)性能的影響不能忽略。圖8波箔12,14,16的滯回曲線Fig.8Hystereticcurvesofbumpfoils12,14,and16圖9轉(zhuǎn)軸最大徑向位移分別為60,70,80μm時(shí)波箔14的滯回曲線Fig.9Hystereticcurvesofbumpfoil14whenmaximumradialdisplacementofshaftis60,70and80μmrespectively轉(zhuǎn)軸最大徑向位移分別為60,70,80μm時(shí)波箔12與波箔16的滯回曲線如圖10所示:波箔12更接近于固定端,在大回路中表現(xiàn)出更高的剛度。這是因?yàn)榉抡鏁r(shí)轉(zhuǎn)軸中心左側(cè)的波箔受固定端響應(yīng)力的影響較大;在卸載過(guò)程中,固定端的響應(yīng)力有助于箔片克服摩擦力,促進(jìn)接觸節(jié)點(diǎn)的滑動(dòng),預(yù)載荷較小時(shí)該現(xiàn)象更為明顯;隨著預(yù)載荷逐漸增大,固定端的響應(yīng)力對(duì)波箔的影響變小,這是因?yàn)檗D(zhuǎn)軸位移增加到一定程度時(shí)反作用力比摩擦力小很多,從而導(dǎo)致預(yù)載荷較小的轉(zhuǎn)軸在卸載過(guò)程中更快進(jìn)入滑動(dòng)狀態(tài)。

圖10波箔12,16的垂直位移與靜載荷的關(guān)系Fig.10Relationshipbetweenverticaldisplacementandstaticloadofbumpfoils12and16轉(zhuǎn)軸徑向位移為2.5C,摩擦因數(shù)為0.1,不同方向上波箔的耗散能量如圖11所示,在圖11的基礎(chǔ)上繼續(xù)研究不同摩擦因數(shù)下箔片的總耗散能量,結(jié)果如圖12所示:影響耗散能量的因素為靜載荷和接觸點(diǎn)的水平位移,能夠看出各波箔耗散能量與圖6中的靜載荷分布