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 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] =surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr); End D1; Wi; If atep<1/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)=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+1 ud12(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中会把错误的覆盖掉 %＝＝＝＝＝＝＝＝＝＝＝＝ 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(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:N＋2 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)=0; SS (2,1,i)=0; SS (2,2,i)=0; end for 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; end end %＝＝＝＝＝＝＝＝＝＝＝＝ %(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)形成最终传递矩阵 %＝＝＝＝＝＝＝＝＝＝＝＝ %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)+us22(:,:,i)*ud21(:,:,i); u22(:,:,i)=us21(:,:,i)*ud12(:,:,i)+us22(:,:,i)*ud22(:,:,i); end u11(:,:,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; end end for i=1:1:N+1 ud11= 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); for i=1:1:N+1 D1=D1*sign(det(Su(:,:,i)); %消奇点 end %＝＝＝＝＝＝(2)不平衡响应值EE＝＝＝＝＝＝ EE（:,:,n+2）=-inv(SS(:,:,N+2)*PP(:,:,N+2); for i=N+1:-1:1 EE(:,:,I)=Sn(:,:,i)*EE(:,:,i+1)-Sn(:,:,i)*UF(:,:,i); end

A.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] clc clear [rub_sign loca loc_rub Famp wi r Famp1 fai]=initial_conditions % 第二步：设置转子系统的参数值 %N： 划分的轴段数 %density: 轴的密度 单位为:[kg/m^3 %Ef: 轴的弹性模量 单位为：[Pa] %L: 每个轴段的长度 单位为：[m]