传递矩阵matlab程序.docx
《传递矩阵matlab程序.docx》由会员分享,可在线阅读,更多相关《传递矩阵matlab程序.docx(38页珍藏版)》请在冰豆网上搜索。
传递矩阵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
fori=1:
1:
M
ptx(i)=0;
pty(i)=0;
end
forii=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’);
gridon
%第四步:
用二分法求固有频率及振型图
%固有频率:
Critical_speed
wi=50;
fori=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;
whilekkk<5000
ifD2*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
ifD1*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;
Ifatep<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
%===============
fori=1:
1:
N+1
forj=1:
1:
2
fork=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
fori=1:
1:
N
forj=1:
1:
2
fork=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
fori=1:
1:
N
forj=1:
1:
2
fork=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)计算质点上传递矩阵―――点矩阵的一部分!
%============
fori=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)计算质点上传递矩阵―――点矩阵的一部分!
%============
fori=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中会把错误的覆盖掉
%============
fori=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
%============
%此处全为计算中间量
%============
fori=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
fori=1:
1:
2
forj=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)形成最终传递矩阵
%============
%Ud11Ud12Ud21Ud22为最终参与计算的传递矩阵
fori=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);
fori=1:
1:
2
forj=1:
1:
2
SS1(i,j)=0;
end
end
fori=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);
fori=1:
1:
N+1
D1=D1*sign(det(Su(:
:
i));%消奇点
end
%======
(2)不平衡响应值EE======
EE(:
:
n+2)=-inv(SS(:
:
N+2)*PP(:
:
N+2);
fori=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_signlocaloc_rubFampwirFamp1fai]=initial_conditions
%第二步:
设置转子系统的参数值
%N:
划分的轴段数
%density:
轴的密度单位为:
[kg/m^3
%Ef:
轴的弹性模量单位为:
[Pa]
%L:
每个轴段的长度单位为:
[m]
%R:
每个轴段的外半径单位为:
[m]
%ro:
每个轴段的内半径单位为:
[m]
%miu:
每个轴段的单元质量[kg/m]
[NdensityEfLRRomiu]=rotor_parameters
%第三步:
设置移动单元质量距阵,移动单元质量距阵,刚度单元质量距阵和陀螺力距距阵
%Mst:
移动单元质量距阵
%Msr:
移动单元质量距阵
%Ks:
刚度单元距阵
%Ge:
陀螺力距单元距阵
[Mst:
MsrKsGe]=Mst_Msr_Ks_Ge(N,density,R,Ro.L.Ef,miu)
%第四步:
距阵组集
%M:
总的质量距阵
%K:
总的刚度距阵
%G:
总的陀螺力矩距阵
[MGK]=M_G_K(N,Ef,R,Ro.Mst,Msr,Ge,Ks,miu,L)
%第五步:
加入支撑刚度和阻尼
[KCDAx]=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;
fori=1:
1:
N3
F(i,1)=0;
end
fori=1:
1;N3
yn(i:
1)=0;
dyn(i:
1)=0;
ddyn(i:
1)=0;
end
t=0;
forn=1:
1:
Fen*80
t=t+ht;
n;
fori=1:
1:
30
newmark_newton_multi
end
ifmod(n,100)==1
n
end
ifn>Fen*60
fori=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;[tu]];
xt1=[xt1;[tx’]];
yt1=[yt1;[ty’]];
end
end
rub_sign
save’xt1.dat’xt1–ascii;
save’yt1.dat’yt1–ascii;
save’ut1.dat’ut1–ascii;
%initial_conditions.m
%该程序主要进行仿真条件设置
Function[rub_signlocaloc_robFampwirFamp1fai]=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=[3030]/180*pi
%roto_parameters.m
%该程序对Jeffcott转子系统进行参数设置
function[NdensityEfLRR0miu]=rotor_parameters
%N:
划分的轴段数
%density:
轴的密度单位为:
[kg/m^3]
%Ef:
轴的弹性模量单位为:
[Pa]
%L每个轴段的长度单位为:
[m]
%R每个轴段的外半径单位为:
[m]
%R0:
每个轴段的内半径单位为:
[m]
%miu:
每个轴段的单元质量[kg/m]
N=11;
Density=7850;
Ef=[2.12.12.12.12.12.12.12.12.12.12.1]*lell;
L=[90.590.580.562.5302022.562.590.590.590.5]/1000;
R=[2020202020902020202020]/2000;
R0=[00000000000]/2000;
fori=1:
1;N
miu(i)=density*pi*(R(i)^2-R0(i)^2)
end
%Mst_Msr_Ks_Ge.m
%该程序设置单元距阵
Function[MstMsrKsGe]=Mst_Msr_Ks_Ge(N,density,R,R0,L,Ef,miu)
%Mst:
移动单元质量距阵
%Msr:
转动单元质量距阵
%Ks:
刚度单元距阵
%Ge:
陀螺力距单元距阵
NN0=N;
NN1=NN0+1
NN2=NN1+1
fori=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)=0Mst(8,7,i)=0
Mst(8,8,i)=4*L(i)^2;
end
fori=1:
1:
NN0
Msr(1,1,i)=36;
Msr(2,1,i)=0;Msr(2,2,i)=36;
Msr(3,1,i)=0Msr(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
fori=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
fori=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
fori=1:
1:
NN0
forj=1:
1:
8
fork=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
fori=1:
1:
NN0
forj=1:
1:
8
fork=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);
G