传递矩阵-matlab程序Word文件下载.doc
《传递矩阵-matlab程序Word文件下载.doc》由会员分享，可在线阅读，更多相关《传递矩阵-matlab程序Word文件下载.doc（23页珍藏版）》请在冰点文库上搜索。
[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<
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;
Critical_speed=wi/2/pi*60
figure;
plot_mode（N,l,SS,Sn）
wi=wi+2;
%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
N+1
forj=1:
2
fork=1:
ud11（j,k.i）=0;
ud12（j,k.i）=0;
ud21（j,k.i）=0;
ud22（j,k.i）=0;
end
end
N
forj=1:
us11（j,k.i）=0;
us12（j,k.i）=0;
us21（j,k.i）=0;
us22（j,k.i）=0;
fork=1:
u11（j,k.i）=0;
u12（j,k.i）=0;
u21（j,k.i）=0;
u22（j,k.i）=0;
%＝＝＝＝＝＝＝＝＝＝＝＝
%（b）计算质点上传递矩阵―――点矩阵的一部分！
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;
%（c）计算质点上传递矩阵―――点矩阵的一部分！
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;
%（d）以下计算的是无质量梁上的传递矩阵―――场矩阵
%计算的锥轴的us是不对的，是随便令的，在后面计算剩余量时，zhui中会把错误的覆盖掉
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;
%此处全为计算中间量
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;
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;
%（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为最终参与计算的传递矩阵
u11（:
:
i）=us11（:
i）*ud11（:
i）+us12（:
i）*ud21（:
i）;
u12（:
i）*ud12（:
i）*ud22（:
u21（:
i）=us21（:
i）+us22（:
u22（:
u11（:
N+1）=ud11（:
N+1）;
u12（:
N+1）=ud12（:
u21（:
N+1）=ud21（:
u22（:
N+1）=ud22（:
SS1（i,j）=0;
ud11=u11（:
ud12=u12（:
ud21=u21（:
ud22=u22（:
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）;
%＝＝＝＝＝＝
（2）计算剩余量＝＝＝＝＝＝
D1=det（SS（:
N+2）;
D1=D1*sign（det（Su（:
i））;
%消奇点
%＝＝＝＝＝＝
（2）不平衡响应值EE＝＝＝＝＝＝
EE（:
n+2）=-inv（SS（:
N+2）*PP（:
fori=N+1:
-1:
1
EE（:
I）=Sn（:
i）*EE（:
i+1）-Sn（:
i）*UF（:
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:
每个轴段的外半径单位为：
%ro:
每个轴段的内半径单位为：
%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;
N3
F（i,1）=0;
1;
yn（i:
1）=0;
dyn（i:
ddyn（i:
t=0;
forn=1:
Fen*80
t=t+ht;
n;
fori=1:
30
newmark_newton_multi
ifmod（n,100）==1
n
end
ifn>
Fen*60
fori=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
u=F（loca*4-4+1,1）;
ut1=[ut1;
[tu]];
xt1=[xt1;
[tx’]];
yt1=[yt1;
[ty’]];
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=0,说明系统无碰摩故障；
否则rub _sign=1
不平衡质量的位置
%loc _rub:
碰摩位置
不平衡质量的大小单位为：
转速单位为：
%r：
偏心半径单位为:
%Famp1：
离心力的大小单位为:
不平衡量的初始相位[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:
划分的轴段数
轴的密度单位为：
[kg/m^3]
轴的弹性模量单位为:
%L每个轴段的长度单位为:
%R每个轴段的外半径单位为:
%R0:
每个轴段的内半径单位为:
每个轴段的单元质量[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;
miu（i）=density*pi*（R（i）^2-R0（i）^2）
%Mst_Msr_Ks_Ge.m
%该程序设置单元距阵
Function[MstMsrKsGe]=Mst_Msr_Ks_Ge（N,density,R,R0,L,Ef,miu）
转动单元质量距阵
NN0=N;
NN1=NN0+1
NN2=NN1+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;
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;
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;
升级会员 