旋翼桨叶固有特性计算

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