传递矩阵-matlab程序

%main_critical.m%该程序使用Riccati传递距阵法计算转子系统的临界转速及振型%本函数中均采用国际单位制%第一步：设置初始条件(调用函数shaft_parameters)%初始值设置包括：轴段数N,搜索次数M%输入轴段参数：内径d,外径D,轴段长度1,支撑刚度K,单元质量mm,极转动惯量Jpp[N,M,d,D,l,K,mm,Jpp]=shaft_parameters;%第二步：计算单元的5个特征值(调用函数shaft_pra_ca1)%单元的5个特征值：%m_k：:质量%Jp_k:极转动惯量%Jd_k:直径转动惯量%EI:弹性模量与截面对中性轴的惯性矩的乘积%rr：剪切影响系数[m_k,Jp_k,EI,rr]=shaft_pra_ca1(N,D,d,1,Jpp,mm);%第三步：计算剩余量(调用函数surp1us_ca1cu1ate)，并绘制剩余量图%剩余量：D1fori=1:1:Mptx(i)=0;pty(i)=0;endforii=1:1:Mwi=ii/1*2+50;[D1,SS,Sn]=surp1us_ca1cu1ate(N,wi,K,m_k,Jp_k,JD_k,1,EI,rr);D1;pty(ii)=D1;ptx(ii)=w1endylabel(‘剩余量')；p1ot(ptx,pty)xlabel(‘角速度red/s');gridon%第四步：用二分法求固有频率及振型图%固有频率：Critica1_speedwi=50;fori=1:1:4order=i[D1,SS,Sn]=surp1us_ca1cu1ate(N,wi,k,m_k,Jp_k,Jd_k,1,EI,rr);Step=1;D2=D1;kkk=1;whi1ekkk<5000ifD2*D1>0wi=wi+step;D2=D1;[D1,SS,Sn]=surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr);endifD1*D2<0wi=wi-step;step=step/2;wi=wi+step;[D1,SS,Sn]=surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr);EndD1;Wi;Ifatep<1/2000Kkk=5000;endendCritical_speed=wi/2/pi*60figure;plot_mode(N,l,SS,Sn)wi=wi+2;end%surplus_calculate,.m%计算剩余量%(1)计算传递矩阵%(2)计算剩余量function[D1,SS,Sn]=surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr);%(1)计算传递矩阵%===============%(a)初值设为0%===============fori=1:1:N+1forj=1:1:2fork=1:1:2ud11(j,k.i)=0;ud12(j,k.i)=0;ud21(j,k.i)=0;ud22(j,k.i)=0;endendendfori=1:1:Nforj=1:1:2fork=1:1:2us11(j,k.i)=0;us12(j,k.i)=0;us21(j,k.i)=0;us22(j,k.i)=0;endendendfori=1:1:Nforj=1:1:2fork=1:1:2u11(j,k.i)=0;u12(j,k.i)=0;u21(j,k.i)=0;u22(j,k.i)=0;endendend%============%(b)计算质点上传递矩阵点矩阵的一部分！%============fori=1:1:N+1ud11(1,1,i)=1;ud11(1,2,i)=0;ud11(2,1,i)=0;ud11(2,2,i)=1;ud21(1,1,i)=0;ud21(1,2,i)=0;ud21(2,1,i)=0;ud21(2,2,i)=0;ud22(1,1,i)=1;ud22(1,2,i)=0;ud22(2,1,i)=0;ud22(2,2,i)=1;end%============%(c)计算质点上传递矩阵点矩阵的一部分！%============fori=1:1:N+1ud12(1,1,i)=0;udl2(l,2,i)=(Jp_k(i)-Jd_k(i))*wW2;%%%考虑陀螺力矩udl2(2,l,i)=m_k(i)*wiA2-k(i);udl2(2,2,i)=0;end%============%(d)以下计算的是无质量梁上的传递矩阵场矩阵%计算的锥轴的us是不对的，是随便令的，在后面计算剩余量时，zhui中会把错误的覆盖掉%============fori=l:l:Nusll(l,l,i)=l;usll(l,2,i)=l(i);usll(2,l,i)=0;usll(2,2,i)=l;usl2(l,l,i)=0;usl2(l,2,i)=0;usl2(2,l,i)=0;usl2(2,2,i)=0;us2l(l,l,i)=l(i)A2/(2*EI(i));us2l(l,2,i)=(l(i)A3*(l-rr(i))/(6*EI(i));us2l(2,l,i)=l(i)/EI(i);us2l(2,2,i)=l(i)A2/(2*EI(I));us22(l,l,i)=l;us22(l,2,i)=l(i);us22(2,l,i)=0;us22(2,2,i)=l;end%============%此处全为计算中间量%============fori=1:1:N＋2Su(1,1,i)=0;Su(1,2,i)=0;Su(2,1,i)=0;Su(2,2,i)=0;Sn(1,1,i)=0;Sn(1,2,i)=0;Sn(2,1,i)=0;Sn(2,2,i)=0;SS(1,1,i)=0;SS(1,2,i)=0;SS(2,1,i)=0;SS(2,2,i)=0;endfori=1:1:2forj=1:1:2SS1(i,j)=0;Ud11(i,j)=0;Ud12(i,j)=0;Ud21(i,j)=0;Ud22(i,j)=0;Us11(i,j)=0;Us12(i,j)=0;Us21(i,j)=0;Us22(i,j)=0;endend%============%(e)调用函数cone_modify修改锥轴的传递矩阵%============cone_modify(4，wi);cone_modify(5，wi);cone_modify(6，wi);cone_modify(7，wi);cone_modify(8，wi);cone_modify(16，wi);cone_modify(17，wi);cone_modify(18，wi);cone_modify(19，wi);cone_modify(22，wi);cone_modify(24，wi);%============%(f)形成最终传递矩阵%============%Ud11Ud12Ud21Ud22为最终参与计算的传递矩阵fori=1:1:Nu11(:,:,i)=us11(:,:,i)*ud11(:,:,i)+us12(:,:,i)*ud21(:,:,i);u12(:,:,i)=us11(:,:,i)*ud12(:,:,i)+us12(:,:,i)*ud22(:,:,i);u21(:,:,i)=us21(:,:,i)*ud11(:,:,i)+us22(:,:,i)*ud21(:,:,i);u22(:,:,i)=us21(:,:,i)*ud12(:,:,i)+us22(:,:,i)*ud22(:,:,i);endu11(:,:,N+1)=ud11(:,:,N+1);u12(:,:,N+1)=ud12(:,:,N+1);u21(:,:,N+1)=ud21(:,:,N+1);u22(:,:,N+1)=ud22(:,:,N+1);fori=1:1:2forj=1:1:2SS1(i,j)=0;endendfori=1:1:N+1ud11=u11(:,:,i);ud12=u12(:,:,i);ud21=u21(:,:,i);ud22=u22(:,:,i);SS(:,:,:i+1)=(ud11*SS1+ud12)*inv(ud21*SS1+ud22);Su(:,:,i)=ud21*SS1+ud22;Sn(:,:,i)=inv(ud21*SS1+ud22);%计算振型时用到SS1=SS(:,:,i+1);end%======(2)计算剩余量======D1=det(SS(:,:,N+2);fori=1:1:N+1D1=D1*sign(det(Su(:,:,i));%消奇点end%======(2)不平衡响应值EE======EE(:,:,n+2)=-inv(SS(:,:,N+2)*PP(:,:,N+2);fori=N+1:-1:1EE(:,:,I)=Sn(:,:,i)*EE(:,:,i+1)-Sn(:,:,i)*UF(:,:,i);endA.2碰摩转子系统计算仿真程序%main.m%该程序主要完成完成jeffcott转子圆周碰摩故障仿真%===========第一步：设置初始条件%rub_sign:碰摩标志，若rub_sign=O,说明系统无碰摩故障；否则rub_sign=l%loca:不平衡质量的位置%loc_rub:碰摩位置%Famp:不平衡质量的大小单位为：[g]%wi:转速单位为：[rad]%r:偏心半径单位为：[mm]%Fampl:离心力的大小单位为：[kg,m]%fai:不平衡量的初始相位[rad]clcclear[rub_signlocaloc_rubFampwirFamplfai]=initial_conditions%第二步：设置转子系统的参数值%N：划分的轴段数%density:轴的密度单位为:[kg/mT%Ef:轴的弹性模量单位为：[Pa]%L:每个轴段的长度单位为：[m]%R:每个轴段的外半径单位为：[m]%ro:每个轴段的内半径单位为：[m]%miu:每个轴段的单元质量[kg/m][NdensityEfLRRomiu]=rotor_parameters%第三步：设置移动单元质量距阵，移动单元质量距阵，刚度单元质量距阵和陀螺力距距阵%Mst:移动单元质量距阵%Msr:移动单元质量距阵%Ks:刚度单元距阵%Ge:陀螺力距单元距阵[Mst:MsrKsGe]=Mst_Msr_Ks_Ge(N,density,R,Ro.L.Ef,miu)%第四步：距阵组集%M：总的质量距阵%K：总的刚度距阵%G:总的陀螺力矩距阵[MGK]=M_G_K(N,Ef,R,Ro.Mst,Msr,Ge,Ks,miu,L)%第五步：加入支撑刚度和阻尼[KCDAx]=K_D(N,K,M,G)%第六步：用Newmark方法进行计算%Fen:每个周期内的步数%ht:每步的长度ut1=[];xt1=[];yt1=[];N3=(N+1)*4;Fen=100;ht=2*pi/wi/Fen;gama=0.5;beita=0.25;a0=1.0/(beita*ht);al=gama/(beita*ht);a2=1.0/(beita*ht);a3=0.5/beita-1.0;a4=gama/beita-1.0;a5=ht/2.0*(gama/beita-2.0);a6=ht*(1.0-gama);a7=gama*ht;fori=1:1:N3F(i,1)=0;endfori=1:1;N3yn(i:1)=0;dyn(i:1)=0;ddyn(i:1)=0;endt=0;forn=1:1:Fen*80t=t+ht;n;fori=1:1:30newmark_newton_multiendifmod(n,100)==1nendifn>Fen*60fori=1:1:N+1x(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*180endu=F(loca*4-4+1,1);ut1=[ut1;[tu]];xt1=[xt1;[tx']];yt1=[yt1;[ty']];endendrub_signsave'xtl.dat'xtl-ascii;save'ytl.dat'ytl-ascii;save'utl.dat'utl-ascii;%initial_conditions.m%该程序主要进行仿真条件设置Function[rub_signlocaloc_robFampwirFamplfai]=initial_conditions%需要设置的初始条件有：%rub_sign:碰摩标志，若rub_sign=O,说明系统无碰摩故障；否则rub_sign=1%loca:不平衡质量的位置%loc_rub:碰摩位置%Famp:不平衡质量的大小单位为：[g]%wi:转速单位为：[rad]%r：偏心半径单位为:[mm]%Famp1：离心力的大小单位为:[kg,m]%fai:不平衡量的初始相位[rad]rub_sign=0;loca=6;loc_rub=8;Famp=[l];wi=3000/60*2*pi;r=30Famp1=Famp(1)/1000*wiA2*r/1000;fai=[3030]/l80*pi%roto_parameters.m%该程序对Jeffcott转子系统进行参数设置function[NdensityEfLRR0miu]=rotor_parameters%N:划分的轴段数%density:轴的密度单位为：[kg/mA3]

%Ef:%L%R%R0:%miu:轴的弹性模量单位为:[Pa]每个轴段的长度单位为:[m]每个轴段的外半径单位为:[m]每个轴段的内半径单位为:[m]每个轴段的单元质量[kg/m]N=11;Density=7850;Ef=[2.12.12.12.12.12.12.12.12.12.12.1]*lell;L=[90.590.580.562.5302022.562.590.590.590.5]/1000;R=[2020202020902020202020]/2000;R0=[00000000000]/2000;fori=1:1;Nmiu(i)=density*pi*(R(i)人2-RO(i)人2)end%Mst_Msr_Ks_Ge.m%该程序设置单元距阵Function[MstMsrKsGe]=Mst_Msr_Ks_Ge(N,density,R,R0,L,Ef,miu)%Mst:移动单元质量距阵%Msr:转动单元质量距阵%Ks:刚度单元距阵%Ge:陀螺力距单元距阵NN0=N;NN1=NN0+1NN2=NN1+1fori=1:1:NN0Mst(1,1,i)=156;Mst(2,1,i)=0;fori=1:1:NN0Mst(1,1,i)=156;Mst(2,1,i)=0;Mst(3,1,i)=0;Mst(4,1,i)22*L(i);Mst(4,4,i)=4*L(i)人2;Mst(5,3,i)=0;Mst(6,2,i)=54;Mst(7,1,i)=0;Mst(7.4.i)=0;Mst(8,3,i)=0;Mst(6,5.i)=0;Mst(7,5,i)=0;Mst(8,5,i)=-22*L(i);Mst(8,8,i)=4*L(i)A2;endfori=1:1:NN0Msr(1,1,i)=36;Msr(2,1,i)=0;Msr(3,1,i)=0Msr(3,3,i)=4*L(i)A2;Mst(2,2,i)=156;Mst(3,2,i)=-22*L(i);Mst(4,2,i)=0;Mst(5,1,i)=54;Mst(5,4,i)=13*L(i);Mst(6,3,i)=13*L(i);Mst(7,2,i)=13*L(i);Mst(8,1,i)=-13*L(i);Mst(8,4,i)=-3*L(i)A2;Mst(6,6,i)=156;Mst(7,6,i)=22*L(i);Mst(8,6,i)=0Msr(2,2,i)=36;Msr(3,2,i)=-3*L(i);Mst(3,3,i)=4*L(i)A2;Mst(4,3,i)=0;Mst(5,2,i)=0Mst(6,1,i)=0;Mst(6,4,i)=0;Mst(7,3,i)=-3*L(i)A2;Mst(8,2,i)=0;Mst(5,5,i)=156;Mst(7,7,i)=4*L(i)A2;Mst(8,7,i)=0Msr(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;Msr(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)A2;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)A2;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)A2;Msr(8,5,i)=-3*L(i);Msr(8,6,i)=0;Msr(8,7,i)=0;Msr(8,8,i)=4*L(i)A2;endfori=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)A2;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);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)A2;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)A2;Ge(8,8,i)=0;endfori=1:1:NN0Ks(1,1,i)=12;Ks(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)A2;Ks(4,1,i)=6*L(i);Ks(4,2,i)=0;Ks(4,3,i)=0;Ks(4,4,i)=4*L(i)A2;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)A2;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)A2;Ks(8,5,i)=-6*L(i);Ks(8,6,i)=0;Ks(8,7,i)=0;Ks(8,8,i)=4*L(i)A2;endfori=1:1:NN0forj=1:1:8fork=1:1:8EI=Ef(i)*pi*(R(i)人4-R0(i)人4-)/4;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)A3;endendendfori=1:1:NN0forj=1:1:8fork=j:1:8Mst(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);endendend%M_G_K.m%该程序进行距阵组集function[MGK]=M_G_K(N,Ef,R,R0,Mst,Msr,Ge,Ks,miu,L)%M:总的质量距阵%K：总的刚度距阵%G：总的陀螺力距距阵NN0=Nfori=1:1:NN0forj=1:1:8forK=1:1:8Ms(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);endendengfori=1:1:Nforj=1:1:8forK=1:1:8M(i*4+j-4,i*4+k-4)=Ms(j,k,i);G(i*4+j-4,i*4+k-4)=Ge(j,k,i);K(i*4+j-4,i*4+k-4)=K(j,k,i);endendendfori=2:1:Nforj=1:1;4fork=1:1:4M(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);endendendKzx(1)=7e7;Kzy(1)=7e7;Kzx(12)=7e7;Kzy(12)=7e7;fori=1:1:N+1K(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%K_D.m%该程序将组尼加入总体距阵function[KCDAx]=K_D(N,K,M,G)Kzx(1)=7e7;Kzy(1)=7e7;Kzx(12)=7e7;Kzy(12)=7e7;fori=1:1:N+1K(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);endformatlongg[AxWW]=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-yital/wil)*(l/wil^2-l/wilA2);beita=2*(yita2*wi2-yital*wil)/(l/wilA2-wilA2);C=alf*M+beita*K;D=C+G;Dzx(l)=7e3;Dzy(l)=7e3;Dzx(l2)=7e3;Dzy(12)=7e3;fori=1:1:N+1D(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_multi.m%该程序为newmark_newton算法xn=yn(:,n);dxn=dyn(:,n);ddxn=ddyn(:,n);ifi==1xn=yn(:,n);endifi>1xn=yn(:,n+1);endK0=K;K1=K*0;F(loca*4-4+1,1)=Famp1*cos(wi*t+fai(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(xxA2+yyA2);ifn>Fen*60&rr>dertrub_sign=1;K1(loc_rub*4-4+1,loc_rub*4-4+1)=k_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+K1;Ifi==1Fn1=K0*xn+Fr;endifi>1Fn1=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);Ifi==1Pj=Pj-Fn1+KT*xn;endifi>1Pj=Pj-Fn1+KT*xn1;endifrr<derti=100;end[QQ,RR]=qr(Kj);'*Pj;ddxn1=a0*(xn1-xn)-a2*dxn-a3*ddxn;dxn1=dxn+a6*ddxn+a7*ddxn1;yn(:,n+1)=xn1;dyn(:,n+1)=dxn1;ddyn(;,n+1)=ddxn1;A.3稳定性分析程序%第一步：设置初始条件[NFenydy2]=initial_conditionsfm=[];forkk=1:1:60w=40+(kk-1)*1/1%频率区间h=1/w/Fen;%每个周期内积分点数为Fen=200w=2*pi*w;%第二步：通过打靶法计算振幅的分岔特性shooting;%第三步：计算floquet稳定性并保存计算数据floquet_stabilityend%shooting.m%该程序主要通过打靶法计算振幅的分岔特征并存储计算数据fori=1:1:30*Fent=0;t=t+(i)*h;y=rkutta(t,h,y,w);ifi>(20*Fen)&mod(i,Fen)==1fprintf(f1,w/2/pi,y(1),y(2),y(3),y(4),y(5),y(6),y(7),y(8));endend%floquet_stability.m%该程序主要计算floquet稳定性并保存计算数据y2=eye(N);s0=y;fori=1:1:Fen*1t=0;t=t+(i)*h;y=rkutta(t,h,y,w);At=fun_At(t,y,w);%根据当前y值求AForj=1:1:N%由dy2=A*y2积分求得y2zzzz=rkutta21(t,h,y2(:,j)',At);endends=y;Rsi=s0-s;S=y2;Drds=eye(N)*1-S;';y=s-dsi';floquet_mul=max(abs(eig(y2)));fm=[fm;w/2/pifloquet_mul];fprintf(f2,'%15.9f,%15.9f',w/2/pi,floquet_mul);%rkutta.mfunctionyout=rkutta(t,h,y,w)N=length(y);fori=1:1:Na(i)=0;d(i)=0;b(i)=0;y0(i)=0;enda(1)=h/2;a(2)=h/2;a(3)=h;a(4)=h;d=fun(t,y,w);b=y;y0=y;fork=1:1:3fori=1:1:Ny(i)=y0(i)+a(k)*d(i);b(i)=b(i)+a(k+1)*d(i)/3;endtt=t+a(k);d=fun(tt,y,w);endfori=1:1:Nyout(i)=b(i)+h*d(i)/6;end%rkutta21.mfunctionyout=ekutta21(t,h,y,A)N=length(y);Fori=1:1:Na(i)=0;d(i)=0;b(i)=0;y0(i)=0;enda(1)=h/2;a(2)=h/2;a(3)=h;a(4)=h;d=fun21(t,y,A);b=y;y0=y;fork=1:1:3fori=1:1:Ny(i)=y0(i)+a(k)*d(i);b(i)=b(i)+a(k+1)*d(i)/3;endtt=t+a(k);d=fun21(t,y,A);endfori=1:1:Nyout(i)=b(i)+h*d(i)/6;end%fun.mfunctiond=fun(t,y,w)N=length(y);%圆盘1的质量m1=0.7902;%圆盘2的质量m2=0.7902;%转轴的直径d0=0.01;l1=0.16;l2=0.16;l3=0.16;e1=5.695e-3*0;%圆盘1的偏心e2=5.695e-5;%圆盘2的偏心E=2.11e11;%转轴的弹性模量delta1=200e-5;delta2=2e-5;kr1=1e5;kr2=1e5;delta1=200e-5;delta2=2e-5;kr1=1e5;kr2=1e5;fr1=0.1;fr2=0.1;%碰摩间隙2%碰摩刚度1%碰摩刚度2%碰摩摩擦系数1%碰摩摩擦系数2I=pi*d0*d0*d0*d0/64;k11=2.3704e5;k22=2.3704e5;k12=-2.0741e5;k21=k12;%刚度c11=53.118385200000006;c12=-37.333800000000004;c21=-37.333800000000004;c22=53.118385200000006;%//阻尼系数omega=w;g=0;%转子动力学方程Ify(5)>=deltalDirac1=1;endify(5)<deltalDirac1=0;endify(7)>delta2Dirac2=1endify(7)>delta2Dirac2=0;endfx1=-kr1*(y(5)-delta1)*Dirac1;fx2=-kr2*(y(7)-delta2)*Dirac2;fy1=-fr1*kr1*(y(5)-delta1)*Dirac1;fy2=-fr2*kr2*(y(7)-delta2)*Dirac2;d(1)=y(2);d(2)=-(c11/m1)*y(2)-(c12/m1)*y(4)-(k11/m1)*y(1)-(k12/m1)*y(3)+e1*omega*omega*sin(omega*t)+fy1/m1-g;d(3)=y(4);d(4)=-(c21/m2)*y(2)-(c22/m2)*y(4)-(k21/m2)*y(1)-(k22/m2)*y(3)+e2*omega*omega*sin(omega*t)+fy2/m2-g;d(5)=y(6);d(6)=-(c11/m1)*y(6)-(c12/m1)*y(8)-(k11/m1)*y(5)-(k12/m1)*y(7)+e1*omega*omega*cos(omega*t)+fy1/m1;d(7)=y(8);d(8)=-(c21/m2)*y(6)-(c22/m2)*y(8)-(k21/m2)*y(5)-(k22/m2)*y(7+e2*omega*omega*cos(omega*t)+fx2/m2;%fun_At.mfunctionAt=fun_At(t,y,w)N=length(y);eps=1e-6;fori=1:1:Nyi=y;yi(i)=y(i)+y(i)*eps;At0=fun(t,y,w);Ati=fun(t,yi,w);Forj=1:1:NAt(j,i)=(Ati(j)-At0(j))/(yi(i)*eps);endendA.4不对中转子系统仿真程序%如下是考虑轴承不对中产生的附加载荷。alfa=pi/12;bita=pi/15;cc=4*cos(alfa)/(3+cos(2*alfa));dd=(1-cos(2*alfa))/(3+cos(2*alfa));T=(II*wiA2/cos(alfa))*(2*cc*dd*sin(2*wi*t)/(l+dd*cos(2*wi*t)));F(local*4-4+4,1)=T*sin(alfa)*cos(bita);%不对中力F(local*4-4+3,l)=T*sin(alfa)*sin(bita);A.5转子系统不对中故障的振动信号小波包分解程序%WPDmisalignment.m%该程序主要研究应用基于小波包的能量比例计算方法提取其故障特征%读取数据x=load(‘misalignment25hz.txt','r+');xl=x(l:l33l2,2);yl=x(l:l33l2,3);dt=0.00078;fs=l/dt;N=length(xl);T=dt*[0:N-l];%时间轴%时城波形plot(T,xl);xlabel(‘Time/sec','fontsize',8);ylabel(‘Disp./um','fontsize',8);%轴心轨迹figure(2)plot(x1,y1,);xlabel('Displacementx/um','fontsize',8);ylabel('Displacementy/um','fontsize',8);%对原信号进行重采样%re-samplingfrom250Hzinto128Hzfs1=1024;dt1=1/fs1;x11=1:N;x2=1:fs/fs1:N;x22=interp1(x11,x1,x2,'spline');Ns=length(x22);T=dt1*[0:Ns-1];figure(2)plot(T,x22)%小波包分解freq=fs1/(Ns)*(0:Ns/2-1);t=wpdec(x22,4,'db44');%基于小波包分解的能量比例计算s(1,:)=wprocef(t,[40]);p(l,:)=(abs(fft(s(l,:),freq))).人2;s(2,:)=wprocef(t,[41]);p(2,:)=(abs(fft(s(2,:),freq))).A2;s(4,:)=wprocef(t,[42]);p(4,:)=(abs(fft(s(4,:),freq))).A2;s(3,:)=wprocef(t,[43]);p(3,:)=(abs(fft(s(3,:),freq))).A2;s(7,:)=wprocef(t,[44]);p(7,:)=(abs(fft(s(7,:),freq))).A2;s(8,:)=wprocef(t,[45]);p(8,:)=(abs(fft(s(8,:),freq))).A2;s(6,:)=wprocef(t,[46]);p(6,:)=(abs(fft(s(6,:),freq))).A2;s(5,:)=wprocef(t,[47]);p(5,:)=(abs(fft(s(5,:),freq))).A2;s(l3,:)=wprocef(t,[48]);p(l3,:)=(abs(fft(s(l3,:),freq))).A2;s(l4,:)=wprocef(t,[49]);p(l4,:)=(abs(fft(s(l4,:),freq))).A2;s(l6,:)=wprocef(t,[4l0]);p(l6,:)=(abs(fft(s(l6,:),freq))).A2;s(l5,:)=wprocef(t,[4ll]);p(l5,:)=(abs(fft(s(l5,:),freq))).A2;s(ll,:)=wprocef(t,[4l2]);p(ll,:)=(abs(fft(s(ll,:),freq))).人2;s(12,:)=wprocef(t,[413]);p(12,:)=(abs(fft(s(12,:),freq))).A2;s(l0,:)=wprocef(t,[4l4]);p(l0,:)=(abs(fft(s(l0,:),freq))).A2;s(9,:)=wprocef(t,[4l5]);p(9,:)=(abs(fft(s(9,:),freq))).A2;ps=sum(p);pp=p/psttt=[l:l6]*32;ttt=[l:l6];bar(ttt,pp,'r');ylabel(‘能量比','fontsize',10,‘fontweight',‘bold');xlabel(‘频带'，‘fontsize',10,‘fontweight',‘bold');A・6HHT变换用于分析碰摩转子振动信号的程序%HHTrubbing.m%读数据x=load(‘healthy.txt',‘r+');X=(x/78.7)*1000;x=X(:,2);y=X(:,3);dt=0.00039;fs=1/dt;N=length(x);t=dt*[0:N-1];%时间轴%EMD分解[c,rn,nIMF]=emd(x');figure(4)title(‘IMFsandresidual')fori=1;nIMFsubplot(nIMF+1,1,i)plot(t,c(:,i));axistight;ylabel(I,'fontsize',8)endsubplot(nIMF+1,1,Nimf+1)plot(t,rn);axistight;ylabel(‘rn',‘fontsiz