考慮彎扭耦合運動的旋轉(zhuǎn)帶冠葉片非線性氣動彈性分析

麻岳敏，曹樹謙,3，郭虎倫

(1. 天津大學(xué) 機力學(xué)系，天津 300354； 2. 天津市非線性動力學(xué)與控制重點實驗室，天津 300354；3. 天津大學(xué) 力學(xué)國家級實驗教學(xué)示范中心，天津 300354)

葉片作為航空發(fā)動機、燃氣輪機等旋轉(zhuǎn)機械的重要部件之一，它的穩(wěn)定性直接關(guān)系到發(fā)動機能否正常工作。它在高的離心荷載、氣動荷載及振動交變荷載等荷載的惡劣工作環(huán)境中，是極易發(fā)生故障的。已有研究表明，葉片故障占發(fā)動機故障的62%以上，而葉片的損壞絕大多數(shù)是由于葉片振動產(chǎn)生的動應(yīng)力過大所致。因此，研究葉片的振動問題具有重要的意義。

從已有文獻看，對旋轉(zhuǎn)葉片的研究主要集中在將不帶葉冠的葉片簡化為旋轉(zhuǎn)懸臂梁來分析。Hodges等[1-2]建立了三自由度懸臂梁模型，研究了振動頻率，分析了顫振速度與極限環(huán)響應(yīng)，并通過實驗與理論結(jié)果進行對比。Bercin等[3-4]考慮剪切和翹曲等因素對旋轉(zhuǎn)梁彎扭耦合振動進行了分析。張偉等[5-6]通過多尺度等方法研究了受活塞氣動力預(yù)扭轉(zhuǎn)薄壁旋轉(zhuǎn)葉片的非線性動力學(xué)問題，指出隨著氣流流速的增加，系統(tǒng)呈現(xiàn)周期運動和混沌運動等復(fù)雜動力學(xué)行為。Eskandary等[7-8]用數(shù)值方法研究了Timoshenko旋轉(zhuǎn)梁的非線性振動。Saxton等[9]建立了復(fù)合材料葉片動力學(xué)方程，分析系統(tǒng)的動力學(xué)行為。姜靜等[10]建立了考慮離心力且由于偏心發(fā)生彎扭耦合的非線性動力學(xué)方程。Kashani等[11]針對Euler-Bernoulli梁分析了梁自由端受軸向載荷和力矩對系統(tǒng)彎曲扭轉(zhuǎn)振動特性和穩(wěn)定性的影響。任勇生等[12]在不考慮葉片幾何和氣動力非線性情況下，對風力機葉片的彎扭耦合線性氣彈穩(wěn)定性進行研究，并用特征值法求穩(wěn)定性。Hansen[13]對風機葉片的兩大類氣動彈性問題進行了綜述，列出了在對葉片氣動彈性穩(wěn)定性分析時所需要關(guān)注的重要問題。王丹等[14]建立彎扭耦合的葉片截面模型，指出葉片顫振過程中振動幅值隨著來流速度的增加而增大。Hosseini[15]研究了旋轉(zhuǎn)葉片在空氣熱彈性載荷下的動力學(xué)行為。Ladge等[16]通過實驗研究了渦輪機葉片在動態(tài)流作用下的穩(wěn)定性。Pourazarm等[17]考慮氣體作用下的變截面旋轉(zhuǎn)葉片的耦合振動。隨著高速旋轉(zhuǎn)機械的發(fā)展，葉冠對葉片的振動特性以及穩(wěn)定性的影響是我們不能忽視的一個問題。

李劍釗等[18]從結(jié)構(gòu)、特性分析和相關(guān)試驗介紹了汽輪機帶冠葉片。周傳月等[19-20]通過有限元方法，研究了帶冠葉片振動特性。Ansari等[21]研究了加質(zhì)量塊的旋轉(zhuǎn)梁的彎扭振動頻率，得到質(zhì)量塊質(zhì)量的增加引起頻率的增加，而末端帶質(zhì)量塊梁長度的增加降低了頻率的結(jié)論。目前，國內(nèi)外學(xué)者針對旋轉(zhuǎn)帶冠葉片受到氣動載荷和離心載荷等多種載荷共同作用時的非線性氣動彈性問題相對較少，有待進一步研究。

本文主要研究旋轉(zhuǎn)帶冠葉片的非線性氣動彈性問題?？紤]發(fā)生彎扭耦合變形，受離心力和氣動力的影響，利用Hamilton原理建立柔性葉片的非線性偏微分方程，對其無量綱化，運用Galerkin離散得到兩自由度耦合的常微分方程。通過諧波平衡法以及數(shù)值模擬研究帶冠旋轉(zhuǎn)葉片的動力學(xué)行為。

1 帶冠葉片動力學(xué)模型

將旋轉(zhuǎn)帶冠葉片簡化為如圖1所示的無預(yù)扭轉(zhuǎn)角、軸向不可伸長、末端加質(zhì)量塊的懸臂梁。考慮葉片安裝在半徑為R，角速度Ω旋轉(zhuǎn)的軸上，受到Y(jié)方向的氣動荷載。其中，梁的質(zhì)心軸與彈性軸距離為e，質(zhì)量塊質(zhì)心與梁末端的剛心重合。

為了簡化建模過程，主要假設(shè)如下：

(1) 葉冠簡化為集中質(zhì)量點，并忽略葉冠對系統(tǒng)結(jié)構(gòu)剛度的影響;

(2) 帶冠葉片為各向同性材料，本構(gòu)關(guān)系滿足胡克定律。

(1)

該點的速度可表示為

(2)

自由端質(zhì)量塊的速度表示為

(3)

系統(tǒng)的動能為

(4)

式中：L為梁的長度；ρ為梁密度；A為梁的橫截面面積；M為質(zhì)量塊質(zhì)量。

系統(tǒng)勢能由葉片變形和旋轉(zhuǎn)產(chǎn)生，表示為

(5)

式中：EI為抗彎剛度；GJ為抗扭剛度；κ(x)為曲率；Fc為葉片所受離心力；Fm為質(zhì)量塊所受離心力，其表達式分別為

(6)

Fm=M(R+L)Ω2

(7)

由材料力學(xué)可知，當懸臂梁的變形較小時，曲率可近似表達為

(8)

由于葉片看作柔性梁，有相對較大的變形，此時曲率的選取就不能近似等于式(8)，而表達為

(9)

設(shè)由彎曲變形引起的橫截面的轉(zhuǎn)角為φ(x,t)，由材料力學(xué)的知識可知，其值是一個相對較小的量，故可近似認為

(10)

(11)

由于φ比較小，將式(9)和式(11)代入式(5)，簡化為

(12)

活塞理論是一種非定常氣動理論，這種理論被廣泛應(yīng)用于超聲速和高超聲速翼面氣動彈性分析中。根據(jù)一階活塞理論氣動力模型[22]，得到二自由度系統(tǒng)受到的氣動力La和氣動力矩Ma為

(13)

(14)

系統(tǒng)虛功為

(15)

運用Hamilton原理

(16)

將式(4)、式(12)、式(15)代入式(16)得動力學(xué)方程

(17)

(18)

邊界條件

當x=0時，

w(0,t)=θ(0,t)=0,w′(0,t)=0

(19)

當x=L時,

(20)

引入無量綱量

式(17)、式(18)變?yōu)?/p>

(21)

(22)

其中，

(23)

(24)

運用Galerkin方法對無量綱方程做一階截斷，選取一階振型離散

(25)

(26)

其中，

(27)

(28)

β=1.875

(29)

將式(25)～式(29)代入式(21)和式(22)中，然后進行Galerkin正交化，整理得到兩自由度的非線性動力學(xué)方程

(30)

系數(shù)Ai和Bj(i=1,2,,7；j=1,2,,6)為系統(tǒng)參數(shù)和轉(zhuǎn)速ω(無量綱化的角速度)的函數(shù)，其表達式詳見附錄。

2 數(shù)值響應(yīng)分析

2.1 臨界顫振速度

分析彎曲位移和扭轉(zhuǎn)位移在不同轉(zhuǎn)速下的時間響應(yīng)曲線，探討系統(tǒng)穩(wěn)定性。圖2為轉(zhuǎn)速ω=5時，彎曲位移和扭轉(zhuǎn)位移時間歷程圖。由圖可以看出葉片振動收斂到零，是穩(wěn)定的。圖3和圖4分別為轉(zhuǎn)速ω=7時，彎曲位移和扭轉(zhuǎn)位移的時間歷程圖和相圖。由它們可以看出此時葉片收斂到穩(wěn)定極限環(huán)(Limit Cycle Oscillation, LCO)，發(fā)生極限環(huán)顫振。這說明在轉(zhuǎn)速5～7存在一個臨界轉(zhuǎn)速，系統(tǒng)發(fā)生Hopf分岔，由穩(wěn)定的零解分岔出穩(wěn)定的極限環(huán)。

圖2 ω=5時，無量綱彎曲位移和扭轉(zhuǎn)位移隨時間 變化的曲線Fig.2 The time history of displacement bending and torsion (dimensionless) at the speed ω=5

圖3 ω=7時，無量綱彎曲位移和扭轉(zhuǎn)位移隨時間 變化的曲線Fig.3 The time history of displacement bending and torsion (dimensionless) at the speed ω=7

圖4 ω=7時，無量綱彎曲和扭轉(zhuǎn)的相圖Fig.4 The phase diagram dimensionless bending and torsion at the speed ω=7

為了研究旋轉(zhuǎn)葉片的顫振臨界解，引入向量

(31)

則式(30)可表示為矩陣形式

(32)

式中：β為系數(shù)矩陣；C5和D5為系數(shù)，具體表達式見附錄。

這樣可以通過式(32)的特征值問題來求系統(tǒng)顫振臨界解。一般特征值都是復(fù)數(shù)，當特征值的實部都是負數(shù)時，旋轉(zhuǎn)葉片的振動是衰減振動，即系統(tǒng)是穩(wěn)定的。當特征值的實部有一個等于零時，葉片作簡諧運動，處于臨界狀態(tài)。當特征值的實部有一個是正數(shù)，則旋轉(zhuǎn)葉片振動發(fā)散。通常情況下，旋轉(zhuǎn)葉片系統(tǒng)的所有特征值實部都是負數(shù)，隨著轉(zhuǎn)速的增加，系統(tǒng)特征值的實部有一個會提前變?yōu)榱?，這時，旋轉(zhuǎn)葉片失穩(wěn)。而使特征值最大實部達到零的轉(zhuǎn)速就是旋轉(zhuǎn)葉片發(fā)生顫振的臨界轉(zhuǎn)速。圖5表示特征值最大的實部隨轉(zhuǎn)速的關(guān)系圖，可以看出轉(zhuǎn)速ω=6.149時，特征值最大實部等于零，這時解有一對純虛根，說明達到臨界條件，發(fā)生顫振。

圖5 特征值最大實部值隨轉(zhuǎn)速ω的變化關(guān)系Fig.5 The relationship between the maximum real value of the eigenvalue and the speed ω

2.2 極限環(huán)響應(yīng)

非線性因素的存在，使得氣動彈性系統(tǒng)在顫振點處發(fā)生了Hopf分岔。轉(zhuǎn)速大于顫振轉(zhuǎn)速時，非線性系統(tǒng)將會出現(xiàn)極限環(huán)，而不是線性系統(tǒng)的氣動彈性失穩(wěn)。根據(jù)Hopf分岔定理，分析極限環(huán)響應(yīng)。圖6為彎曲運動和扭轉(zhuǎn)運動的LCO幅值隨轉(zhuǎn)速的變化圖，系統(tǒng)在K點(ω=6.149)發(fā)生Hopf分岔后，出現(xiàn)了極限環(huán)，并且極限環(huán)的幅值隨轉(zhuǎn)速ω的增大而增大。

圖6 彎曲運動和扭轉(zhuǎn)運動LCO幅值(Aq和 Ap)隨轉(zhuǎn)速ω的變化關(guān)系Fig.6 The relationship between bending and torsional amplitude of LCO (Aq and Ap)and the rotational speed ω

諧波平衡法是求解非線性振動問題常用的一種近似解析法，它的求解過程歸結(jié)為代數(shù)方程組的求解。本文運用諧波平衡法研究葉片的極限環(huán)響應(yīng)?？紤]方程是自由振動方程，無外激勵，所以可以設(shè)解的模式為

q=a1cos?t+b1sin?tp=a2cos?t

(33)

將式(33)代入式(30)得

(34)

(-B1a1?2-B2a1?2+B4a2+B5b1?+B6a1)cos?t+ (-B2b1?2-B3a2?-B5a1?+B6b1)sin?t=0

(35)

略去3次諧波，由cos?t及sin?t項的系數(shù)等于零可得

(36)

(37)

-B1a2?2-B2a1?2+B4a2+B5b1?+B6a1=0

(38)

-B2b1?2-B3a2?-B5a1?+B6b1=0

(39)

由式(36)～式(39)得到a1,b1,a2,?的近似解。

對解進行分析，與數(shù)值解進行對比。如圖7所示，諧波平衡法的結(jié)果與數(shù)值結(jié)果吻合較好。兩種方法都說明在轉(zhuǎn)速ω=6.149時葉片發(fā)生顫振，出現(xiàn)極限環(huán)。而且隨著轉(zhuǎn)速的增加，極限環(huán)幅值A(chǔ)q的增長趨勢也相似。

圖7 解析解與數(shù)值解的結(jié)果對比Fig.7 Comparison of analytical solution and numerical solution

3 結(jié)構(gòu)參數(shù)對臨界顫振轉(zhuǎn)速及極限環(huán)響應(yīng)的影響

3.1 對臨界顫振轉(zhuǎn)速的影響

取l=1.5,n=1時，分析質(zhì)量塊質(zhì)量對顫振轉(zhuǎn)速的影響。從圖8可知，由于質(zhì)量塊的存在，改變了原有系統(tǒng)的穩(wěn)定性。0.1≤S≤0.2時，顫振轉(zhuǎn)速ωf急劇下降，說明當質(zhì)量塊質(zhì)量占梁質(zhì)量達到某一值時，開始影響系統(tǒng)的穩(wěn)定性。隨著S的增大導(dǎo)致顫振轉(zhuǎn)速ωf減小，系統(tǒng)容易發(fā)生顫振。

圖9為顫振轉(zhuǎn)速隨質(zhì)量塊位置的變化曲線。取l=1.5,S=0.16分析。當0.5≤n≤0.6時，發(fā)生顫振的轉(zhuǎn)速ωf有一個大的變化，質(zhì)量塊開始影響系統(tǒng)的穩(wěn)定性； 當n≥0.7時，臨界顫振轉(zhuǎn)速隨n的增大而減小，近似呈現(xiàn)線性關(guān)系。不難看出隨著質(zhì)量塊越靠近梁自由端，越容易發(fā)生顫振。

圖8 顫振轉(zhuǎn)速ωf隨質(zhì)量塊質(zhì)量S(無量綱)的變化關(guān)系Fig.8 The relationship between the flutter speed ωf and the mass S(dimensionless) of the mass

圖9 顫振轉(zhuǎn)速ωf隨質(zhì)量塊位置n(無量綱)的變化關(guān)系Fig.9 The relationship between the flutter speed ωf and the mass position n

圖10表示無質(zhì)量塊時，改變梁的長度對顫振速度的影響。一開始長度較短時，不會發(fā)生顫振。達到一臨界值后，開始發(fā)生顫振。之后隨著梁長度的增大，發(fā)生顫振的幾率就隨之增大。Resor等[23]也指出隨著葉片長度增長，顫振速度會逐漸減小。分別取S=0.16和S=0.40(相對于梁L=0.15 m時)，n=1，分析有質(zhì)量塊時，梁長度對顫振速度的影響，如圖11所示。由于質(zhì)量塊的存在，在葉片較短時會發(fā)生顫振，并且隨著梁長度的增加，顫振速度也隨之增加。但是，當葉片達到一定長度，質(zhì)量塊的影響就基本不存在了。這時和無質(zhì)量塊時的情況一樣，會隨著葉片長度的增大，而導(dǎo)致容易發(fā)生顫振。說明由于質(zhì)量塊存在，短梁會受到一定的影響。

圖10 無質(zhì)量塊時，顫振轉(zhuǎn)速ωf隨梁長度l的變化關(guān)系Fig.10 No mass, the relationship between the flutter speed ωf and the beam length l

圖11 有質(zhì)量塊時，顫振轉(zhuǎn)速ωf隨梁長度l的變化關(guān)系Fig.11 There is mass, the relationship between the flutter speed ωf and the beam length l

圖12 顫振轉(zhuǎn)速ωf隨和的變化關(guān)系Fig.12 The relationship between the flutter speed ωf

3.2 對極限環(huán)響應(yīng)的影響

圖13和圖14是分析質(zhì)量塊對極限環(huán)響應(yīng)的幅值影響。固定l和n，分別取S為0.16, 0.25, 0.50, 0.75，得到圖13所示的四條曲線，可以看出質(zhì)量塊重量越大，越容易發(fā)生Hopf分岔，極限環(huán)幅值越大。同理，固定l和S，分別取n為0.85, 0.90, 0.95, 1.00，圖14說明質(zhì)量塊位置越靠近末端，越容易發(fā)生Hopf分岔，同時極限環(huán)幅值越大。

圖13 質(zhì)量塊不同質(zhì)量下的LCO幅值A(chǔ)q 隨轉(zhuǎn)速ω的變化關(guān)系Fig.13 The relationship between the LCO amplitude Aq and the speed ω at different quality of the mass

圖14 質(zhì)量塊不同位置下的LCO幅值隨轉(zhuǎn)速ω的變化關(guān)系Fig.14 The relationship between the LCO amplitude Aq and the speed ω at different positions of the mass

4 結(jié) 論

本文建立了帶冠葉片的兩自由度耦合非線性動力學(xué)方程。由特征值理論得到了顫振臨界轉(zhuǎn)速，用諧波平衡法和數(shù)值方法對比分析了極限環(huán)響應(yīng)，結(jié)果發(fā)現(xiàn)兩者吻合的很好。研究了旋轉(zhuǎn)狀態(tài)下，質(zhì)量塊的質(zhì)量和位置、梁的長度和寬度等參數(shù)對振動響應(yīng)的影響規(guī)律。結(jié)果表明，質(zhì)量塊對顫振速度和極限環(huán)幅值有明顯的影響，質(zhì)量塊的質(zhì)量越大，位置越靠近末端，越容易發(fā)生顫振，極限環(huán)幅值越大。同時質(zhì)量塊的存在會對短梁有一定影響。本文主要針對不考慮葉片之間耦合的單個帶冠葉片的超音速、高超音速非線性氣動彈性的穩(wěn)定性分析，對多個葉片之間的耦合問題將在后續(xù)研究中給與考慮。