故障转子系统的非线性振动分析与诊断方法附录A matlab程序.doc

《故障转子系统的非线性振动分析与诊断方法附录A matlab程序.doc》由会员分享，可在线阅读，更多相关《故障转子系统的非线性振动分析与诊断方法附录A matlab程序.doc（23页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、A.1 传递距阵法分析程序%该程序使用Riccati传递距阵法计算转子系统的临界转速及振型%本函数中均采用国际单位制% 第一步：设置初始条件（调用函数shaft_parameters）%初始值设置包括：轴段数N，搜索次数M%输入轴段参数：内径d,外径D，轴段长度l，支撑刚度K，单元质量mm，极转动惯量JppN,M,d,D,l,K,mm,Jpp=shaft_parameters;% 第二步：计算单元的5个特征值（调用函数shaft_pra_cal）%单元的5个特征值：%m_k：:质量%Jp_k：极转动惯量%Jd_k：直径转动惯量%EI：弹性模量与截面对中性轴的惯性矩的乘积%rr：剪切影响系数m_

2、k,Jp_k,EI,rr=shaft_pra_cal(N,D,d,l,Jpp,mm);% 第三步：计算剩余量（调用函数surplus_calculate）,并绘制剩余量图%剩余量：D1for i=1:1:M ptx(i)=0; pty(i)=0;endfor ii=1:1:M wi=ii/1*2+50; D1,SS,Sn=surplus_calculate(N,wi,K,m_k,Jp_k,JD_k,l,EI,rr); D1;pty(ii)=D1;ptx(ii)=w1endylabel(剩余量)；plot(ptx,pty)xlabel(角速度red/s);grid on% 第四步：用二分法求固有

3、频率及振型图%固有频率：Critical_speedwi=50;for i=1:1:4order=i D1,SS,Sn=surplus_calculate(N,wi,k,m_k,Jp_k,Jd_k,l,EI,rr); Step=1;D2=D1;kkk=1;while kkk0 wi=wi+step; D2=D1; D1,SS,Sn=surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr); end if D1*D20 wi=wi-step; step=step/2; wi=wi+step; D1,SS,Sn =surplus_calculate(N,wi,

4、K,m_k,Jp_k,Jd_k,l,EI,rr); EndD1;Wi;If atepFen*60 for i=1:1:N+1 x(i,1)=yn(i*4-3,n)*le6; y(i,1)=yn(i*4-2,n)*le6; sitax(I,1)=yn(i+4-0,n)/pi*180 end u=F(loca*4-4+1,1); ut1=ut1;t u; xt1=xt1;t x; yt1=yt1;t y; endendrub_signsave xt1 ascii;save yt1 ascii;save ut1 ascii;%该程序主要进行仿真条件设置Function rub_sign loca l

5、oc_rob Famp wi r Famp1 fai=initial_conditions%需要设置的初始条件有：%rub_sign: 碰摩标志，若rub_sign=0,说明系统无碰摩故障；否则 rub_sign=1%loca: 不平衡质量的位置%loc_rub: 碰摩位置%Famp: 不平衡质量的大小 单位为：g%wi: 转速 单位为：rad% r： 偏心半径 单位为:mm%Famp1： 离心力的大小 单位为:kg,m%fai: 不平衡量的初始相位 radrub_sign=0;loca=6;loc_rub=8;Famp=1;wi=3000/60*2*pi;r=30Famp1=Famp(1)/

6、1000*wi2*r/1000;fai=30 30/180*pi %ro%该程序对Jeffcott转子系统进行参数设置 function N density Ef L R R0 miu=rotor_parameters%N: 划分的轴段数%density: 轴的密度 单位为：kg/m3%Ef: 轴的弹性模量 单位为:Pa%L 每个轴段的长度 单位为:m%R 每个轴段的外半径 单位为:m%R0: 每个轴段的内半径 单位为:m%miu: 每个轴段的单元质量 kg/mN=11;Density=7850;Ef=2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1*le

7、ll;L=90.5 90.5 80.5 62.5 30 20 22.5 62.5 90.5 90.5 90.5/1000;R=20 20 20 20 20 90 20 20 20 20 20/2000;R0=0 0 0 0 0 0 0 0 0 0 0/2000;for i=1:1;N miu(i)=density*pi*(R(i)2-R0(i)2)end%Mst%该程序设置单元距阵Function Mst Msr Ks Ge= Mst_Msr_Ks_Ge(N,density,R,R0,L,Ef,miu)%Mst: 移动单元质量距阵%Msr: 转动单元质量距阵%Ks: 刚度单元距阵%Ge: 陀螺

8、力距单元距阵NN0=N;NN1=NN0+1NN2=NN1+1for i=1:1:NN0 Mst(1,1,i)=156;Mst(2,1,i)=0; Mst(2,2,i)=156;Mst(3,1,i)=0; Mst(3,2,i)=-22*L(i); Mst(3,3,i)= 4*L(i)2; Mst(4,1,i)22*L(i); Mst(4,2,i)=0; Mst(4,3,i)=0;Mst(4,4,i)=4*L(i)2; Mst(5,1,i)=54; Mst(5,2,i)=0Mst(5,3,i)=0; Mst(5,4,i)=13*L(i); Mst(6,1,i)=0;Mst(6,2,i)=54; M

9、st(6,3,i)=13*L(i); Mst(6,4,i)=0;Mst(7,1,i)=0; Mst(7,2,i)=13*L(i); Mst(7,3,i)=-3*L(i)2;Mst(7.4.i)=0; Mst(8,1,i)=-13*L(i); Mst(8,2,i)=0;Mst(8,3,i)=0; Mst(8,4,i)=-3*L(i)2; Mst(5,5,i)=156;Mst(6,5.i)=0; Mst(6,6,i)=156; Mst(7,5,i)=0; Mst(7,6,i)=22*L(i); Mst(7,7,i)=4*L(i)2;Mst(8,5,i)=-22*L(i); Mst(8,6,i)=0

10、 Mst(8,7,i)=0Mst(8,8,i)=4*L(i)2;endfor i=1:1:NN0Msr(1,1,i)=36;Msr(2,1,i)=0; Msr(2,2,i)=36;Msr(3,1,i)=0 Msr(3,2,i)=-3*L(i); Msr(3,3,i)=4*L(i)2;Msr(4,1,i)=3*L(i); Msr(4,2,i)=0; Msr(4,3,i)=0;Msr(4,4,i)=4*L(i)2; Msr(5,1,i)=36; Msr(5,2,i)=0;Msr(5,3,i)=0; Msr(5,4,i)=-3*L(i); Msr(6,1,i)=0;Msr(6,2,i)=-36; M

11、sr(6,3,i)= 3*L(i); Msr(6,4,i)=0;Msr(7,1,i)=0; Msr(7,2,i)=-3*L(i); Msr(7,3,i)=-L(i)2;Msr(7,4,i)=0; Msr(8,1,i)=3*L(i); Msr(8,2,i)=0;Msr(8,3,i)=0; Msr(8,4,i)=-L(i)2; Msr(5,5,i)=36;Msr(6,5,i)=0; Msr(6,6,i)=36; Msr(7,5,i)=0;Msr(7,6,i)=3*L(i); Msr(7,7,i)=4*L(i)2; Msr(8,5,i)=-3*L(i);Msr(8,6,i)=0; Msr(8,7,i

12、)=0; Msr(8,8,i)=4*L(i)2;endfor i=1:1:NN0Ge(1,1,i)=0;Ge(2,1,i)=36; Ge(2,2,i)=0;Ge(3,1,i)=-3*l(i); Ge(3,2,i)=0; Ge(3,3,i)=0;Ge(4,1,i)=0; Ge(4,2,i)=-3*L(i); Ge(4,3,i)=4*L(i)2;Ge(4,4,i)=0; Ge(5,1,i)=0; Ge(5,2,i)=36;Ge(5,3,i)=-3*L(i); Ge(5,4,i)=0; Ge(6,1,i)=-36;Ge(6,2,i)=0; Ge(6,3,i)=0; Ge(6,4,i)=-3*L(i0

13、);Ge(7,1,i)=-3*L(i); Ge(7,2,i)=0; Ge(7,3,i)=0;Ge(7,4,i)=L(i)62; Ge(8,1,i)=0; Ge(8,2,i)=-3*L(i);Ge(8,3,i)=-L(i)2; Ge(8,4,i)=0; Ge(5,5,i)=0;Ge(6,5,i)=36; Ge(6,6,i)=0; Ge(7,5,i)=3*L(i);Ge(7,6,i)=0; Ge(7,7,i)=0; Ge(8,5,i)=0;Ge(8,6,i)=3*L(i); Ge(8,7,i)=4*L(i)2; Ge(8,8,i)=0;endfor i=1:1:NN0Ks(1,1,i)=12;Ks

14、(2,1,i)=0; Ks(2,2,i)=12;Ks(3,1,i)=0; Ks(3,2,i)=-6*L(i); Ks(3,3,i)=4*L(i)2;Ks(4,1,i)=6*L(i); Ks(4,2,i)=0; Ks(4,3,i)=0;Ks(4,4,i)=4*L(i)2; Ks(5,1,i)=-12; Ks(5,2,i)=0;Ks(5,3,i)=0; Ks(5,4,i)=-6*L(i); Ks(6,1,i)=0;Ks(6,2,i)=-12; Ks(6,3,i)=6*L(i); Ks(6.4.i)=0;Ks(7,1,i)=0; Ks(7,2,i)=-6*L(i); Ks(7,3,i)=2*L(i)

15、2;Ks(7,4,i)=0; Ks(8.1,i)=6*L(i); Ks(8,2,i)=0;Ks(8,3,i)=0; Ks(8,4,i)=2*L(i)62; Ks(5,5,i)=12;Ks(6,5,i)=0; Ks(6,6,i)=12; Ks(7,5,i)=0;Ks(7,6,i)=6*L(i); Ks(7,7,i)=4*L(i)2; Ks(8,5,i)=-6*L(i);Ks(8,6,i)=0; Ks(8,7,i)=0; Ks(8,8,i)=4*L(i)2;endfor i=1:1:NN0 for j=1:1:8 for k=1:1:8 EI=Ef(i)*pi*(R(i)4-R0(i)4-)/4;

16、 Mst(j,k,i)=Mst(j,k,i)*miu(i)*L(i)/420; Msr(j,k,i)=Msr(j,k,i)*miu(i)*R(i)2/120/L(i); Ge(j,k,i)=-Ger(j,k,i)*2*miu(i)*R(i)2/120/L(i); Ks(j,k,i)=Ks(j,k,i)*EI/L(i)3; end endendfor i=1:1:NN0 for j=1:1:8 for k=j:1:8 Mst(j,k,i)=Mst(k,j,i); Msr(j,k,i)=Msr(k,j,i); Ks(j,k,i)=Ks(k,j,i); Ge(j,k,i)=-Ge(k,j,i); e

17、ndendend%该程序进行距阵组集function M G K=M_G_K(N,Ef,R,R0,Mst,Msr,Ge,Ks,miu,L)% M: 总的质量距阵% K：总的刚度距阵% G：总的陀螺力距距阵NN0=Nfor i=1:1:NN0 for j=1:1:8 for K=1:1:8 Ms(j,k,i)=Mst(j,k,i)+Msr(j,k,i); Ks(j,k,i)=Ks(j,k,i); Ge(j,k,i)=-Ge(j,k,i); end endengfor i=1:1:N for j=1:1:8 for K=1:1:8 M(i*4+j-4,i*4+k-4)=Ms(j,k,i); G(i

18、*4+j-4,i*4+k-4)=Ge(j,k,i); K(i*4+j-4,i*4+k-4)= K (j,k,i); end endendfor i=2:1:Nfor j=1:1;4 for k=1:1:4 M(i*4+j-4,i*4+k-4)= M(i*4+j-4,i*4+k-4)+Ms(j+4,k+4,i-1); G(i*4+j-4,i*4+k-4)= G(i*4+j-4,i*4+k-4)+Ge(j+4,k+4,i-1); K(i*4+j-4,i*4+k-4)= K(i*4+j-4,i*4+k-4)+Ks(j+4,k+4,i-1); end endendKzx(1)=7e7;Kzy(1)=7

19、e7;Kzx(12)=7e7;Kzy(12)=7e7;for i=1:1:N+1 K(i*4+1-4,i*4+1-4)= K(i*4+1-4,i*4+1-4)+Kzx(i); K(i*4+2-4,i*4+2-4)= K(i*4+2-4,i*4+2-4)+Kzy(i);end%该程序将组尼加入总体距阵function K C D Ax=K_D(N,K,M,G)Kzx(1)=7e7;Kzy(1)=7e7;Kzx(12)=7e7;Kzy(12)=7e7;for i=1:1:N+1 K(i*4+1-4,i*4+1-4)= K(i*4+1-4,i*4+1-4)+Kzx(i); K(i*4+2-4,i*4

20、+2-4)= K(i*4+2-4,i*4+2-4)+Kzy(i);endformat long gAx WW=eig(inv(M)*K);f=sqrt(WW)/(2*pi);f0=diag(f);f00=abs(sort(f0)fn123=f00(1) f00(3) f00(5)wi1=54*2*pi;wi2=284*2*pi;yita1=0.1;yita2=0.2;alf=2*(yita2/wi2-yita1/wi1)*(1/wi12-1/wi12);beita=2*(yita2*wi2-yita1*wi1)/(1/wi12-wi12);C=alf*M+beita*K;D=C+G;Dzx(1

21、)=7e3;Dzy(1)=7e3;Dzx(12)=7e3;Dzy(12)=7e3;for i=1:1:N+1 D(i*4+1-4,i*4+1-4)=D(i*4+1-4,i*4+1-4)+Dzx(i); D(i*4+2-4,i*4+2-4)=D(i*4+2-4,i*4+2-4)+Dzy(i);endD=D*1;%该程序为newmark_newton算法xn=yn(:,n);dxn=dyn(:,n);ddxn=ddyn(:,n);if i=1xn=yn(:,n);endif i1xn=yn(:,n+1);endK0=K;K1=K*0;F(loca*4-4+1,1)=Famp1*cos(wi*t+f

22、ai(1); F(loca*4-4+2,1)=Famp1*sin(wi*t+fai(1); F(loca*4+1,1)=Famp1*cos(wi*t+fai(1);F(loca*4+2,1)=Famp1*sin(wi*t+fai(1);Pn1=F;Fr=Pn1*0;dert=20/1e6;k_rub=5e4;miu_rub=0.2;xx=yn(loc_rub*4-4+1,n)+5/1e6;yy=yn(loc_rub*4-4+2,n);rr=sqrt(xx2+yy2);if nFen*60&rrdert rub_sign=1; K1(loc_rub*4-4+1,loc_rub*4-4+1)=k_

23、rub*(rr-dert); K1(loc_rub*4-4+2,loc_rub*4-4+1)=k_rub*(rr-dert)*miu_rub; K1(loc_rub*4-4+1,loc_rub*4-4+2)=-k_rub*(rr-dert)*miu_rub; K1(loc_rub*4-4+2,loc_rub*4-4+2)=k_rub*(rr-dert); Fr(loc_rub*4-4+1,l)=k_rub*(rr-dert)/rr*(xx-miu_rub*yy); Fr(loc_rub*4-4+2,l)=k_rub*(rr-dert)/rr*(yy+miu_rub*xx);endKT=K0+K

24、1;If i=1 Fn1=K0*xn+Fr;endif i1 Fn1=K0*xn1+Fr;endKj=KT+a0*M+a1*C;Pj=Pn1+M*(a0*xn+a2*dxn+a3*ddxn)+C*(a1*xn+a4*dxn+a5*ddxn);If i=1 Pj=Pj-Fn1+KT*xn;endif i1 Pj=Pj-Fn1+KT*xn1;end if rr(20*Fen) & mod(i,Fen)=1 fprintf(f1,w/2/pi,y(1),y(2),y(3),y(4), y(5),y(6),y(7),y(8); end end%该程序主要计算 floquet稳定性并保存计算数据y2=eye(N); s0=y; for i=1:1:Fen*1 t=0; t=t+(i)*h; y=rkutta(t,h,y,w); At=fun_At(t,y,w); % 根据当前y值求A For j=1:1:N % 由dy2=A*y2积分求得y2 zzzz=rkutta21(t