1、旋翼槳葉固有特性計算E3摘要：本文基于彈性槳葉動力學知識，將分別對揮舞，擺振，扭轉進行模型建立：然后采用變星分離方法結合 模態疊加法，計算額定轉速狀態，槳葉的揮舞前3階、擺振前2階和扭轉1階的固有特性：采用MATLAB編程帶 入數據進行求解分析。關鍵詞：彈性槳葉分，離變量法，模態疊加一*、引言對于槳葉剖面重心無偏心的這種常規釵接式旋翼，槳葉不存在各自由度之間的耦合，或者各自 由度間的耦合很小，故在動力特性分析時，可以將揮舞、擺振、扭轉單獨分析，分析結果的精度度 誤差不大。二、模型建立及求解彈性槳葉揮舞固有特性計算1.1揮舞運動微分方程采用積分形式的牛頓法求解：半徑為1處外段上各力及力矩平衡，如

2、圖1.1取外段處一微段dp, 建立對1處的力矩平衡方程：慣性力Ft = mZ(p)dp力臂(p - r)離心力玲=?nO.2pdp力臂(Z(p) - Z(r)氣動力凡力臂(p-r)(1. 1)M(r) = 凡一 mZ)(p -r)dp - 7nO2p(Z(p) - Z(r) dp由工程梁理論，結構彎矩可以表示為：(1.2)M(r) = EIyZf,其中Z =穿,Z表示梁的撓度。 or將(1.2)式帶入(1.1)式，整理得,REIyZu +RR+ J mO2p(Z(p) -Z(r) dp + mZ(p -r)dp =rrrr)dp(1. 3)方程(1.3)兩邊對r求兩次導數得到mZ +(7),Z

3、) - Q2 J mpdp = Fz(1. 4)r.即槳葉的揮舞運動微分方程。1. 2模型求解這是一個變微分方程，通過分離變量法求解。Z(r,t)的解可以表示為= Wn(rTn(t(1. 5)Z“(r,t)表示第n階固有振動，(r)表示第n階固有振型函數，G(t)可以看做是模態對應的廣義坐 標。W(r)滿足邊界條件：(1.6)(1. 7)葉(0 =。 以(。)=0g(R) = 0Sn(R) = 0那么滿足邊界條件(1.6)、(L7)連續系統的通解，可以采用模態疊加法，表示為全部模態函數的線性組合0O(1.8)Z(r,t) = W(r)7；(t)n=l將(1.8)式帶入(1.4)式中得0O2 m

4、 叫 Tn + n=ln=l2(Eiw；y -Q2(1.9)令凡=0,同時令(1. 10)(引Q2J mpd j = mWnn則由(1.9)式可得到自由振動方程0O(1. 11)2 m叫(九 + 3京7；) = 0 n=l將(1.9)、(1.10)式等號兩邊作用算子#()WQr，根據固有振型的正交性可知,WtWndr= lfi=n 0, i W(1. 12)得到(1. 13)(1.14)Mn(Tn + nTn) = 0Kn = nMn = EIyW；2 + D2 (徂2 J式中Mn =dr第n階主質量：心為第n階主剛度。自由振動方程可寫成MnTn+KnTn = 0固有頻率ln阻Cln其中M =

5、(n = 1,2,3)為質量矩陣,K =K2n(n = 1,2,3)為剛度矩陣，C =2nMJCan-(n = 1,2,3)為特征向量。2 feEIyn2 + n2 (w；2mpdp)3 W-dr對于(1.8)式取前N個有限項作為近似解，則有(1. 16)NZ(r,t) = 叫偵)7；(t)n= 1,2,3-TVn=l(1. 17)由于實際中無法得到真實振型函數噸(r),故在計算中代以滿足分邊界條件的N個假設模態函數W(r)、W2(r)、w3(r) - wN(r)c對于本例取N=3,那么第n階模態函數I4(r)可近似表達為3吼()=+ C2Z2G)+ C3Z3G) = GmG)1=1(1. 1

6、8)其中Cizv C2n,C3n是待定系數。r e 成)=(1. 19)10 /r e310 zr ex4 zr ex!w2(r) = ( ) - ( ) + ()3以 3 Re， R r eB zr e6 商)+2(商)r e 日-6設g) = sin(，nt + +)JM (1. 15)式可以寫成(K 娣 M)C= 0的(尸)=5(1.20)(1.21)(1.22)(1.23)I mwf(r)wn(r)dr+ Q2R 1mpdp dr(1.24)槳葉揮舞固有振動歸結為求解(1-22)式廣義特征值問題，可求得3個特征值口角即固有頻率和3 個特征向量G則前三階振型函數可以表示為噸(r) =阮(

7、r) w2(r) w3(r)C(1.25)1.3計算結果將己知參數帶入上述模型，利用Matlab很容易求得槳葉揮舞前三階固有頻率，和振型函數。固有頻率為：co。 = 1.0152Q, 3陀=3.0737Q, 3陀=7.3514(1振型函數為：14(r) =0.0677182*1 + 0.00838172*(0.204082*1 - 0.0204082)A3 - 0.0090672*(0.204082*r - 0.0204082)A4+ 0.00333708*(0.204082*1 - 0.0204082)% 0.000274189*(0.204082*1 - 0.0204082)A6 0.00

8、677182lV2(r) = 1.12416*(0.204082*1 -0.0204082)勺 0.190552*r + 3.85872*(0.204082*1 - 0.0204082)A4-5.6422l*(0.204082*r-0.0204082)A5 + 1.99315*(0.204082*1- 0.0204082)A6 + 0.0190552I4(r) =0.528103*1- 42.1677*(0.204082*1- - 0.0204082)A3 + 94.1201 *(0.204082*1- - 0.0204082)A4-74.9931 *(0.204082*r-0.0204082

9、)A5+ 20.7809*(0.204082*1- - 0.0204082)A6 - 0.0528103圖1.2彈性槳葉揮舞振動前三階振型擺振固有特性結算2.1擺振運動微分方程圖2.1槳葉擺振運動受力圖同理采用積分形式牛頓法，建立擺振運動微分方程慣性力& = mx(p)dp力臂(p -r)離心力Fn = mQ2pdp 力臂- x(r)；氣動力Fx，力臂(p-r)(2. 1)Mz(r) = J (& 一 mX)(p - r) - niO.2p (x(p) ： - X(r) dp r整理得r y(2.2)mX + (EIzX) _ Q2 X1 I mpdp - QrmX = Fz2.2方程求解&r

10、,t)的解可以表示為Xn(r,t) = Wn(r)Tn(t)(2.3)Xn(r,t)表示第n階固有振動，l(r)表示第n階固有振型函數，匕。)可以看做是模態對應的廣義坐 標。Wn(r)滿足邊界條件：(e) =0(e) = 0M(R) = 0I W(R) = 0那么滿足邊界條件(2.4)、(2.5)連續系統的通解，可以采用模態疊加法，表示為全部模態函數的線性組合(2.4)(2. 5)X(r,t) = W(r)Tn(t)n=l(2.6)采用相同的處理方法，得到擺振自由振動方程為MnXn+KnTn = 0 其中Mn = m WdrKn =磁雙” =JEIZW2 + Q2 (J mpdp + O.2m

11、W dr為第n階主質量；K”為第n階主剛度。同樣可以假設”n()=ClMiO) + C2n2()= cj(r)i=l其中cln, C2n是待定系數。r ewjr)=R e10 /r e310 zr ex4 zr ex5w2(r) = ( ) - ( ) + ()3 &e， 3 R旅 、R 設W) = sin( +。),則(2.7)式可以寫成(K 一娣M)C= 0其中Af =：偵(n = 1,2,3)為質量矩陣，K =偵(”=1,2,3)為剛度矩陣，C =凱 mwf(r)wn(r)dr LM2nJLA2nJLc2n.ElyWi wn + fl2 I wtwn氏mpdp dr(2. 7)(2.8)

12、(2.9)(2. 10)(2. 11)(2. 12)(2. 13)(2. 14)(n = 1,2)為特征向量?？梢郧蟪鰯[振固有特性。2.3計算結果計算得到槳葉擺振前兩階固有頻率為= 0.1749Q, 3陀=3.64620前兩階振型函數為必(r) = 0.00523153*(0.204082*i 0.0204082)A4 0.00523153*(0.204082*r - 0.0204082)A3 0.0678462*r - 0.00156946*(0.204082*1 - 0.0204082)% + 0.0067846214)=5.14495*(0.204082*1 - 0.0204082)A3

13、 - 0.240067*r - 5.14495*(0.204082*1 - 0.0204082)A4 + 1.54349*(0.204082*1- 0.0204082)人5 + 0.0240067圖2.2槳葉扭轉前兩階振型扭轉固有特性3.1扭轉振動微分方程Q 虬考慮1處剖面受力情況如圖3.1,微段dr力矩方程為M(r) =a2Iee -Ma +I9e(3. 1)式中q2%3表示螺槳力矩，螺槳力矩是由離心力引起的，I剖面處質量微元dm上的離心力在弦向方 向會有一個分量，如圖3.2。設質量微元dm距變距軸線距離為x,那么弦向離心力分量為 (O2dmVr2 += XQ2dmQ圖3.2質量微元上的離心

14、力假設槳葉扭轉角為。，該質量微元dm離心力分量對變距軸線力矩為（址）（知2命）,那么整個剖面 進行積分即可得到螺槳力矩J 伽）印2命）=紂2section(3.2)圖3.3離心力弦向分雖對變距軸的力矩J x2dm = 0Q.2Ie section令Mq = 0,同時考慮到M（X）= GU （GJ為抗扭剛度）則由上式可得到自由扭轉振動方程 drGjey -qtIqO -iGe = o(3. 3)3.2方程求解采用分離變量法，設。（r,t） = f（r）p（t）帶入上式得扭轉振動微分方程 （G女）+ /灑（r）3 一 Q2）=。(3.4)其中她心 +”（t）= o(3. 5)f（r）滿足邊界條件G

15、次(0) = K疝(0)f(R) = 0f (R) = 1(3. 6)K為操縱線系剛度。由微分方程可以看出，槳葉旋轉和不旋轉時的振型相同，則O2 Q2 = coo敏為不旋轉時槳葉扭轉固有頻率。令Q = 0得到GJC +說我=0(3.8)該方程的解可.以表示為E(r) = A cos pr + B sin pr f(p =gT(3.9)將-r）帶入邊界條件（3.6）,得到告=5泌E(r) = cos(pR 伊尸)(3. 10)(3. 11)求解方程(3.10)得到的值，那么口音=筍，進而計算旋轉時固有頻率 le3 = J- + n2(3.12)3.3計算結果考慮操縱線系剛度=0,10000.20

16、000. 00四種情況，計算不同操縱線系剛度下的扭轉一階固有頻率, 結果如下： 瓦。=o, co = i.oooonK0 = 10000 f 3 = 3.80000Ke = 20000, co = 4.59860Kq = 8, co = 6.21S2Q一階扭轉振型函數圖如E圖3.4不同操縱線系剛度下的一階扭轉振型三、總結從以上的計算可以看到，同樣轉速卜擺振固有頻率比同階揮舞固有頻率小，主要是揮舞和擺振 時，離心力作用方式不一樣，一個是平行力系，另一個是中心力系這就導致揮舞運動比擺振運動的 離心剛度大；對于扭轉振動，離心力分量會產生螺槳力矩，但離心力剛度所占的比例也是微不足道。 從圖3.4中可以

17、看出，對于一階扭轉振型，操縱線系剛度影響明顯，當電=8時，根部扭轉角=0相 當根部固支，而其他幾種情況可以反映根部扭轉角在整個槳葉扭轉角里占比較大。在揮舞、擺振、扭轉固有特性計算中不約而同地采用了假設模態，即先要假設一個(組)振型 函數，前提是滿足邊界條件，然后再帶入方程中求固有頻率和振型函數，可見預設振型函數也是很 關鍵的，本文未涉及探究過程。以上是在理想理論模型，且相互獨立情況下計算的彈性槳葉揮舞、扭轉、擺振的固有特性，而 實際中三種運動的耦合以及旋翼與支持系統之間的耦合會帶來顯著的差別和影響。通過此門課程的學習及后續的知識補充，從無到有，熟悉和理解了旋翼主要動力學問題：振動 響應(氣彈耦

18、合)、動應力和振動、動力不穩定性。雖然實際問題復雜，掌握基礎簡化假設的理論和 處理問題的方法手段對以后深化分析有著不可或缺的作用。四、程序fiinction outl ,out2 = Flap(R,m,e,EIy,RPM)%計算槳葉揮舞前三階固有特性%Instmction%outl-固有頻率無因此化,out2-振型函數%R-旋翼半徑m,R=5;%m-槳葉線密度kg/m,m=5.5;%RPM-旋翼轉速 rpm,RPM=400;%e-餃偏置量m,e=0.1;%EIy .抗彎剛度,EIy=6.9*10M;syins r%定義符號變量OMEGA=RPM*2 *pi/60;%rad/sFN=0.5*m*

19、OMEGAA2*(RA2-iA2);% 離 心力%假設模態函數w=w 1 ,w2,w3 wl=(r-e)/(R-e);w2=( 10/3)*(r-e)/(R-e)A3 -(10/3) *(r.e)/(R Y)MH(r.e)/(R.e)A5;w3=5*(r.e)/(R-e)A4-6*(r-e)/(R-e)人 5+2*(r.e)/(R.e)人 6;w=wl w2 w3;dw=diff(w);%對1,一階導數ddv=di任(w,2);%對 r 二階導數%初始化質量陣M,剛度陣KM=zeros(3);K=zeros(3);%計算M, Kfor i=l:3for j=l:3dM=m*w(i)M(iJ)=

20、int(dM,e,R);% 質量陣dK=(EIy*ddw(i) * ddv(j 汁 dw(i) *dw(j) *FN);K(ij)=int(dKeR);% 剛度陣endend%求特征值矩陣V和特征向量矩陣CC,q=eig(K,M);for i=l:3wb(i)=r(i,i)A0.5/OMEGA;%第 i 階固有頻率Wim,e,EIz,RPM)%計算槳葉擺振前兩階固有特性%Instmction%outl-固有頻率無因此化,。ut2-振型函數%R-旋翼半徑m,R=5;%m-槳葉線密度kg/m,m=5.5;%RPM-旋翼轉速 ipm,RPM=400;%e-擺振餃偏置量m,e=0.1;%EIz .抗彎

21、剛度,EIz= 1.7 * 10抗;OMEGA=RPM*2 *pi/60;%radzssyins rFN=0.5*m*OMEGAA2*(RA2-iA2);% 離 心力%假設模態函數w=wl,w2wl=(r-e)/(R-e);w2=( 10/3)*(r-e)/(R-e)A3 -(10/3) *(r.e)/(R Y)MH(i.e)/(R.e)A5;w=wl w2;d-diff(w);ddv=difi(w2)；M=zeros(2);K=zeros(2);for i=l:2for j=l:2dM=m*w(i) *w(j);M(i,j )=mt(dM,e,R);% 質量陣dK=(EIz* ddw(i) *ddw(j)+dw *dw(j) *FN) -OMEG A A2 *m* w(i) *w(j);K(ij)=int(dK,e,R);% 剛度陣endendC,Ar=eig(K,M);%求特征值矩陣V和特征向量矩陣Cfor i=l:2wb(i)=Xi,i)A0.5/OMEGA;%第 i 階固有頻率WiImz.RPM) %Instniction%扭轉一階固有特性計算固有頻率%R-旋翼半徑m%Imy ,，Imz-剖面質量慣性矩%R