%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 R0 miu]=rotor_parameters;
%第三步：设置移动单元质量距阵，移动单元质量距阵，刚度单元质量距阵和陀螺力距距阵
%Mst: 移动单元质量距阵
%Msr: 移动单元质量距阵
%Ks: 刚度单元距阵
%Ge: 陀螺力距单元距阵
[Mst Msr Ks Ge]=Mst_Msr_Ks_Ge(N,density,R,R0,L,Ef,miu);
%第四步：距阵组集
%M： 总的质量距阵
%K： 总的刚度距阵
%G: 总的陀螺力矩距阵
[M G K]=M_G_K(N,Ef,R,R0,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*ht);
a1=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)*1e6;
y(i,1)=yn(i*4-2,n)*1e6;
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_rub 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;
***************************************************************************
%rotor_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]*1e11;
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(i);
Ge(7,1,i)=-3*L(i); Ge(7,2,i)=0; Ge(7,3,i)=0;
Ge(7,4,i)=L(i)^2; 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)^2; 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)=-Ge(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
end
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
