传递矩阵-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）,并绘制剩余量图

%剩余量：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]

%R: 每个轴段的外半径单位为：[m]

%ro: 每个轴段的内半径单位为：[m]

%miu: 每个轴段的单元质量[kg/m]

[N density Ef L R Ro miu]=rotor_parameters

% 第三步：设置移动单元质量距阵，移动单元质量距阵，刚度单元质量

距阵和陀螺力距距阵

%Mst: 移动单元质量距阵

%Msr: 移动单元质量距阵

%Ks: 刚度单元距阵

%Ge: 陀螺力距单元距阵

[Mst: Msr Ks Ge]=Mst_Msr_Ks_Ge(N,density,R,Ro.L.Ef,miu)

% 第四步：距阵组集

%M：总的质量距阵

%K：总的刚度距阵

%G: 总的陀螺力矩距阵

[M G K]=M_G_K(N,Ef,R,Ro.Mst,Msr,Ge,Ks,miu,L)

% 第五步：加入支撑刚度和阻尼

[K C D Ax]=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;

for i=1:1:N3

F(i,1)=0;

end

for i=1:1;N3

yn(i:1)=0;

dyn(i:1)=0;

ddyn(i:1)=0;

end

t=0;

for n=1:1:Fen*80

t=t+ht;

n;

for i=1:1:30

newmark_newton_multi

end

if mod(n,100)==1

n

end

if n>Fen*60

for i=1:1:N+1

x(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*180

end

u=F(loca*4-4+1,1);

ut1=[ut1;[t u]];

xt1=[xt1;[t x’]];

yt1=[yt1;[t y’]];

end

end

rub_sign

save ’xt1.dat’ xt1 –ascii;

save ’yt1.dat’ yt1 –ascii;

save ’ut1.dat’ ut1 –ascii;

%initial_conditions.m

%该程序主要进行仿真条件设置

Function [rub_sign loca loc_rob Famp wi r Famp1 fai]=initial_conditions

%需要设置的初始条件有：

%rub_sign: 碰摩标志，若rub_sign=0,说明系统无碰摩故障；否则rub_sign=1 %loca: 不平衡质量的位置

%loc_rub: 碰摩位置

%Famp: 不平衡质量的大小单位为：[g]

%wi: 转速单位为：[rad]

% r：偏心半径单位为:[mm]

%Famp1：离心力的大小单位为:[kg,m]

%fai: 不平衡量的初始相位[rad]

rub_sign=0;

loca=6;

loc_rub=8;

Famp=[1];

wi=3000/60*2*pi;

r=30

Famp1=Famp(1)/1000*wi^2*r/1000;

fai=[30 30]/180*pi

%roto_parameters.m

%该程序对Jeffcott转子系统进行参数设置

function [N density Ef L R R0 miu]=rotor_parameters

%N: 划分的轴段数

%density: 轴的密度单位为：[kg/m^3]

%Ef: 轴的弹性模量单位为:[Pa]

%L 每个轴段的长度单位为:[m]

%R 每个轴段的外半径单位为:[m]

%R0: 每个轴段的内半径单位为:[m]

%miu: 每个轴段的单元质量[kg/m]

N=11;

Density=7850;

Ef=[2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1]*lell;

L=[90.5 90.5 80.5 62.5 30 20 22.5 62.5 90.5 90.5 90.5]/1000;

R=[20 20 20 20 20 90 20 20 20 20 20]/2000;

R0=[0 0 0 0 0 0 0 0 0 0 0]/2000;

for i=1:1;N

miu(i)=density*pi*(R(i)^2-R0(i)^2)

end

%Mst_Msr_Ks_Ge.m

%该程序设置单元距阵

Function [Mst Msr Ks Ge]= Mst_Msr_Ks_Ge(N,density,R,R0,L,Ef,miu)

%Mst: 移动单元质量距阵

%Msr: 转动单元质量距阵

%Ks: 刚度单元距阵

%Ge: 陀螺力距单元距阵

NN0=N;

NN1=NN0+1

NN2=NN1+1

for i=1:1:NN0

Mst(1,1,i)=156;

Mst(2,1,i)=0; Mst(2,2,i)=156;

Mst(3,1,i)=0; Mst(3,2,i)=-22*L(i); Mst(3,3,i)= 4*L(i)^2;

Mst(4,1,i)22*L(i); Mst(4,2,i)=0; Mst(4,3,i)=0;

Mst(4,4,i)=4*L(i)^2; Mst(5,1,i)=54; Mst(5,2,i)=0

Mst(5,3,i)=0; Mst(5,4,i)=13*L(i); Mst(6,1,i)=0;

Mst(6,2,i)=54; Mst(6,3,i)=13*L(i); Mst(6,4,i)=0;

Mst(7,1,i)=0; Mst(7,2,i)=13*L(i); Mst(7,3,i)=-3*L(i)^2;

Mst(7.4.i)=0; Mst(8,1,i)=-13*L(i); Mst(8,2,i)=0;

Mst(8,3,i)=0; Mst(8,4,i)=-3*L(i)^2; Mst(5,5,i)=156;

Mst(6,5.i)=0; Mst(6,6,i)=156;

Mst(7,5,i)=0; Mst(7,6,i)=22*L(i); Mst(7,7,i)=4*L(i)^2;

Mst(8,5,i)=-22*L(i); Mst(8,6,i)=0 Mst(8,7,i)=0

Mst(8,8,i)=4*L(i)^2;

end

for i=1:1:NN0

Msr(1,1,i)=36;

Msr(2,1,i)=0; Msr(2,2,i)=36;

Msr(3,1,i)=0 Msr(3,2,i)=-3*L(i); Msr(3,3,i)=4*L(i)^2;

Msr(4,1,i)=3*L(i); Msr(4,2,i)=0; Msr(4,3,i)=0;

Msr(4,4,i)=4*L(i)^2; Msr(5,1,i)=36; Msr(5,2,i)=0;

Msr(5,3,i)=0; Msr(5,4,i)=-3*L(i); Msr(6,1,i)=0;

Msr(6,2,i)=-36; Msr(6,3,i)= 3*L(i); Msr(6,4,i)=0;

Msr(7,1,i)=0; Msr(7,2,i)=-3*L(i); Msr(7,3,i)=-L(i)^2;

Msr(7,4,i)=0; Msr(8,1,i)=3*L(i); Msr(8,2,i)=0;

Msr(8,3,i)=0; Msr(8,4,i)=-L(i)^2; Msr(5,5,i)=36;

Msr(6,5,i)=0; Msr(6,6,i)=36; Msr(7,5,i)=0;

Msr(7,6,i)=3*L(i); Msr(7,7,i)=4*L(i)^2; Msr(8,5,i)=-3*L(i);

Msr(8,6,i)=0; Msr(8,7,i)=0; Msr(8,8,i)=4*L(i)^2; end

for i=1:1:NN0

Ge(1,1,i)=0;

Ge(2,1,i)=36; Ge(2,2,i)=0;

Ge(3,1,i)=-3*l(i); Ge(3,2,i)=0; Ge(3,3,i)=0;

Ge(4,1,i)=0; Ge(4,2,i)=-3*L(i); Ge(4,3,i)=4*L(i)^2;

Ge(4,4,i)=0; Ge(5,1,i)=0; Ge(5,2,i)=36;

Ge(5,3,i)=-3*L(i); Ge(5,4,i)=0; Ge(6,1,i)=-36;

Ge(6,2,i)=0; Ge(6,3,i)=0; Ge(6,4,i)=-3*L(i0);

Ge(7,1,i)=-3*L(i); Ge(7,2,i)=0; Ge(7,3,i)=0;

Ge(7,4,i)=L(i)62; Ge(8,1,i)=0; Ge(8,2,i)=-3*L(i);

Ge(8,3,i)=-L(i)^2; Ge(8,4,i)=0; Ge(5,5,i)=0;

Ge(6,5,i)=36; Ge(6,6,i)=0; Ge(7,5,i)=3*L(i);

Ge(7,6,i)=0; Ge(7,7,i)=0; Ge(8,5,i)=0;

Ge(8,6,i)=3*L(i); Ge(8,7,i)=4*L(i)^2; Ge(8,8,i)=0;

end

for i=1:1:NN0

Ks(1,1,i)=12;

Ks(2,1,i)=0; Ks(2,2,i)=12;

Ks(3,1,i)=0; Ks(3,2,i)=-6*L(i); Ks(3,3,i)=4*L(i)^2;

Ks(4,1,i)=6*L(i); Ks(4,2,i)=0; Ks(4,3,i)=0;

Ks(4,4,i)=4*L(i)^2; Ks(5,1,i)=-12; Ks(5,2,i)=0;

Ks(5,3,i)=0; Ks(5,4,i)=-6*L(i); Ks(6,1,i)=0;

Ks(6,2,i)=-12; Ks(6,3,i)=6*L(i); Ks(6.4.i)=0;

Ks(7,1,i)=0; Ks(7,2,i)=-6*L(i); Ks(7,3,i)=2*L(i)^2;

Ks(7,4,i)=0; Ks(8.1,i)=6*L(i); Ks(8,2,i)=0;

Ks(8,3,i)=0; Ks(8,4,i)=2*L(i)62; Ks(5,5,i)=12;

Ks(6,5,i)=0; Ks(6,6,i)=12; Ks(7,5,i)=0;

Ks(7,6,i)=6*L(i); Ks(7,7,i)=4*L(i)^2; Ks(8,5,i)=-6*L(i);

Ks(8,6,i)=0; Ks(8,7,i)=0; Ks(8,8,i)=4*L(i)^2; end

for i=1:1:NN0

for j=1:1:8

for k=1:1:8

EI=Ef(i)*pi*(R(i)^4-R0(i)^4-)/4;

Mst(j,k,i)=Mst(j,k,i)*miu(i)*L(i)/420;

Msr(j,k,i)=Msr(j,k,i)*miu(i)*R(i)^2/120/L(i);

Ge(j,k,i)=-Ger(j,k,i)*2*miu(i)*R(i)^2/120/L(i);

Ks(j,k,i)=Ks(j,k,i)*EI/L(i)^3;

end

end

end

for i=1:1:NN0

for j=1:1:8

for k=j:1:8

Mst(j,k,i)=Mst(k,j,i);

Msr(j,k,i)=Msr(k,j,i);

Ks(j,k,i)=Ks(k,j,i);

Ge(j,k,i)=-Ge(k,j,i);

end

end

end

%M_G_K.m

%该程序进行距阵组集

function [M G K]=M_G_K(N,Ef,R,R0,Mst,Msr,Ge,Ks,miu,L)

% M: 总的质量距阵

% K：总的刚度距阵

% G：总的陀螺力距距阵

NN0=N

for 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(j,k,i)=Ks(j,k,i);

Ge(j,k,i)=-Ge(j,k,i);

end

end

eng

for i=1:1:N

for j=1:1:8

for K=1:1:8

M(i*4+j-4,i*4+k-4)=Ms(j,k,i);

G(i*4+j-4,i*4+k-4)=Ge(j,k,i);

K(i*4+j-4,i*4+k-4)= K (j,k,i);

end

end

end

for i=2:1:N

for k=1:1:4

M(i*4+j-4,i*4+k-4)= M(i*4+j-4,i*4+k-4)+Ms(j+4,k+4,i-1);

G(i*4+j-4,i*4+k-4)= G(i*4+j-4,i*4+k-4)+Ge(j+4,k+4,i-1);

K(i*4+j-4,i*4+k-4)= K(i*4+j-4,i*4+k-4)+Ks(j+4,k+4,i-1);

end

end

end

Kzx(1)=7e7;

Kzy(1)=7e7;

Kzx(12)=7e7;

Kzy(12)=7e7;

for i=1:1:N+1

K(i*4+1-4,i*4+1-4)= K(i*4+1-4,i*4+1-4)+Kzx(i);

K(i*4+2-4,i*4+2-4)= K(i*4+2-4,i*4+2-4)+Kzy(i);

end

%K_D.m

%该程序将组尼加入总体距阵

function [K C D Ax]=K_D(N,K,M,G)

Kzx(1)=7e7;

Kzy(1)=7e7;

Kzx(12)=7e7;

Kzy(12)=7e7;

for i=1:1:N+1

K(i*4+1-4,i*4+1-4)= K(i*4+1-4,i*4+1-4)+Kzx(i);

K(i*4+2-4,i*4+2-4)= K(i*4+2-4,i*4+2-4)+Kzy(i);

end

format long g

[Ax WW]=eig(inv(M)*K);

f=sqrt(WW)/(2*pi);

f0=diag(f);

f00=abs(sort(f0))

fn123=[f00(1) f00(3) f00(5)]

wi1=54*2*pi;

wi2=284*2*pi;

yita1=0.1;

yita2=0.2;

alf=2*(yita2/wi2-yita1/wi1)*(1/wi1^2-1/wi1^2);

beita=2*(yita2*wi2-yita1*wi1)/(1/wi1^2-wi1^2);

C=alf*M+beita*K;

D=C+G;

Dzx(1)=7e3;

Dzy(1)=7e3;

Dzx(12)=7e3;

for i=1:1:N+1

D(i*4+1-4,i*4+1-4)=D(i*4+1-4,i*4+1-4)+Dzx(i);

D(i*4+2-4,i*4+2-4)=D(i*4+2-4,i*4+2-4)+Dzy(i);

end

D=D*1;

%newmark_newton_multi.m

%该程序为newmark_newton算法

xn=yn(:,n);

dxn=dyn(:,n);

ddxn=ddyn(:,n);

if i==1

xn=yn(:,n);

end

if i>1

xn=yn(:,n+1);

end

K0=K;

K1=K*0;

F(loca*4-4+1,1)=Famp1*cos(wi*t+fai(1));

F(loca*4-4+2,1)=Famp1*sin(wi*t+fai(1));

F(loca*4+1,1)=Famp1*cos(wi*t+fai(1));

F(loca*4+2,1)=Famp1*sin(wi*t+fai(1));

Pn1=F;

Fr=Pn1*0;

dert=20/1e6;

k_rub=5e4;

miu_rub=0.2;

xx=yn(loc_rub*4-4+1,n)+5/1e6;

yy=yn(loc_rub*4-4+2,n);

rr=sqrt(xx^2+yy^2);

if n>Fen*60&rr>dert

rub_sign=1;

K1(loc_rub*4-4+1,loc_rub*4-4+1)=k_rub*(rr-dert);

K1(loc_rub*4-4+2,loc_rub*4-4+1)=k_rub*(rr-dert)*miu_rub;

K1(loc_rub*4-4+1,loc_rub*4-4+2)=-k_rub*(rr-dert)*miu_rub;

K1(loc_rub*4-4+2,loc_rub*4-4+2)=k_rub*(rr-dert);

Fr(loc_rub*4-4+1,l)=k_rub*(rr-dert)/rr*(xx-miu_rub*yy);

Fr(loc_rub*4-4+2,l)=k_rub*(rr-dert)/rr*(yy+miu_rub*xx); end

KT=K0+K1;

If i==1

Fn1=K0*xn+Fr;

end

if i>1

Fn1=K0*xn1+Fr;

end

Kj=KT+a0*M+a1*C;

Pj=Pn1+M*(a0*xn+a2*dxn+a3*ddxn)+C*(a1*xn+a4*dxn+a5*ddxn); If i==1

Pj=Pj-Fn1+KT*xn;

end

if i>1

Pj=Pj-Fn1+KT*xn1;

end

if rr

i=100;

end

[QQ,RR]=qr(Kj);

xn1=RR\QQ’ *Pj;

ddxn1=a0*(xn1-xn)-a2*dxn-a3*ddxn;

dxn1=dxn+a6*ddxn+a7*ddxn1;

yn(:,n+1)=xn1;

dyn(:,n+1)=dxn1;

ddyn(;,n+1)=ddxn1;

A.3 稳定性分析程序

% 第一步：设置初始条件

[N Fen y d y2]=initial_conditions

fm=[];

for kk=1:1:60

w=40+(kk-1)*1/1 %频率区间

h=1/w/Fen; %每个周期内积分点数为Fen=200 w=2*pi*w;

% 第二步：通过打靶法计算振幅的分岔特性

shooting;

% 第三步：计算floquet稳定性并保存计算数据

floquet_stability

end

%shooting.m

%该程序主要通过打靶法计算振幅的分岔特征并存储计算数据

for i=1:1:30*Fen

t=0;

t=t+(i)*h;

y=rkutta(t,h,y,w);

if i>(20*Fen) & mod(i,Fen)==1

fprintf(f1,w/2/pi,y(1),y(2),y(3),y(4), y(5),y(6),y(7),y(8));

end

end

%floquet_stability.m

%该程序主要计算floquet稳定性并保存计算数据

y2=eye(N);

s0=y;

for i=1:1:Fen*1

t=0;

t=t+(i)*h;

y=rkutta(t,h,y,w);

At=fun_At(t,y,w); % 根据当前y值求A

For j=1:1:N % 由dy2=A*y2积分求得y2

zzzz=rkutta21(t,h,y2(:,j)’,At);

end

end

s=y;

Rsi=s0-s;

S=y2;

Drds=eye(N)*1-S;

dsi=Drds\(-Rsi)’ ;

y=s-dsi’ ;

floquet_mul=max(abs(eig(y2)));

fm=[fm;w/2/pi floquet_mul];

fprintf(f2,’ %15.9f, %15.9f\n’ ,w/2/pi, floquet_mul);

% rkutta.m

function yout=rkutta(t,h,y,w)

N=length(y);

for i=1:1:N

a(i)=0;

d(i)=0;

b(i)=0;

y0(i)=0;

end

a(1)=h/2; a(2)=h/2;

a(3)=h; a(4)=h;

d=fun(t,y,w);

b=y;

y0=y;

for k=1:1:3

for i=1:1:N

y(i)=y0(i)+a(k)*d(i);

b(i)=b(i)+a(k+1)*d(i)/3;

end

tt=t+a(k);

d=fun(tt,y,w);

end

for i=1:1:N

yout(i)=b(i)+h*d(i)/6;

end

%rkutta21.m

function yout=ekutta21(t,h,y,A)

N=length(y);

For i=1:1:N

a(i)=0;

d(i)=0;

b(i)=0;

y0(i)=0;

end

a(1)=h/2; a(2)=h/2;

a(3)=h; a(4)=h;

d=fun21(t,y,A);

b=y;

y0=y;

for k=1:1:3

for i=1:1:N

y(i)=y0(i)+a(k)*d(i);

b(i)=b(i)+a(k+1)*d(i)/3;

end

tt=t+a(k);

d=fun21(t,y,A);

end

for i=1:1:N

yout(i)=b(i)+h*d(i)/6;

end

%fun.m

function d=fun(t,y,w)

N=length(y); %圆盘1的质量m1=0.7902; %圆盘2的质量m2=0.7902; %转轴的直径

d0=0.01;

l1=0.16;

l2=0.16;

l3=0.16;

e1=5.695e-3*0; %圆盘1的偏心e2=5.695e-5; %圆盘2的偏心E=2.11e11; %转轴的弹性模量

delta1=200e-5; %碰摩间隙1

delta2=2e-5; %碰摩间隙2

kr1=1e5; %碰摩刚度1

kr2=1e5; %碰摩刚度2

fr1=0.1; %碰摩摩擦系数1

fr2=0.1; %碰摩摩擦系数2

I=pi*d0*d0*d0*d0/64;

k11=2.3704e5;

k22=2.3704e5;

k12=-2.0741e5;

k21=k12; % 刚度

c11=53.118385200000006; c12=-37.333800000000004;

c21=-37.333800000000004; c22=53.118385200000006; % //阻尼系数omega=w;

g=0;

% 转子动力学方程

If y(5)>=deltal

Dirac1=1;

end

if y(5)

Dirac1=0;

end

if y(7)>delta2

Dirac2=1

end

if y(7)>delta2

Dirac2=0;

end

fx1=-kr1*(y(5)-delta1)*Dirac1;

fx2=-kr2*(y(7)-delta2)*Dirac2;

fy1=-fr1*kr1*(y(5)-delta1)*Dirac1;

fy2=-fr2*kr2*(y(7)-delta2)*Dirac2;

d(1)=y(2);

d(2)=-(c11/m1)*y(2)-(c12/m1)*y(4)-(k11/m1)*y(1)-(k12/m1)*y(3) +e1*omega*omega*sin(omega*t)+fy1/m1-g;

d(3)=y(4);

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/m2-g;

d(5)=y(6);

d(6)=-(c11/m1)*y(6)-(c12/m1)*y(8)-(k11/m1)*y(5)-(k12/m1)*y(7) +e1*omega*omega*cos(omega*t)+fy1/m1;

d(7)=y(8);

d(8)=-(c21/m2)*y(6)-(c22/m2)*y(8)-(k21/m2)*y(5)-(k22/m2)*y(7 +e2*omega*omega*cos(omega*t)+fx2/m2;

function At=fun_At(t,y,w)

N=length(y);

eps=1e-6;

for i=1:1:N

yi=y;

yi(i)=y(i)+y(i)*eps;

At0=fun(t,y,w);

Ati=fun(t,yi,w);

For j=1:1:N

At(j,i)=(Ati(j)-At0(j))/(yi(i)*eps);

end

end

A.4 不对中转子系统仿真程序

%如下是考虑轴承不对中产生的附加载荷。

alfa=pi/12;

bita=pi/15;

cc=4*cos(alfa)/(3+cos(2*alfa));

dd=(1-cos(2*alfa))/ (3+cos(2*alfa));

T=(II*wi^2/cos(alfa))*(2*cc*dd*sin(2*wi*t)/(1+dd*cos(2*wi*t)));

F (local*4-4+4,1)=T*sin(alfa)*cos(bita); % 不对中力

F (local*4-4+3,1)=T*sin(alfa)*sin(bita);

A.5 转子系统不对中故障的振动信号小波包分解程序

%WPDmisalignment.m

%该程序主要研究应用基于小波包的能量比例计算方法提取其故障特征

%读取数据

x=load(‘misalignment25hz.txt’,’r+’);

x1=x(1:13312,2);

y1= x(1:13312,3);

dt=0.00078;

fs=1/dt;

N=length(x1);

T=dt*[0:N-1]; %时间轴

% 时城波形

plot(T,x1);

xlabel(‘Time/sec’,’fontsize’,8);

ylabel(‘Disp. /um’,’fontsize’,8);

figure(2)

plot(x1,y1,);

xlabel(’Displacement x/um’,’fontsize’,8); ylabel (’ Displacement y/um’,’fontsize’,8);

%对原信号进行重采样

%re-sampling from 250Hz into 128Hz

fs1=1024;

dt1=1/fs1;

x11=1:N;

x2=1:fs/fs1:N;

x22=interp1(x11,x1,x2,’spline’); Ns=length(x22); T=dt1*[0:Ns-1];

figure(2)

plot(T,x22)

%小波包分解

freq=fs1/(Ns)*(0:Ns/2-1);

t=wpdec(x22,4,’db44’);

%基于小波包分解的能量比例计算

s(1,:)=wprocef(t,[4 0]);

p(1,:)=(abs(fft(s(1,:),freq))).^2;

s(2,:)=wprocef(t,[4 1]);

p(2,:)=(abs(fft(s(2,:),freq))).^2;

s(4,:)=wprocef(t,[4 2]);

p(4,:)=(abs(fft(s(4,:),freq))).^2;

s(3,:)=wprocef(t,[4 3]);

p(3,:)=(abs(fft(s(3,:),freq))).^2;

s(7,:)=wprocef(t,[4 4]);

p(7,:)=(abs(fft(s(7,:),freq))).^2;

s(8,:)=wprocef(t,[4 5]);

p(8,:)=(abs(fft(s(8,:),freq))).^2;

s(6,:)=wprocef(t,[4 6]);

p(6,:)=(abs(fft(s(6,:),freq))).^2;

s(5,:)=wprocef(t,[4 7]);

p(5,:)=(abs(fft(s(5,:),freq))).^2;

s(13,:)=wprocef(t,[4 8]);

p(13,:)=(abs(fft(s(13,:),freq))).^2;

s(14,:)=wprocef(t,[4 9]);

p(14,:)=(abs(fft(s(14,:),freq))).^2;

s(16,:)=wprocef(t,[4 10]);

p(16,:)=(abs(fft(s(16,:),freq))).^2;

s(15,:)=wprocef(t,[4 11]);

p(15,:)=(abs(fft(s(15,:),freq))).^2;

s(11,:)=wprocef(t,[4 12]);

p(11,:)=(abs(fft(s(11,:),freq))).^2;

s(12,:)=wprocef(t,[4 13]);

p(12,:)=(abs(fft(s(12,:),freq))).^2;

s(10,:)=wprocef(t,[4 14]);

p(10,:)=(abs(fft(s(10,:),freq))).^2;

s(9,:)=wprocef(t,[4 15]);

p(9,:)=(abs(fft(s(9,:),freq))).^2;

ps=sum(p);

pp=p/ps

ttt=[1:16]*32;

ttt=[1:16];

bar(ttt,pp,’r’);

ylabel(‘能量比’,’fontsize’, 10, ‘fontweight’, ‘bold’);

xlabel(‘频带’, ‘fontsize’, 10, ‘fontweight’, ‘bold’);

A.6 HHT变换用于分析碰摩转子振动信号的程序

%HHTrubbing.m

%读数据

x=load(‘healthy.txt’, ‘r+’);

X=(x/78.7)*1000;

x=X(:,2);

y=X(:,3);

dt=0.00039;

fs=1/dt;

N=length(x);

t=dt*[0:N-1]; %时间轴

% EMD分解

[c,rn,nIMF]=emd(x’);

figure(4)

title(‘IMFs and residual’)

for i=1;nIMF

subplot(nIMF+1,1,i)

plot(t,c(:,i));axis tight;

ylabel (I,’fontsize’,8)

end

subplot(nIMF+1,1,Nimf+1)

plot(t,rn);

axis tight;

ylabel (‘rn’, ‘fontsize’,8);

xlabel (’Time/sec’,’fontsize’, 8’);

%瞬时频率谱

Fai=[];

w=[];

for i=1:nIMF

Hc1=imag(Hilbert(c(:,i)));

Fai1=atan(Hc1./c(:,i));

Fai=[Fai Fai1];

w1=[diff(Fai1);0];

w=[w,w1];

end

figure(5)

for i=1:nIMF

subplot(nIMF,1,i)

plot(t,w(:,i))

ylabel (I, ‘fontsize’, 8);

end

xlabel (‘Time/sec’, ‘fontsize’, 8);

%Hilbert时频谱

c=c’ ;

[A,fa,tt]=hhspectrum(c);

[im,tt]=toimage(A,fa);

disp_hhs(im,[],fs);

%Hilbert边际谱

NN=size(im,1);

for k=1:NN

bjp(k)=sum(im(k,:))*1/fs;

end

f=(0:NN-1)/NN*(fs/2);

figure(7)

plot(f,bjp);

xlabel(‘Frenquency/Hz’, ‘Fontsize’, 8); ylabel(‘Amplitude/um’, ‘Fontsize’, 8);

% emd.m

%该函数用来完成EMD分解

function[c,rn,Nimf]=emd(y,Ns,fs)

s1=y; % 数据

Ns=length(y);

fs=6.4*1000*2.56;

s2=s1;

for i=1:Ns

s2(Ns-i+1)=s1(i);

end

s=[s2 s1 s2]; %原始信号左右延拓