!!故障转子系统的非线性振动分析与诊断方法附录A-matlab程序.docVIP

A.1 传递距阵法分析程序 %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）,并绘制剩余量图 %剩余量：D1 for i=1:1:M ptx(i)=0; pty(i)=0; end for 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)=w1 end ylabel(‘剩余量’)； plot(ptx,pty) xlabel(‘角速度red/s’); grid on % 第四步：用二分法求固有频率及振型图 %固有频率：Critical_speed wi=50; for i=1:1:4 order=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 kkk5000 if D2*D10 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*D20 wi=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); End D1; Wi; If atep1/2000 Kkk=5000; end end Critical_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 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 end end for i=1:1:N for 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 end end for i=1:1:N for 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 end end %＝＝＝＝＝＝＝＝＝＝＝＝ %(b)计算质点上传递矩阵―――点矩阵的一部分！ %＝＝＝＝＝＝＝＝＝＝＝＝ for i=1:1:N+1 ud11(1,1,i)=1; ud11(1,2,i)