求助:求传递矩阵法计算临界转速的MATLAB程序
各位朋友,小弟做毕设急需"用传递矩阵法计算临界转速的MATLAB程序",对象是内燃机的主轴系统(单转子系统)
希望各位朋友大哥大姐帮帮忙,小弟在此不胜感激
QQ:383778438 EMAIL: blandise@163.com 自己写吧，不难，有半天时间就差不多了 原帖由 blandise 于 2007-5-14 17:15 发表
各位朋友,小弟做毕设急需"用传递矩阵法计算临界转速的MATLAB程序",
对象是内燃机的主轴系统(单转子系统)
希望各位朋友大哥大姐帮帮忙,小弟在此不胜感激
QQ:383778438 EMAIL: blandise@163.com
我这有个计算的程序，你拿去参考一下吧，你要改参数和模型的
%求解转子系统前三个临界转速和主振型的传递矩阵法
clc
clear
%等截面轴参数
l1=０.１５;
d=0.04;
A=pi*d*d/4;
%轮盘参数
D=0.5;
h=0.025;
%盘轴材料参数（忽略轴的质量）
a=1;
u=0.3;
rou=7800;
E=2.0e11;
G=E/(2*(1+u));
I=pi*(d^4)/64;
K1=2.0e7;
v1=6*E*I/(a*G*A*l1*l1);
mi=rou*pi*D^2/4*h;%轮盘的集质量
Jp=mi*D^2/8;
Jd=Jp/2;
Ji=Jp-Jd;
%参数的数组形式
L=;
M=;
K=;
v=;
J=;
k=0;
Tit=['第一阶频率的振型和弯矩图';'第二阶频率的振型和弯矩图';'第三阶频率的振型和弯矩图'];
for w=0:0.01:4000;
for i=1:15;
T(:,:,i)=
(L(i)^2)*(M(i)*w^2-K(i))/(2*E*I) 1+L(i)*J(i)*w^2/(E*I) L(i)/(E*I) L(i)^2/(2*E*I);
L(i)*(M(i)*w^2-K(i)) J(i)*w^2 1 L(i);
M(i)*w^2-K(i) 0 0 1];
end
H=T(:,:,1);
for i2=2:15;
H=T(:,:,i2)*H;
end
F=H(3,1)*H(4,2)-H(3,2)*H(4,1);
if F*(-1)^k < 0 %求解临界转速
k=k+1;
wi(k)=w;
w=wi(k)
ni(k)=wi(k)*30/pi;
end
end
for i1=1:3;
w=wi(i1);
for j=1:14;
T(:,:,j)=
(L(j)^2)*(M(j)*w^2-K(j))/(2*E*I) 1+L(j)*J(j)*w^2/(E*I) L(j)/(E*I) L(j)^2/(2*E*I);
L(j)*(M(j)*w^2-K(j)) J(j)*w^2 1 L(j);
M(j)*w^2-K(j) 0 0 1];
end
H=T(:,:,1);
for j=2:15;
H=T(:,:,j)*H;
end
b=-H(4,1)/H(4,2);
X(:,1)=(');
for n=2:16;
X(:,n)=T(:,:,n-1)*X(:,n-1);
%相邻两质点右边的传递关系
end
for j1=1:15;
y(j1)=X(1,j1);
z(j1)=X(3,j1);
x(j1)=(j1-1)*l1;
end
y(16)=X(1,16);
x(16)=1.56;
z(16)=X(3,16);
y=y/max(abs(y));%归一化
z=z/max(abs(z));
subplot(3,1,i1)
plot(x,y,'b-',x,z,'r:')
title(Tit(i1,:))
xlabel('轴长'),ylabel('不平衡值')
axis()
grid on
z;
end
legend('振型','弯矩')
ni
wi
回复 #3 hector1982 的帖子
楼上是一强人，不顶对不起天下！ thanks!! 佩服啊，顶一个ba 这个论坛卧虎藏龙
是个值得期待的论坛
望大家交流，共同进步:@) 确实厉害我也需要 谢谢，真是高手啊 这个程序需要设定什么参数呢？ 原帖由 leqinwei 于 2009-2-19 20:21 发表
这个程序需要设定什么参数呢？
从程序上没有，只要你给出转子结构参数和物性参数就可以 确实是好东西，振动论坛会越来越火的。 是很强，偶学习学习。 clc
clear
%等截面轴参数
ll1=0.08;ll2=0.079;ll3=0.088;ll4=0.053;ll5=0.1;ll6=0.02;ll7=0.02;ll8=0.08;ll9=0.08;ll10=0.02;ll11=0.02;ll12=0.1;ll13=0.053;ll14=0.088;ll15=0.081;ll16=0.066;
lg1=0.03;lg2=0.045;lg3=0.045;lg4=0.03;lg5=0.045;lg6=0.08;lg7=0.08;lg8=0.1;lg9=0.076;lg10=0.122;lg11=0.122;lg12=0.122;lg13=0.122;lg14=0.085;lg15=0.085;
lh1=0.058;lh2=0.066;lh3=0.08;lh4=0.053;lh5=0.1;lh6=0.02;lh7=0.02;lh8=0.08;lh9=0.08;lh10=0.02;lh11=0.02;lh12=0.1;lh13=0.035;
dl1=0.08;dl2=dl1;dl3=dl1;dl4=dl1;dl5=dl1;dl6=0.133;dl7=0.113;dl8=0.143;dl9=0.143;dl10=0.133;dl11=0.133;dl12=dl1;dl13=dl1;dl14=dl1;dl15=dl1;dl16=dl1;
dg1=0.2;dg2=0.125;dg3=0.125;dg4=0.2;dg5=0.225;dg6=0.2;dg7=0.2;dg8=0.125;dg9=0.115;dg10=dg9;dg11=dg9;dg12=dg9;dg13=dg9;dg14=dg9;dg15=dg9;
dh1=0.08;dh2=dh1;dh3=dh1;dh4=dh1;dh5=dh1;dh6=0.128;dh7=0.128;dh8=0.132;dh9=0.132;dh10=0.128;dh11=0.128;dh12=dh1;dh13=dh1;
Al1=pi*dl1^2/4;Al2=pi*dl2^2/4;Al3=pi*dl3^2/4;Al4=pi*dl4^2/4;Al5=pi*dl5^2/4;Al6=pi*dl6^2/4;Al7=pi*dl7^2/4;Al8=pi*dl8^2/4;Al9=pi*dl9^2/4;Al10=pi*dl10^2/4;Al11=pi*dl11^2/4;Al12=pi*dl12^2/4;Al13=pi*dl13^2/4;Al14=pi*dl14^2/4;Al15=pi*dl15^2/4;Al16=pi*dl16^2/4;
Ag1=pi*dg1^2/4;Ag2=pi*dg2^2/4;Ag3=pi*dg3^2/4;Ag4=pi*dg4^2/4;Ag5=pi*dg5^2/4;Ag6=pi*dg6^2/4;Ag7=pi*dg7^2/4;Ag8=pi*dg8^2/4;Ag9=pi*dg9^2/4;Ag10=pi*dg10^2/4;Ag11=pi*dg11^2/4;Ag12=pi*dg12^2/4;Ag13=pi*dg13^2/4;Ag14=pi*dg14^2/4;Ag15=pi*dg15^2/4;
Ah1=pi*dh1^2/4;Ah2=pi*dh2^2/4;Ah3=pi*dh3^2/4;Ah4=pi*dh4^2/4;Ah5=pi*dh5^2/4;Ah6=pi*dh6^2/4;Ah7=pi*dh7^2/4;Ah8=pi*dh8^2/4;Ah9=pi*dh9^2/4;Ah10=pi*dh10^2/4;Ah11=pi*dh11^2/4;Ah12=pi*dh12^2/4;Ah13=pi*dh13^2/4;
%盘轴材料参数（忽略轴的质量）
a=0.886;
Kc=1.96e6;%油膜刚度
Kb=2.7048e9;%轴承座刚度
mb=80;%轴承座参振质量
u=0.3;
E=2.068e11;
G=E/(2*(1+u));
Il1=pi*(dl1^4)/64;Il2=pi*(dl2^4)/64;Il3=pi*(dl3^4)/64;Il4=pi*(dl4^4)/64;Il5=pi*(dl5^4)/64;Il6=pi*(dl6^4)/64;Il7=pi*(dl7^4)/64;Il8=pi*(dl8^4)/64;Il9=pi*(dl9^4)/64;Il10=pi*(dl10^4)/64;Il11=pi*(dl11^4)/64;Il12=pi*(dl12^4)/64;Il13=pi*(dl13^4)/64;Il14=pi*(dl14^4)/64;Il15=pi*(dl15^4)/64;Il16=pi*(dl16^4)/64;
Ig1=pi*(dg1^4)/64;Ig2=pi*(dg2^4)/64;Ig3=pi*(dg3^4)/64;Ig4=pi*(dg4^4)/64;Ig5=pi*(dg5^4)/64;Ig6=pi*(dg6^4)/64;Ig7=pi*(dg7^4)/64;Ig8=pi*(dg8^4)/64;Ig9=pi*(dg9^4)/64;Ig10=pi*(dg10^4)/64;Ig11=pi*(dg11^4)/64;Ig12=pi*(dg12^4)/64;Ig13=pi*(dg13^4)/64;Ig14=pi*(dg14^4)/64;Ig15=pi*(dg15^4)/64;
Ih1=pi*(dh1^4)/64;Ih2=pi*(dh2^4)/64;Ih3=pi*(dh3^4)/64;Ih4=pi*(dh4^4)/64;Ih5=pi*(dh5^4)/64;Ih6=pi*(dh6^4)/64;Ih7=pi*(dh7^4)/64;Ih8=pi*(dh8^4)/64;Ih9=pi*(dh9^4)/64;Ih10=pi*(dh10^4)/64;Ih11=pi*(dh11^4)/64;Ih12=pi*(dh12^4)/64;Ih13=pi*(dh13^4)/64;
%齿轮参数
zg=226;zl=43;zh=39;b=0.153;%齿宽
mg=589.52;ml=21.3381;mh=17.5529;%齿轮重量
m=3.5;%齿轮模数
Rl=m*zl/2/1000;%L轴上齿轮半径（单位m）
Rh=m*zh/2/1000;%H轴上齿轮半径（单位m）
Rg=m*zg/2/1000;%G轴上齿轮半径（单位m）
znb=0.05;%阻尼比
e=1.3;%重合度
kgl=1/(0.04723+0.1551/zg+0.25791/zl)*(0.75*e+0.25)*b*1e6;%齿轮GL总啮合刚度
cgl=2*znb*sqrt(kgl/(1/mg+1/ml));%齿轮GL有效啮合阻尼
kgh=1/(0.04723+0.1551/zg+0.25791/zh)*(0.75*e+0.25)*b*1e6;%齿轮GH总啮合刚度
cgh=2*znb*sqrt(kgh/(1/mg+1/mh));%齿轮GH有效啮合阻尼
vl1=6*E*Il1/(a*G*Al1*ll1^2);vl2=6*E*Il2/(a*G*Al2*ll2^2);vl3=6*E*Il3/(a*G*Al3*ll3^2);vl4=6*E*Il4/(a*G*Al4*ll4^2);vl5=6*E*Il5/(a*G*Al5*ll5^2);vl6=6*E*Il6/(a*G*Al6*ll6^2);vl7=6*E*Il7/(a*G*Al7*ll7^2);vl8=6*E*Il8/(a*G*Al8*ll8^2);vl9=6*E*Il9/(a*G*Al9*ll9^2);vl10=6*E*Il10/(a*G*Al10*ll10^2);vl11=6*E*Il11/(a*G*Al11*ll11^2);vl12=6*E*Il12/(a*G*Al12*ll12^2);vl13=6*E*Il13/(a*G*Al13*ll13^2);vl14=6*E*Il14/(a*G*Al14*ll14^2);vl15=6*E*Il15/(a*G*Al15*ll15^2);vl16=6*E*Il16/(a*G*Al16*ll16^2);
vg1=6*E*Ig1/(a*G*Ag1*lg1^2);vg2=6*E*Ig2/(a*G*Ag2*lg2^2);vg3=6*E*Ig3/(a*G*Ag3*lg3^2);vg4=6*E*Ig4/(a*G*Ag4*lg4^2);vg5=6*E*Ig5/(a*G*Ag5*lg5^2);vg6=6*E*Ig6/(a*G*Ag6*lg6^2);vg7=6*E*Ig7/(a*G*Ag7*lg7^2);vg8=6*E*Ig8/(a*G*Ag8*lg8^2);vg9=6*E*Ig9/(a*G*Ag9*lg9^2);vg10=6*E*Ig10/(a*G*Ag10*lg10^2);vg11=6*E*Ig11/(a*G*Ag11*lg11^2);vg12=6*E*Ig12/(a*G*Ag12*lg12^2);vg13=6*E*Ig13/(a*G*Ag13*lg13^2);vg14=6*E*Ig14/(a*G*Ag14*lg14^2);vg15=6*E*Ig15/(a*G*Ag15*lg15^2);
vh1=6*E*Ih1/(a*G*Ah1*lh1^2);vh2=6*E*Ih2/(a*G*Ah2*lh2^2);vh3=6*E*Ih3/(a*G*Ah3*lh3^2);vh4=6*E*Ih4/(a*G*Ah4*lh4^2);vh5=6*E*Ih5/(a*G*Ah5*lh5^2);vh6=6*E*Ih6/(a*G*Ah6*lh6^2);vh7=6*E*Ih7/(a*G*Ah7*lh7^2);vh8=6*E*Ih8/(a*G*Ah8*lh8^2);vh9=6*E*Ih9/(a*G*Ah9*lh9^2);vh10=6*E*Ih10/(a*G*Ah10*lh10^2);vh11=6*E*Ih11/(a*G*Ah11*lh11^2);vh12=6*E*Ih12/(a*G*Ah12*lh12^2);vh13=6*E*Ih13/(a*G*Ah13*lh13^2);
%轮盘的集质量
ml1=34.76;ml2=3.12;ml3=3.47;ml4=2.09;ml5=3.95;ml6=2.18;ml7=1.57;ml8=431.6;ml9=10.09;ml10=2.18;ml11=2.18;ml12=3.95;ml13=2.09;ml14=3.16;ml15=37.9;ml16=2.6;
mg1=7.4;mg2=4.34;mg3=4.34;mg4=7.4;mg5=14.05;mg6=609.25;mg7=19.73;mg8=0.96;mg9=6.2;mg10=9.95;mg11=9.95;mg12=9.95;mg13=9.95;mg14=6.93;mg15=6.93;
mh1=22.39;mh2=2.6;mh3=3.16;mh4=2.09;mh5=3.95;mh6=2.02;mh7=2.02;mh8=397.67;mh9=8.59;mh10=2.02;mh11=2.02;mh12=3.95;mh13=1.38;
Jpl1=3.3025557e-1; Jpl2=ml2*dl2^2/8; Jpl3=ml3*dl3^2/8; Jpl4=ml4*dl4^2/8; Jpl5=ml5*dl5^2/8; Jpl6=ml6*dl6^2/8; Jpl7=ml7*dl7^2/8; Jpl8=ml8*(2*Rl)^2/8; Jpl9=ml9*dl9^2/8; Jpl10=ml10*dl10^2/8; Jpl11=ml11*dl11^2/8; Jpl12=ml12*dl12^2/8; Jpl13=ml13*dl13^2/8; Jpl14=ml14*dl14^2/8; Jpl15=3.453623e-1; Jpl16=0.444;
Jpg1=mg1*dg1^2/8;Jpg2=mg2*dg2^2/8;Jpg3=mg3*dg3^2/8;Jpg4=mg4*dg4^2/8;Jpg5=mg5*dg5^2/8;Jpg6= 5.027e1;Jpg7=mg7*dg7^2/8;Jpg8=mg8*dg8^2/8;Jpg9=mg9*dg9^2/8;Jpg10=mg10*dg10^2/8;Jpg11=mg11*dg11^2/8;Jpg12=mg12*dg12^2/8;Jpg13=mg13*dg13^2/8;Jpg14=mg14*dg14^2/8;Jpg15=mg15*dg15^2/8;
Jph1=1.7080093e-1;Jph2=mh2*dh2^2/8;Jph3=mh3*dh3^2/8;Jph4=mh4*dh4^2/8;Jph5=mh5*dh5^2/8;Jph6=mh6*dh6^2/8;Jph7=mh7*dh7^2/8;Jph8=mh8*(Rh*2)^2/8;Jph9=mh9*dh9^2/8;Jph10=mh10*dh10^2/8;Jph11=mh11*dh11^2/8;Jph12=mh12*dh12^2/8;Jph13=mh13*dh13^2/8;
Jdl1=2.3812813e-1;Jdl2=Jpl2/2;Jdl3=Jpl3/2;Jdl4=Jpl4/2;Jdl5=Jpl5/2;Jdl6=Jpl6/2;Jdl7=Jpl7/2;Jdl8=Jpl8/2;Jdl9=Jpl9/2;Jdl10=Jpl10/2;Jdl11=Jpl11/2;Jdl12=Jpl12/2;Jdl13=Jpl13/2;Jdl14=Jpl14/2;Jdl15=1.9737902e-1 ;Jdl16=0.258;
Jdg1=Jpg1/2;Jdg2=Jpg2/2;Jdg3=Jpg3/2;Jdg4=Jpg4/2;Jdg5=Jpg5/2;Jdg6=2.64e1;Jdg7=Jpg7/2;Jdg8=Jpg8/2;Jdg9=Jpg9/2;Jdg10=Jpg10/2;Jdg11=Jpg11/2;Jdg12=Jpg12/2;Jdg13=Jpg13/2;Jdg14=Jpg14/2;Jdg15=Jpg15/2;
Jdh1=9.22548e-2;Jdh2=Jph2/2;Jdh3=Jph3/2;Jdh4=Jph4/2;Jdh5=Jph5/2;Jdh6=Jph6/2;Jdh7=Jph7/2;Jdh8=Jph8/2;Jdh9=Jph9/2;Jdh10=Jph10/2;Jdh11=Jph11/2;Jdh12=Jph12/2;Jdh13=Jph13/2;
%参数的数组形式
Ll= ;
Lg=;
Lh=;
vl=;
vg=;
vh=;
Jdl=;
Jdh=;
Jdg=;
Jpl=;
Jph=;
Jpg=;
Ml=;
Mg=;
Mh=;
Il=;
Ig=;
Ih=;
k=0;
for w=0:1:4000;
p=cos(20/180*pi);q=sin(20/180*pi);
K11=Kc*(Kb-mb*w^2)/(Kc+Kb-mb*w^2);%支撑总刚度
Kx=K11*p;Ky=K11*q;%支撑刚度x和y向的分量
Kxxl=;
Kyyl=;
Kxxg=;
Kyyg=;
Kxxh=;
Kyyh=;
for n=1:15;
T(:,:,n)=[1-(Lg(n)^3)*(1-vg(n))*(-Mg(n)*w^2+Kxxg(n))/(6*E*Ig(n)) Lg(n)-Lg(n)^2*Jdg(n)*w^2/(2*E*Ig(n)) 0 Lg(n)^2*Jpg(n)*w^2*i/(2*E*Ig(n))Lg(n)^2/(2*E*Ig(n)) Lg(n)^3*(1-vg(n))/(6*E*Ig(n)) 0 0;
(Lg(n)^2)*(Mg(n)*w^2-Kxxg(n))/(2*E*Ig(n)) 1-Lg(n)*Jdg(n)*w^2/(E*Ig(n)) 0 Lg(n)*Jpg(n)*w^2*i/(E*Ig(n))Lg(n)/(E*Ig(n))Lg(n)^2/(2*E*Ig(n)) 0 0;
0 -Lg(n)^2*Jpg(n)*w^2*i/(2*E*Ig(n)) 1+(Lg(n)^3)*(1-vg(n))*(Mg(n)*w^2-Kyyg(n))/(6*E*Ig(n)) Lg(n)-Lg(n)^2*Jdg(n)*w^2/(2*E*Ig(n)) 0 0 Lg(n)^2/(2*E*Ig(n)) Lg(n)^3*(1-vg(n))/(6*E*Ig(n));
0 -Lg(n)*Jpg(n)*w^2*i/(E*Ig(n))Lg(n)^2*(Mg(n)*w^2-Kyyg(n))/(2*E*Ig(n)) 1-Lg(n)*Jdg(n)*w^2/(E*Ig(n)) 0 0Lg(n)/(E*Ig(n))Lg(n)^2/(2*E*Ig(n));
Lg(n)*(Mg(n)*w^2-Kxxg(n)) -Jdg(n)*w^2 0 i*Jpg(n)*w^2 1 Lg(n) 0 0;
Mg(n)*w^2-Kxxg(n) 0 0 0 0 1 0 0;
0 -i*Jpg(n)*w^2 Lg(n)*(Mg(n)*w^2-Kyyg(n)) -Jdg(n)*w^2 0 0 1 Lg(n);
0 0 Mg(n)*w^2-Kyyg(n) 00 0 0 1] ;
end
H=T(:,:,1);
for n2=2:15;%矩阵传递
H=T(:,:,n2)*H;
end
F=det([H(5,1) H(5,2) H(5,3) H(5,4);
H(6,1) H(6,2) H(6,3) H(6,4);
H(7,1) H(7,2) H(7,3) H(7,4) ;
H(8,1) H(8,2) H(8,3) H(8,4)]);
if F*(-1)^k < 0 %求解临界转速
k=k+1;
wi(k)=w;
w=wi(k);
ni(k)=wi(k)*30/pi;
end
end
wi
ni LS的厉害，不知可否解释一下，传递矩阵状态参量怎么有八个？