達(dá)琦,袁奇,李浦

(1.西安交通大學(xué)能源與動力工程學(xué)院,710049,西安;2.西安交通大學(xué)陜西省葉輪機(jī)械及動力裝備工程實驗室,710049,西安)

拉桿轉(zhuǎn)子常用于重型燃?xì)廨啓C(jī)和航空發(fā)動機(jī),通過在中心或周向拉桿螺栓上施加預(yù)緊力來保證多級輪盤之間的相對靜止并使輪盤段具有一定剛度。然而,即使受到較小的彎曲外載荷,轉(zhuǎn)子輪盤接觸面仍會部分區(qū)域脫開[1],導(dǎo)致接觸面積減小進(jìn)而造成轉(zhuǎn)子整體剛度下降。文獻(xiàn)[2]對某航空發(fā)動機(jī)的鼓筒式轉(zhuǎn)子接觸段剛度進(jìn)行研究,結(jié)果表明輪盤之間的相對轉(zhuǎn)角隨彎矩線性增加呈現(xiàn)出非線性的變化規(guī)律。在建立具有螺栓連接結(jié)構(gòu)的航空發(fā)動機(jī)轉(zhuǎn)子力學(xué)模型時,常采用角剛度分段線性模型來近似連接部位的剛度變化規(guī)律[3]。與連續(xù)轉(zhuǎn)子相比,角剛度分段線性模型的動力學(xué)響應(yīng)與線性模型的動力學(xué)響應(yīng)存在較大差異。為了解釋組合轉(zhuǎn)子升速和降速中存在的振幅突跳現(xiàn)象,有學(xué)者綜合考慮輪盤接觸面和拉桿剛度的共同作用,建立具有三次方剛度的轉(zhuǎn)子模型[4-5]。拉桿轉(zhuǎn)子的特殊結(jié)構(gòu)使得連接段剛度具有非線性的變化規(guī)律,且組合式轉(zhuǎn)子多個輪盤之間的接觸模式可能不同,因而其動力學(xué)響應(yīng)十分復(fù)雜[6];在特定條件下,幅頻響應(yīng)中出現(xiàn)的線性系統(tǒng)所不具有的雙穩(wěn)態(tài)現(xiàn)象使得轉(zhuǎn)子故障頻發(fā),造成了巨大的經(jīng)濟(jì)損失。對于非線性轉(zhuǎn)子系統(tǒng)出現(xiàn)的振幅突跳現(xiàn)象,國內(nèi)外學(xué)者從不同角度進(jìn)行研究[4,7-10],將其視為轉(zhuǎn)子系統(tǒng)的一種故障,并盡量設(shè)法避免。目前國內(nèi)對于拉桿轉(zhuǎn)子非線性振動特性研究較少,以往的文獻(xiàn)均采用剛度硬化模型,即整體剛度隨轉(zhuǎn)速增加而增加,對于幅頻響應(yīng)中出現(xiàn)的振幅突跳現(xiàn)象無一致結(jié)論,且現(xiàn)有分析模型缺乏依據(jù)。

針對以上問題,本文首先利用三維非線性有限元方法計算得到拉桿與接觸面等效剛度隨輪盤橫向間距的變化規(guī)律;然后分別采用雙線性模型和三次方模型擬合非線性剛度曲線,建立兩種拉桿轉(zhuǎn)子力學(xué)模型;最后通過數(shù)值仿真分析,對兩組模型幅頻響應(yīng)中出現(xiàn)的非線性振動現(xiàn)象進(jìn)行了對比研究,進(jìn)而分析不同系統(tǒng)參數(shù)對轉(zhuǎn)子振動特性的影響。

1 拉桿轉(zhuǎn)子橫向剛度分析

圖1為周向拉桿組合式轉(zhuǎn)子模型,該模型由軸頭、輪轂、輪盤、拉桿和接觸面組成,l3、l4分別為輪轂和輪盤的軸向長度,2l1表示輪轂和輪盤總長度,l2為軸頭長,模型總長為2l。將輪盤接觸面和拉桿的剛度一并考慮,Ks、Kh、Kds、Kr分別表示軸頭、輪轂、輪盤和拉桿的抗彎剛度。

(a)拉桿轉(zhuǎn)子示意圖

(b)轉(zhuǎn)子彎曲剛度模型圖1 周向拉桿組合式轉(zhuǎn)子

轉(zhuǎn)子輪盤段的剛度可表示為

(1)

利用彎矩-轉(zhuǎn)角關(guān)系和有限元軟件計算輪盤段的抗彎剛度

(2)

Kfij表示在j方向上引起轉(zhuǎn)子在i方向上產(chǎn)生單位轉(zhuǎn)角所施加的彎矩。

為簡化分析,認(rèn)為轉(zhuǎn)子各向同性,并不計交叉方向上的剛度,彎曲剛度可以表示為

(3)

(4)

式中δij為克羅內(nèi)克函數(shù)。利用式(3)可以在三維有限元軟件ANSYS中求出輪盤段的抗彎剛度Kc,然后將Kc代入式(1)計算得到Kr。

用力矩-面積法求輪盤段的橫向剛度

(5)

利用力矩面積法分別得到輪轂和輪盤的橫向剛度Khl、Kdsl,然后根據(jù)Kcl、Khl、Kdsl和Krl之間的關(guān)系,即可得到拉桿和接觸面的等效橫向剛度

(6)

由于輪盤相對轉(zhuǎn)角很小,認(rèn)為sinφ≈φ,因此兩輪盤間距X近似表示為φl,其中l(wèi)為接觸面到軸頭的距離,橫向位移與轉(zhuǎn)角關(guān)系示意圖如圖2所示。

圖2 轉(zhuǎn)子輪盤橫向位移與轉(zhuǎn)角關(guān)系示意圖

等效剛度對應(yīng)的橫向回復(fù)力可表示為

Frl=KrlX

(7)

利用有限元軟件ANSYS建立上述三維有限元模型,如圖3所示,模型的參數(shù)取值如表1所示。

表1 模型轉(zhuǎn)子物理參數(shù)

圖3 拉桿轉(zhuǎn)子三維有限元模型

參考文獻(xiàn)[11]的方法,利用式(3)計算得到輪盤相對轉(zhuǎn)角φ隨彎矩M變化關(guān)系,如圖4所示,并根據(jù)式(5)(6)計算得到等效橫向剛度值|Krl|隨兩輪盤間距X的變化曲線,如圖5所示。

圖4 輪盤間相對轉(zhuǎn)角隨彎矩變化關(guān)系

圖5 等效橫向剛度隨輪盤間距變化曲線

2 轉(zhuǎn)子動力學(xué)方程建模

將兩輪盤分別簡化為質(zhì)量為m1、m2的剛性圓盤,兩盤質(zhì)量偏心量取相同值且均用e表示,偏心量矢量夾角為θ。將拉桿和接觸面等效為不計質(zhì)量的抗彎彈簧,回復(fù)力用F(q)表示,其中q為輪盤位移矢量。輪盤段兩端分別與不計質(zhì)量的彈性軸相連,兩盤橫向剛度分別為k1、k2。左、右端軸和抗彎彈簧的阻尼系數(shù)分別表示為c1、c2、c3。轉(zhuǎn)子系統(tǒng)模型如圖6所示。

圖6 轉(zhuǎn)子系統(tǒng)模型示意圖

模型中單個輪盤的位移矢量表示為

qi=[xiyi]T

(8)

式中x、y為單個輪盤的自由度。

根據(jù)達(dá)朗貝爾原理建立系統(tǒng)運動方程

(9)

(10)

(11)

(12)

式中:Ω為轉(zhuǎn)子角速度;F(x1-x2)為非線性回復(fù)力。記X=[x1x2y1y2]T,式(9)～(12)的矩陣形式為

MX″+CX′+KX=f(t,X)

(13)

式中

3 三次方非線性剛度模型

由三維有限元計算結(jié)果如圖4所示,可知隨著接觸段載荷的增加,拉桿和接觸面等效橫向剛度呈非線性變化,其變化趨勢可以用三次多項式函數(shù)擬合,將回復(fù)力表示為

F(x)=k3x-k4x3

(14)

式中:k3為接觸段的線性剛度;k4為輪盤接觸面積減小造成兩輪盤橫向剛度削弱后的剛度。

3.1 模型動力學(xué)方程

動力學(xué)方程在x和y方向無耦合,根據(jù)控制方程的對稱性,只研究x方向的橫向振動。

令τ=Ωt,x3=x1-x2,k1=β2k2=β3k3=β4k4,將式(9)～(12)改寫為無量綱化形式

(15)

利用諧波平衡法近似求解式(13),動態(tài)響應(yīng)形式解為

(16)

(17)

利用迭代算法求解式(17),可得方程組的近似解。

3.2 非線性代數(shù)方程組優(yōu)化求解算法

3.2.1Broyden迭代算法Broyden算法具有超線性收斂速度,在一定條件下收斂速度快[12-13]。迭代過程可表示為

r(k+1)=r(k)-Ak-1f(r(k))

(18)

Aksk=yk

(19)

(20)

當(dāng)殘差小于誤差限時求解收斂。

(21)

圖7 模型轉(zhuǎn)子幅頻曲線計算流程圖

3.3 計算結(jié)果與分析

計算時使用系統(tǒng)參數(shù)如表2所示。根據(jù)計算結(jié)果繪制拉桿轉(zhuǎn)子在升速和降速過程中的無量綱幅頻響應(yīng)曲線,進(jìn)而改變系統(tǒng)參數(shù),分析不同系統(tǒng)參數(shù)對響應(yīng)曲線的影響。

表2 拉桿轉(zhuǎn)子非線性系統(tǒng)參數(shù)取值

3.3.1 轉(zhuǎn)子升速和降速曲線 取表2中第1組參數(shù),其余參數(shù)α=0.75,β2=1,β3=1,β4=5×105m-2,繪制輪盤質(zhì)心x方向的相對幅值r隨轉(zhuǎn)子相對轉(zhuǎn)速λ變化曲線,如圖8所示,圖中rij表示第i個轉(zhuǎn)盤的j階諧波響應(yīng)。由圖8可知:無量綱幅頻曲線具有一階、三階諧波,二階諧波幅值過小;兩輪盤無量綱幅頻響應(yīng)曲線隨無量綱轉(zhuǎn)速變化趨勢相同,但幅值不同。

(a)輪盤1升速和降速一階諧波曲線

(b)輪盤1升速和降速三階諧波曲線

(c)輪盤2升速和降速一階諧波曲線

(d)輪盤2升速和降速三階諧波曲線圖8 拉桿轉(zhuǎn)子輪盤升速/降速幅頻特性曲線

比較β4=0(即無非線性剛度項)與β4=5×105m-2時輪盤1的一階諧波曲線,如圖9所示,可知非線性剛度使得轉(zhuǎn)子二階共振轉(zhuǎn)速發(fā)生了明顯的變化,但一階共振轉(zhuǎn)速變化相對不大。

圖9 拉桿轉(zhuǎn)子輪盤1一階諧波對比曲線

圖8所示的無量綱幅頻曲線在λ<1時,三階諧波曲線的兩個共振峰分別對應(yīng)一階、二階固有頻率的1/3次諧波。次諧波的出現(xiàn)表明三次方剛度模型受到非線性剛度項的影響表現(xiàn)出與線性模型不同的振動響應(yīng)特性。一階、三階諧波曲線均可觀察到振幅突跳現(xiàn)象,但升速、降速曲線發(fā)生突跳的位置不同。這是因為動力學(xué)方程中存在非線性項,使同一轉(zhuǎn)速下的方程存在多解。

根據(jù)解的穩(wěn)定性判斷方法[15],將解分為穩(wěn)定解和不穩(wěn)定解。不穩(wěn)定解不能長期存在,因此隨著系統(tǒng)參數(shù)的變化,無量綱幅值作為微分方程組的解從單解區(qū)域進(jìn)入不穩(wěn)定區(qū)域(同樣也是多解區(qū)域)時,會發(fā)生從一個穩(wěn)定解突跳至另一個穩(wěn)定解的行為。由于系統(tǒng)控制方程只在特定的頻率范圍存在多解區(qū)域,在模擬升速或者降速時,由于進(jìn)入多解區(qū)域時對應(yīng)的頻率不同,使得振幅突跳位置不同。

3.3.2 結(jié)構(gòu)參數(shù)對振動特性的影響 為分析系統(tǒng)響應(yīng)出現(xiàn)的振幅突跳現(xiàn)象,有必要研究系統(tǒng)參數(shù)對解響應(yīng)曲線的影響。取表2中參數(shù),分別比較偏心量e、不平衡量夾角θ和無量綱阻尼系數(shù)ζ在不同取值下所對應(yīng)的升速幅頻響應(yīng)曲線,其余參數(shù)β2=β3=1,β4=5×105m-2?？紤]到左右輪盤振動規(guī)律類似,非線性現(xiàn)象主要發(fā)生在系統(tǒng)二階固有頻率處,且二階、三階諧波相比一階諧波幅值過小,因此只繪制輪盤1的一階諧波在其第2階共振頻率附近的圖像。

模型轉(zhuǎn)子不同參數(shù)對響應(yīng)曲線影響對比如圖10所示,可知增加偏心量e或兩輪盤間不平衡量夾角θ,減小無量綱阻尼系數(shù)ζ都會使得第2階共振峰向左方傾斜,并使多解存在的范圍擴(kuò)大。

(a)偏心量的影響

(b)兩盤偏心距夾角的影響

(c)無量綱阻尼的影響圖10 模型轉(zhuǎn)子不同參數(shù)對響應(yīng)曲線的影響

不同質(zhì)量比時輪盤1的幅頻響應(yīng)曲線如圖11所示,計算參數(shù)取值見表3,其余參數(shù)θ=90°,β2=1,β3=1。由圖11中可知,增加質(zhì)量比會使轉(zhuǎn)子系統(tǒng)第1階、第2階固有頻率減小。輪盤1第2階無量綱共振幅值隨質(zhì)量比增加而增大。輪盤2對應(yīng)無量綱幅頻曲線由于篇幅所限未展示,但其共振時無量綱幅值變化規(guī)律與輪盤1相反,第2階無量綱共振幅值隨質(zhì)量比的增加而減小。

圖11 不同質(zhì)量比時輪盤1的幅頻響應(yīng)曲線

組號αe/μmζi/10-3β4/105 m-210.45105.7520.60105.7530.70105.7540.90105.75

圖12 非線性剛度占比對輪盤1幅頻響應(yīng)曲線的影響

共振峰的偏移使系統(tǒng)存在多解的頻率范圍擴(kuò)大,進(jìn)而使振幅不穩(wěn)定區(qū)域擴(kuò)大。減小輪盤偏心量,兩盤偏心量夾角和非線性剛度占比,或者增加無量綱阻尼均能縮減雙穩(wěn)態(tài)的存在區(qū)域;而改變兩輪盤之間的質(zhì)量比,則能夠改變振幅突跳發(fā)生的位置。

表4 對比ξ時轉(zhuǎn)子系統(tǒng)的參數(shù)取值

3.3.3 結(jié)構(gòu)參數(shù)影響分析 結(jié)構(gòu)參數(shù)對系統(tǒng)響應(yīng)的影響結(jié)果可從能量和非線性剛度影響作用強(qiáng)弱的角度進(jìn)行分析。

(1)由式(15)可知,增加偏心量e相當(dāng)于增加非線性剛度占比ξ,同時具有減小拉桿剛度的作用。隨著偏心量增加,非線性剛度影響增大但總剛度減小,輪盤1的響應(yīng)曲線不穩(wěn)定區(qū)域即解的雙穩(wěn)態(tài)區(qū)域變大,第2階固有頻率減小。

(2)增加不平衡量夾角會使兩輪盤的位移差即(r1-r2)增大,進(jìn)而使得非線性剛度k4(r1-r2)3增大,使輪盤1二階固有頻率處解存在的雙穩(wěn)態(tài)區(qū)域擴(kuò)大。

(3)在外部輸入能量不變的情況下,增加無量綱阻尼系數(shù)ζ使系統(tǒng)在一個周期內(nèi)被阻尼消耗的能量增加,因而系統(tǒng)響應(yīng)幅值減小。

(4)由式(15)可知,增加質(zhì)量比α相當(dāng)于增加輪盤1的質(zhì)量比,解曲線的振幅增加而各階固有頻率減小。

4 雙線性剛度模型

隨著輪盤相對位移的變化,圖4中關(guān)系曲線可分成斜率不同的兩個部分,變化規(guī)律可近似用兩段斜率不同的直線表示。假設(shè)拉桿和接觸面的等效剛度隨拉桿的伸長呈線性變化,則對應(yīng)的回復(fù)力為

F(s)=κs

(22)

式中κ、s分別為拉桿、接觸面的等效剛度和拉桿的伸長量。

固定輪盤1,接觸段的橫向剛度為

(23)

式中:sx為輪盤2的位移在x方向上的投影;F(s)為拉桿受到的拉力。

轉(zhuǎn)子受到載荷發(fā)生彎曲,拉桿變形的形態(tài)如圖13所示,用一條折線近似描述拉桿的變形。拉桿伸長在x方向上的投影為

(24)

式中:θ1、θ2分別為彎折前、后拉桿與z軸的夾角;scr為發(fā)生彎折時拉桿的伸長量。

(a)未受載荷時接觸部位

(b)受外載荷時拉桿彎曲示意圖圖13 拉桿受載荷變形示意圖

將式(22)(24)代入式(23)可得

(25)

用x代替sx表示輪盤的橫向位移,由式(24)可知等效橫向剛度呈雙線性,結(jié)果如圖14所示。與有限元剛度計算結(jié)果對比可知,兩個模型等效回復(fù)力曲線相符合,雙線性模型可近似反映剛度的變化趨勢,圖14中兩斜率關(guān)系為kl2=0.8kl1。

圖14 雙線性剛度及與有限元模擬剛度對比

4.1 雙線性剛度模型動力學(xué)方程

設(shè)xcr為拉桿發(fā)生臨界彎折時在x方向上的伸長量,則

xcr=scrsinθ1

(26)

記單個圓盤各方向的橫向剛度為

(27)

將式(23)(24)代入式(14),得到雙線性模型的控制方程組

(28)

定義函數(shù)

(29)

式中:rcr=xcr/e為拉桿彎折時的無量綱閾值;βli=kli/k1(i=1,2)為等效剛度與輪盤1剛度的比值;γ=βl2/βl1為彎折前、后剛度的比值。

4.2 轉(zhuǎn)子穩(wěn)態(tài)響應(yīng)行為分析

取不同轉(zhuǎn)速Ω,并使用4階龍格-庫塔方法求解式(28),記錄穩(wěn)態(tài)響應(yīng)幅值。取無量綱阻尼系數(shù)ζ1=ζ2=ζ3=0.01,α=0.75,γ=0.7,rcr=10,β2=βl1=1,計算得到無量綱幅頻曲線,如圖15所示,可知雙線性模型無量綱幅頻曲線中同樣出現(xiàn)了振幅突跳現(xiàn)象。輪盤2的無量綱幅頻響應(yīng)曲線與輪盤1的類似。

圖15 轉(zhuǎn)子雙線性模型輪盤1幅頻響應(yīng)曲線

無量綱幅頻響應(yīng)曲線中出現(xiàn)的振幅突跳現(xiàn)象由剛度突變造成。當(dāng)r1與r2的差值接近rcr并且不斷增大時,響應(yīng)曲線先是經(jīng)歷了一個振蕩區(qū)域,在此范圍內(nèi),無量綱振幅變化無明顯規(guī)律;然后無量綱振幅突然增加,通過第二節(jié)臨界轉(zhuǎn)速后,隨著轉(zhuǎn)速增大,無量綱幅值表現(xiàn)出規(guī)律性的減小,如圖16所示,λ取值為1.76～2.20。兩輪盤位移差值隨轉(zhuǎn)速變化圖如圖17所示,可知當(dāng)無量綱差值為9.17接近xcr時,響應(yīng)曲線開始變得不穩(wěn)定。

(a)模型1雙穩(wěn)態(tài)區(qū)域 (b)模型2不穩(wěn)定區(qū)域圖16 兩模型不穩(wěn)定區(qū)域?qū)Ρ?/p>

圖17 兩輪盤位移差值隨轉(zhuǎn)速變化圖

4.3 三次方剛度模型與雙線性剛度模型對比

采用三次方剛度模型和雙線性剛度模型對有限元計算結(jié)果進(jìn)行近似并分析其非線性響應(yīng),計算結(jié)果均表現(xiàn)出振幅突跳現(xiàn)象,且突跳位置基本相同。兩模型差別和聯(lián)系在于:三次方剛度模型的升速和降速幅頻曲線振幅突跳發(fā)生的位置不同,而兩個突跳之間的頻段是振幅的雙穩(wěn)態(tài)區(qū)域。對比兩個模型的計算結(jié)果得出結(jié)論:轉(zhuǎn)子系統(tǒng)剛度非線性變化是影響振幅突跳和致使系統(tǒng)不穩(wěn)定的原因。

5 結(jié) 論

本文利用三次方剛度模型和雙線性剛度模型近似拉桿和接觸面等效橫向剛度,分別建立了轉(zhuǎn)子動力學(xué)模型,并對兩種模型的計算結(jié)果進(jìn)行對比,得到如下結(jié)論。

(1)兩種近似模型均能模擬拉桿轉(zhuǎn)子的非線性動力學(xué)特性,并反映出組合式轉(zhuǎn)子在工作中發(fā)生的振幅突跳現(xiàn)象。

(2)三次方剛度模型的幅頻曲線,振幅突跳原因是動力學(xué)方程在一定頻段存在多解,且存在諧波響應(yīng);但其升速和降速幅頻曲線發(fā)生振幅突跳的位置不同,兩突跳之間的區(qū)域即雙穩(wěn)態(tài)區(qū)域。雙線性模型中,振幅發(fā)生突跳在于剛度發(fā)生了突變,并在此之前經(jīng)歷一個不穩(wěn)定區(qū)域,該區(qū)域?qū)?yīng)前一模型的雙穩(wěn)態(tài)區(qū)域。

(3)減小輪盤偏心量e、兩盤偏心量間的夾角θ,或者增加無量綱阻尼系數(shù)ζ,均可使解響應(yīng)曲線的不穩(wěn)定區(qū)域減小,而改變輪盤質(zhì)量比α,則能夠改變此不穩(wěn)定區(qū)域的位置。通過調(diào)整系統(tǒng)參數(shù),能有效避免轉(zhuǎn)子系統(tǒng)在不穩(wěn)定的頻率范圍內(nèi)工作。