传递矩阵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
升级会员 