【正文】 k,i)。for i=1:1:N+1 K(i*4+14,i*4+14)= K(i*4+14,i*4+14)+Kzx(i)。Dzy(12)=7e3。Fr=Pn1*0。If i==1 Fn1=K0*xn+Fr。for kk=1:1:60w=40+(kk1)*1/1 %頻率區(qū)間h=1/w/Fen。Rsi=s0s。a(3)=h。 a(2)=h/2。l2=。 % 剛度c11=。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/m2g。 dd=(1cos(2*alfa))/ (3+cos(2*alfa))。 ylabel(’ Displacement y/um’,’fontsize’,8)。s(3,:)=wprocef(t,[4 3])。s(15,:)=wprocef(t,[4 11])。 HHT變換用于分析碰摩轉(zhuǎn)子振動信號的程序%%讀數(shù)據(jù)x=load(‘’, ‘r+’)。%瞬時頻率譜Fai=[]。 endf=(0:NN1)/NN*(fs/2)。 IMFflag=1。c=[c c1’]。[i X(i)]]。 f(i,:)=instfreq(an(i.:)’,tt,1)’。Xa=[]。 breakendh1=h1k。X=X0。[im,tt]=toimage(A,fa)。 ylabel (I,’fontsize’,8)endsubplot(nIMF+1,1,Nimf+1) plot(t,rn)。ttt=[1:16]。s(14,:)=wprocef(t,[4 9])。s(2,:)=wprocef(t,[4 1])。xlabel(‘Time/sec’,’fontsize’,8)。endend 不對中轉(zhuǎn)子系統(tǒng)仿真程序%如下是考慮軸承不對中產(chǎn)生的附加載荷。fy2=fr2*kr2*(y(7)delta2)*Dirac2。k11=。 %圓盤1的質(zhì)量m1=。 d(i)=0。 b(i)=0。 y=rkutta(t,h,y,w)。dyn(:,n+1)=dxn1。 K1(loc_rub*44+2,loc_rub*44+2)=k_rub*(rrdert)。 F(loca*44+2,1)=Famp1*sin(wi*t+fai(1))。D=C+G。end%%該程序?qū)⒔M尼加入總體距陣function [K C D Ax]=K_D(N,K,M,G)Kzx(1)=7e7。 endendend%%該程序進(jìn)行距陣組集function [M G K]=M_G_K(N,Ef,R,R0,Mst,Msr,Ge,Ks,miu,L)% M: 總的質(zhì)量距陣% 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(7,5,i)=0。 Ks(5,4,i)=6*L(i)。Ge(8,6,i)=3*L(i)。 Ge(6,4,i)=3*L(i0)。Ge(2,1,i)=36。Msr(7,4,i)=0。 Msr(3,3,i)=4*L(i)^2。 Mst(8,1,i)=13*L(i)。 Mst(3,3,i)= 4*L(i)^2。Famp=[1]。 n。 beita=。 ud22= u22(:,:,i)。cone_modify（18，wi）。endfor i=1:1:2 for j=1:1:2 SS1(i,j)=0。 us22(1,1,i)=1。end%＝＝＝＝＝＝＝＝＝＝＝＝%(c)計(jì)算質(zhì)點(diǎn)上傳遞矩陣―――點(diǎn)矩陣的一部分！%＝＝＝＝＝＝＝＝＝＝＝＝for i=1:1:N+1ud12(1,1,i)=0。 end endendfor i=1:1:Nfor j=1:1:2 for k=1:1:2 u11(j,)=0。 wi=wi+step。 pty(i)=0。D2=D1。 ud12(j,)=0。 ud21(1,2,i)=0。 us12(1,2,i)=0。 Sn (1,2,i)=0。 endend%＝＝＝＝＝＝＝＝＝＝＝＝%(e)調(diào)用函數(shù)cone_modify修改錐軸的傳遞矩陣%＝＝＝＝＝＝＝＝＝＝＝＝cone_modify（4，wi）。 u12(:,:,N+1)=ud12(:,:,N+1)。end 碰摩轉(zhuǎn)子系統(tǒng)計(jì)算仿真程序%%該程序主要完成完成jeffcott轉(zhuǎn)子圓周碰摩故障仿真%＝＝＝＝＝＝＝＝＝＝＝第一步：設(shè)置初始條件%rub_sign:碰摩標(biāo)志，若rub_sign=0,說明系統(tǒng)無碰摩故障；否則rub_sign=1%loca: 不平衡質(zhì)量的位置%loc_rub: 碰摩位置%Famp: 不平衡質(zhì)量的大小單位為：[g]%wi: 轉(zhuǎn)速 單位為：[rad]%r: 偏心半徑 單位為：[mm]%Fampl: 離心力的大小 單位為：[kg,m]%fai: 不平衡量的初始相位 [rad]clcclear[rub_sign loca loc_rub Famp wi r Famp1 fai]=initial_conditions% 第二步：設(shè)置轉(zhuǎn)子系統(tǒng)的參數(shù)值%N： 劃分的軸段數(shù)%density: 軸的密度 單位為:[kg/m^3%Ef: 軸的彈性模量 單位為：[Pa]%L: 每個軸段的長度 單位為：[m]%R: 每個軸段的外半徑 單位為：[m]%ro: 每個軸段的內(nèi)半徑 單位為：[m]%miu: 每個軸段的單元質(zhì)量 [kg/m][N density Ef L R Ro miu]=rotor_parameters% 第三步：設(shè)置移動單元質(zhì)量距陣，移動單元質(zhì)量距陣，剛度單元質(zhì)量距陣和陀螺力距距陣%Mst: 移動單元質(zhì)量距陣%Msr: 移動單元質(zhì)量距陣%Ks: 剛度單元距陣%Ge: 陀螺力距單元距陣[Mst: Msr Ks Ge]=Mst_Msr_Ks_Ge(N,density,R,miu)% 第四步：距陣組集%M： 總的質(zhì)量距陣%K： 總的剛度距陣%G: 總的陀螺力矩距陣[M G K]=M_G_K(N,Ef,R,Msr,Ge,Ks,miu,L)% 第五步：加入支撐剛度和阻尼[K C D Ax]=K_D(N,K,M,G)% 第六步：用Newmark方法進(jìn)行計(jì)算%Fen: 每個周期內(nèi)的步數(shù)%ht: 每步的長度ut1=[]。for i=1:1:N3 F(i,1)=0。[t y’]]。for i=1:1。Mst(6,2,i)=54。 Mst(7,7,i)=4*L(i)^2。 Msr(6,1,i)=0。Msr(7,6,i)=3*L(i)。 Ge(5,1,i)=0。 Ge(5,5,i)=0。Ks(4,1,i)=6*L(i)。 Ks(,i)=6*L(i)。 Msr(j,k,i)=Msr(j,k,i)*miu(i)*R(i)^2/120/L(i)。 K(i*4+j4,i*4+k4)= K(i*4+j4,i*4+k4)+Ks(j+4,k+4,i1)。f00=abs(sort(f0))fn123=[f00(1) f00(3) f00(5)]wi1=54*2*pi。dxn=dyn(:,n)。yy=yn(loc_rub*44+2,n)。endif i1 Pj=PjFn1+KT*xn1。 y=rkutta(t,h,y,w)。floquet_mul=max(abs(eig(y2)))。for k=1:1:3 for i=1:1:N y(i)=y0(i)+a(k)*d(i)。y0=y。 %轉(zhuǎn)軸的彈性模量delta1=200e5。g=0。% function At=fun_At(t,y,w)N=length(y)。x1=x(1:13312,2)。x22=interp1(x11,x1,x2,’spline’)。p(8,:)=(abs(fft(s(8,:),freq))).^2。p(12,:)=(abs(fft(s(12,:),freq))).^2。fs=1/dt。 w1=[diff(Fai1)。 % 數(shù)據(jù) Ns=length(y)。breakendYa=spline(Xa(:,1)/fs,Xa(:,2),t)。 %殘余t=0:(1/fs) :((Ns1)/fs)。if min(size(x))==1 if size(x,2)==1 x=x’ 。 end Nmodes=1。rn=rn(Ns+1:2*Ns)。 m1=(Ya+Yb)/2。s2=s1。 w=[w,w1]。