故障转子系统的非线性振动分析与诊断方法附录Amatlab程序.pdfVIP

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