傳遞矩陣-matlab程序

%main_critical.m%該程序使用Riccati傳遞距陣法計算轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速及振型%本函數(shù)中均采用國際單位制%第一步：設(shè)置初始條件(調(diào)用函數(shù)shaft_parameters)%初始值設(shè)置包括：軸段數(shù)N,搜索次數(shù)M%輸入軸段參數(shù)：內(nèi)徑d,外徑D,軸段長度1,支撐剛度K,單元質(zhì)量mm,極轉(zhuǎn)動慣量Jpp[N,M,d,D,l,K,mm,Jpp]=shaft_parameters;%第二步：計算單元的5個特征值(調(diào)用函數(shù)shaft_pra_ca1)%單元的5個特征值：%m_k：:質(zhì)量%Jp_k:極轉(zhuǎn)動慣量%Jd_k:直徑轉(zhuǎn)動慣量%EI:彈性模量與截面對中性軸的慣性矩的乘積%rr：剪切影響系數(shù)[m_k,Jp_k,EI,rr]=shaft_pra_ca1(N,D,d,1,Jpp,mm);%第三步：計算剩余量(調(diào)用函數(shù)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)初值設(shè)為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)計算質(zhì)點(diǎn)上傳遞矩陣點(diǎn)矩陣的一部分！%============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)計算質(zhì)點(diǎn)上傳遞矩陣點(diǎn)矩陣的一部分！%============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)以下計算的是無質(zhì)量梁上的傳遞矩陣場矩陣%計算的錐軸的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)調(diào)用函數(shù)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));%消奇點(diǎn)end%======(2)不平衡響應(yīng)值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碰摩轉(zhuǎn)子系統(tǒng)計算仿真程序%main.m%該程序主要完成完成jeffcott轉(zhuǎn)子圓周碰摩故障仿真%===========第一步：設(shè)置初始條件%rub_sign:碰摩標(biāo)志，若rub_sign=O,說明系統(tǒng)無碰摩故障；否則rub_sign=l%loca:不平衡質(zhì)量的位置%loc_rub:碰摩位置%Famp:不平衡質(zhì)量的大小單位為：[g]%wi:轉(zhuǎn)速單位為：[rad]%r:偏心半徑單位為：[mm]%Fampl:離心力的大小單位為：[kg,m]%fai:不平衡量的初始相位[rad]clcclear[rub_signlocaloc_rubFampwirFamplfai]=initial_conditions%第二步：設(shè)置轉(zhuǎn)子系統(tǒng)的參數(shù)值%N：劃分的軸段數(shù)%density:軸的密度單位為:[kg/mT%Ef:軸的彈性模量單位為：[Pa]%L:每個軸段的長度單位為：[m]%R:每個軸段的外半徑單位為：[m]%ro:每個軸段的內(nèi)半徑單位為：[m]%miu:每個軸段的單元質(zhì)量[kg/m][NdensityEfLRRomiu]=rotor_parameters%第三步：設(shè)置移動單元質(zhì)量距陣，移動單元質(zhì)量距陣，剛度單元質(zhì)量距陣和陀螺力距距陣%Mst:移動單元質(zhì)量距陣%Msr:移動單元質(zhì)量距陣%Ks:剛度單元距陣%Ge:陀螺力距單元距陣[Mst:MsrKsGe]=Mst_Msr_Ks_Ge(N,density,R,Ro.L.Ef,miu)%第四步：距陣組集%M：總的質(zhì)量距陣%K：總的剛度距陣%G:總的陀螺力矩距陣[MGK]=M_G_K(N,Ef,R,Ro.Mst,Msr,Ge,Ks,miu,L)%第五步：加入支撐剛度和阻尼[KCDAx]=K_D(N,K,M,G)%第六步：用Newmark方法進(jìn)行計算%Fen:每個周期內(nèi)的步數(shù)%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%該程序主要進(jìn)行仿真條件設(shè)置Function[rub_signlocaloc_robFampwirFamplfai]=initial_conditions%需要設(shè)置的初始條件有：%rub_sign:碰摩標(biāo)志，若rub_sign=O,說明系統(tǒng)無碰摩故障；否則rub_sign=1%loca:不平衡質(zhì)量的位置%loc_rub:碰摩位置%Famp:不平衡質(zhì)量的大小單位為：[g]%wi:轉(zhuǎn)速單位為：[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%該程序?qū)effcott轉(zhuǎn)子系統(tǒng)進(jìn)行參數(shù)設(shè)置function[NdensityEfLRR0miu]=rotor_parameters%N:劃分的軸段數(shù)%density:軸的密度單位為：[kg/mA3]

%Ef:%L%R%R0:%miu:軸的彈性模量單位為:[Pa]每個軸段的長度單位為:[m]每個軸段的外半徑單位為:[m]每個軸段的內(nèi)半徑單位為:[m]每個軸段的單元質(zhì)量[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%該程序設(shè)置單元距陣Function[MstMsrKsGe]=Mst_Msr_Ks_Ge(N,density,R,R0,L,Ef,miu)%Mst:移動單元質(zhì)量距陣%Msr:轉(zhuǎn)動單元質(zhì)量距陣%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%該程序進(jìn)行距陣組集function[MGK]=M_G_K(N,Ef,R,R0,Mst,Msr,Ge,Ks,miu,L)%M:總的質(zhì)量距陣%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%該程序?qū)⒔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穩(wěn)定性分析程序%第一步：設(shè)置初始條件[NFenydy2]=initial_conditionsfm=[];forkk=1:1:60w=40+(kk-1)*1/1%頻率區(qū)間h=1/w/Fen;%每個周期內(nèi)積分點(diǎn)數(shù)為Fen=200w=2*pi*w;%第二步：通過打靶法計算振幅的分岔特性shooting;%第三步：計算floquet穩(wěn)定性并保存計算數(shù)據(jù)floquet_stabilityend%shooting.m%該程序主要通過打靶法計算振幅的分岔特征并存儲計算數(shù)據(jù)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穩(wěn)定性并保存計算數(shù)據(jù)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);%根據(jù)當(dāng)前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的質(zhì)量m1=0.7902;%圓盤2的質(zhì)量m2=0.7902;%轉(zhuǎn)軸的直徑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;%轉(zhuǎn)軸的彈性模量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%碰摩摩擦系數(shù)1%碰摩摩擦系數(shù)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;%//阻尼系數(shù)omega=w;g=0;%轉(zhuǎn)子動力學(xué)方程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不對中轉(zhuǎn)子系統(tǒng)仿真程序%如下是考慮軸承不對中產(chǎn)生的附加載荷。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轉(zhuǎn)子系統(tǒng)不對中故障的振動信號小波包分解程序%WPDmisalignment.m%該程序主要研究應(yīng)用基于小波包的能量比例計算方法提取其故障特征%讀取數(shù)據(jù)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);%對原信號進(jìn)行重采樣%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變換用于分析碰摩轉(zhuǎn)子振動信號的程序%HHTrubbing.m%讀數(shù)據(jù)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