传递矩阵-matlab程序Word文档下载推荐.doc
《传递矩阵-matlab程序Word文档下载推荐.doc》由会员分享，可在线阅读，更多相关《传递矩阵-matlab程序Word文档下载推荐.doc（23页珍藏版）》请在冰豆网上搜索。
%剩余量：
D1
fori=1:
1:
M
ptx（i）=0;
pty（i）=0;
end
forii=1:
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
ylabel（‘剩余量’）；
plot（ptx,pty）
xlabel（‘角速度red/s’）;
gridon
%第四步：
用二分法求固有频率及振型图
%固有频率：
Critical_speed
wi=50;
fori=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<
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_rubFam
升级会员 