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_k,Jp

2、_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 kkk<5000 if D2*D1>0 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*D2<0 wi=wi-step; step=step/2; wi=wi+step; D1,SS,Sn =surp

4、lus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr); EndD1;Wi;If atep<1/2000 Kkk=5000;end endCritical_speed=wi/2/pi*60 figure;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%for i=1:1:N+1 for

5、 j=1:1:2 for k=1:1:2 ud11(j,k.i)=0; ud12(j,k.i)=0; ud21(j,k.i)=0; ud22(j,k.i)=0; end endendfor i=1:1:Nfor j=1:1:2 for k=1:1:2 us11(j,k.i)=0; us12(j,k.i)=0; us21(j,k.i)=0; us22(j,k.i)=0; end endendfor i=1:1:Nfor j=1:1:2 for k=1:1:2 u11(j,k.i)=0; u12(j,k.i)=0; u21(j,k.i)=0; u22(j,k.i)=0; end endend%(b

6、)計算質點上傳遞矩陣點矩陣的一局部！%for i=1:1:N+1 ud11(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)計算質點上傳遞矩陣點矩陣的一局部！%for i=1:1:N+1ud12(1,1,i)=0;ud12(1,2,i)=(Jp_k(i)-Jd_k(i)*wi2;

7、 %考慮陀螺力矩ud12(2,1,i)=m_k(i)*wi2-k(i);ud12(2,2,i)=0;end%(d)以下計算的是無質量梁上的傳遞矩陣場矩陣%計算的錐軸的us是不對的，是隨便令的，在后面計算剩余量時，zhui中會把錯誤的覆蓋掉%for i=1:1:N us11(1,1,i)=1; us11(1,2,i)=1(i); us11 (2,1,i)=0; us11(2,2,i)=1; us12(1,1,i)=0; us12(1,2,i)=0; us12 (2,1,i)=0; us12(2,2,i)=0; us21(1,1,i)=1(i)2/(2*EI(i); us21(1,2,i)=(1(

8、i)3*(1-rr(i)/(6*EI(i); us21(2,1,i)=1(i)/EI(i); us21(2,2,i)=1(i)2/(2*EI(I); us22(1,1,i)=1; us22(1,2,i)=1(i); us22 (2,1,i)=0; us22(2,2,i)=1;end%此處全為計算中間量%for i=1:1:N2 Su (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)=

9、0; SS (2,1,i)=0; SS (2,2,i)=0;endfor i=1:1:2 for j=1:1:2 SS1(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_modify4，wi;cone_modify5，wi;cone_modify6，wi;cone_modify7，wi;cone_modify8，wi;cone_modify1

10、6，wi;cone_modify17，wi;cone_modify18，wi;cone_modify19，wi;cone_modify22，wi;cone_modify24，wi;%(f)形成最終傳遞矩陣%Ud11 Ud12 Ud21 Ud22 為最終參與計算的傳遞矩陣for i=1:1:N u11(:,:,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

11、)+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);for i=1:1:2 for j=1:1:2 SS1(i,j)=0; endendfor i=1:1:N+1ud11= u11(:,:,i); ud12= u12(:,:,i); ud21= u21(

12、:,:,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);for i=1:1:N+1 D1=D1*sign(det(Su(:,:,i); %消奇點end%(2)不平衡響應值EEEE:,:,n+2=-inv(SS(:,:,N+2)*PP(:,:,N+2);for i=N+1

13、:-1:1 EE(:,:,I)=Sn(:,:,i)*EE(:,:,i+1)-Sn(:,:,i)*UF(:,:,i);end A.2 碰摩轉子系統計算仿真程序%該程序主要完成完成jeffcott轉子圓周碰摩故障仿真%第一步：設置初始條件%rub_sign:碰摩標志，假設rub_sign=0,說明系統無碰摩故障；否那么rub_sign=1%loca: 不平衡質量的位置%loc_rub: 碰摩位置%Famp: 不平衡質量的大小單位為：g%wi: 轉速 單位為：rad%r: 偏心半徑 單位為：mm%Fampl: 離心力的大小 單位為：kg,m%fai: 不平衡量的初始相位 radclcclearrub

14、_sign loca loc_rub Famp wi r Famp1 fai=initial_conditions% 第二步：設置轉子系統的參數值%N： 劃分的軸段數%density: 軸的密度 單位為:kg/m3%Ef: 軸的彈性模量 單位為：Pa%L: 每個軸段的長度 單位為：m%R: 每個軸段的外半徑 單位為：m%ro: 每個軸段的內半徑 單位為：m%miu: 每個軸段的單元質量 kg/mN density Ef L R Ro miu=rotor_parameters% 第三步：設置移動單元質量距陣，移動單元質量距陣，剛度單元質量距陣和陀螺力距距陣%Mst: 移動單元質量距陣%Msr:

15、移動單元質量距陣%Ks: 剛度單元距陣%Ge: 陀螺力距單元距陣Mst: Msr Ks Ge=Mst_Msr_Ks_Ge(N,density,R,Ro.L.Ef,miu)% 第四步：距陣組集%M： 總的質量距陣%K： 總的剛度距陣%G: 總的陀螺力矩距陣M G K=M_G_K(N,Ef,R,Ro.Mst,Msr,Ge,Ks,miu,L)% 第五步：參加支撐剛度和阻尼K C D Ax=K_D(N,K,M,G)% 第六步：用Newmark方法進行計算%Fen: 每個周期內的步數%ht: 每步的長度ut1=;xt1=;yt1=;N3=(N+1)*4;Fen=100;ht=2*pi/wi/Fen;ga

16、ma=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;for i=1:1:N3 F(i,1)=0;endfor i=1:1;N3 yn(i:1)=0; dyn(i:1)=0; ddyn(i:1)=0;endt=0;for n=1:1:Fen*80 t=t+ht; n; for i=1:1:30 newmark_newton_

17、multiend if mod(n,100)=1 n end if n>Fen*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 loc_

18、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)/1

19、000*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*lel

20、l;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: 陀螺力

21、距單元距陣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; Ms

22、t(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

23、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; Ms

24、r(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)

25、=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)

26、;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(

27、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)2

28、;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;

29、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); en

30、dendend%該程序進行距陣組集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*

31、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)=7e

32、7;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+

33、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)

34、=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 i>1xn=yn(:,n+1);endK0=K;K1=K*0;F(loca*4-4+1,1)=Famp1*cos(wi*

35、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(xx2+yy2);if n>Fen*60&rr>dert rub_sign=1; K1(loc_rub*4-4+1,lo

36、c_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

37、*xx);endKT=K0+K1;If i=1 Fn1=K0*xn+Fr;endif i>1 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 i>1 Pj=Pj-Fn1+KT*xn1;end if rr<dert i=100;endQQ,RR=qr(Kj);xn1=RRQQ *Pj;ddxn1=a0*(xn1-xn)-a2*dxn-a3*ddxn;dxn1=dxn+a6*ddxn+a7*d

38、dxn1;yn(:,n+1)=xn1;dyn(:,n+1)=dxn1;ddyn(;,n+1)=ddxn1;A.3 穩定性分析程序% 第一步：設置初始條件N Fen y d y2=initial_conditionsfm=;for kk=1:1:60w=40+(kk-1)*1/1 %頻率區間h=1/w/Fen; %每個周期內積分點數為Fen=200w=2*pi*w;% 第二步：通過打靶法計算振幅的分岔特性shooting;% 第三步：計算 floquet穩定性并保存計算數據floquet_stabilityend%該程序主要通過打靶法計算振幅的分岔特征并存儲計算數據for i=1:1:30*Fe

39、n t=0; t=t+(i)*h; y=rkutta(t,h,y,w); if i>(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=rk

40、utta21(t,h,y2(:,j) ,At);endends=y;Rsi=s0-s;S=y2;Drds=eye(N)*1-S;dsi=Drds(-Rsi) ;y=s-dsi ;floquet_mul=max(abs(eig(y2);fm=fm;w/2/pi floquet_mul;fprintf(f2, %15.9f, %15.9fn ,w/2/pi, floquet_mul);function yout=rkutta(t,h,y,w)N=length(y);for i=1:1:N a(i)=0; d(i)=0; b(i)=0; y0(i)=0;enda(1)=h/2; a(2)=h/2;a

41、(3)=h; a(4)=h;d=fun(t,y,w);b=y;y0=y;for k=1:1:3 for i=1:1:N y(i)=y0(i)+a(k)*d(i); b(i)=b(i)+a(k+1)*d(i)/3; end tt=t+a(k); d=fun(tt,y,w);endfor i=1:1:N yout(i)=b(i)+h*d(i)/6;endfunction yout=ekutta21(t,h,y,A)N=length(y);For i=1:1:N a(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;

42、d=fun21(t,y,A);b=y;y0=y;for k=1:1:3 for i=1:1:N y(i)=y0(i)+a(k)*d(i); b(i)=b(i)+a(k+1)*d(i)/3; end tt=t+a(k); d=fun21(t,y,A);endfor i=1:1:N yout(i)=b(i)+h*d(i)/6;endfunction d=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.

43、695e-5; %圓盤2的偏心E=2.11e11; %轉軸的彈性模量delta1=200e-5; %碰摩間隙1delta2=2e-5; %碰摩間隙2kr1=1e5; %碰摩剛度1kr2=1e5; %碰摩剛度2fr1=0.1; %碰摩摩擦系數1fr2=0.1; %碰摩摩擦系數2I=pi*d0*d0*d0*d0/64;k11=2.3704e5;k22=2.3704e5;k12=-2.0741e5;k21=k12; % 剛度omega=w;g=0;% 轉子動力學方程If y(5)>=deltal Dirac1=1;endif y(5)<deltalDirac1=0;endif y(7)&

44、gt;delta2 Dirac2=1endif y(7)>delta2 Dirac2=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(

45、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(omeg

46、a*t)+fx2/m2;function At=fun_At(t,y,w)N=length(y);eps=1e-6;for i=1:1:N yi=y; yi(i)=y(i)+y(i)*eps; At0=fun(t,y,w); Ati=fun(t,yi,w); For j=1:1:N At(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

47、*alfa); T=(II*wi2/cos(alfa)*(2*cc*dd*sin(2*wi*t)/(1+dd*cos(2*wi*t); F (local*4-4+4,1)=T*sin(alfa)*cos(bita); % 不對中力 F (local*4-4+3,1)=T*sin(alfa)*sin(bita);A.5 轉子系統不對中故障的振動信號小波包分解程序%該程序主要研究應用基于小波包的能量比例計算方法提取其故障特征%讀取數據x=load(,r+);x1=x(1:13312,2);y1= x(1:13312,3);dt=0.00078;fs=1/dt;N=length(x1);T=dt*0:N-1; %時間軸% 時城波形plot(T,x1);xlabel(Time/sec,fontsize,8);ylabel(Disp. /um,fontsize,8);%軸心軌跡figur