传递矩阵-matlab程序

%main_critical.m%该程序使用Riccati传递距阵法计算转子系统的临界转速及振型%本函数中均采用国际单位制%第一步：设置初始条件（调用函数shaft_parameters）%初始值设置包括：轴段数N，搜索次数M%输入轴段参数：内径d,外径D，轴段长度l，支撑刚度K，单元质量mm，极转动惯量Jpp[N,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_k,EI,rr]=shaft_pra_cal(N,D,d,l,Jpp,mm);%第三步：计算剩余量（调用函数surplus_calculate）,并绘制剩余量图%剩余量：D1fori=1:1:Mptx(i)=0;pty(i)=0;endforii=1:1:Mwi=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’);gridon%第四步：用二分法求固有频率及振型图%固有频率：Critical_speedwi=50;fori=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;whilekkk5000ifD2*D10wi=wi+step;D2=D1;[D1,SS,Sn]=surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr);endifD1*D20wi=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;Ifatep1/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;ud12(1,2,i)=(Jp_k(i)-Jd_k(i))*wi^2;%%%考虑陀螺力矩ud12(2,1,i)=m_k(i)*wi^2-k(i);ud12(2,2,i)=0;end%＝＝＝＝＝＝＝＝＝＝＝＝%(d)以下计算的是无质量梁上的传递矩阵―――场矩阵%计算的锥轴的us是不对的，是随便令的，在后面计算剩余量时，zhui中会把错误的覆盖掉%＝＝＝＝＝＝＝＝＝＝＝＝fori=1:1:Nus11(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(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%＝＝＝＝＝＝＝＝＝＝＝＝%此处全为计算中间量%＝＝＝＝＝＝＝＝＝＝＝＝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=0,说明系统无碰摩故障；否则rub_sign=1%loca:不平衡质量的位置%loc_rub:碰摩位置%Famp:不平衡质量的大小单位为：[g]%wi:转速单位为：[rad]%r:偏心半径单位为：[mm]%Fampl:离心力的大小单位为：[kg,m]%fai:不平衡量的初始相位[rad]clcclear[rub_signlocaloc_rubFampwirFamp1fai]=initial_conditions%第二步：设置转子系统的参数值%N：划分的轴段数%density:轴的密度单位为:[kg/m^3%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)==1nendifnFen*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(l