用matlab编写的传递矩阵法求主轴的临界转速

zhubeidou · 发表于 2007-6-2 17:05

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

ml=[0,0,0,0,0,0,0,0,0,0,0,0,0];
jl=[0,0,0,0,0,0,0,0,0,0,0,0,0];
eml=[0,0,0,0,0,0,0,0,0,0,0,0,0];
d=[0,0,0,0,0,0,0,0,0,0,0,0,0];
l=[0,0,0,0,0,0,0,0,0,0,0,0,0];
I=[0 0 0 0 0 0 0 0 0 0 0 0 0];
k5=0;
k10=0;
c5=0;
c10=0;
ro=7.85e3;
pi=3.1416;
k5=2e9;
k10=2e9;
c5=1.2e5;
c10=1.2e5;
d=[1.2 1.8 2.0 2.0 2.0 3.0 3.5 3.0 2.0 2.0 2.0 1.8 1.2]/100;
l=[2.0 0.5 0.5 2.2 0.6 1.0 5.0 1.0 0.6 2.2 0.5 0.5 2.0]/100;
for i=1:13
ml(1,i)=ro*pi*d(1,i)*d(1,i)*l(1,i)/4;
end
m=[0,0,0,0,0,0,0,0,0,0,0,0,0,0];
eml=[0.47 0 0 0 0 0 0 0 0 0 0 0 0];
m(1,1)=ml(1,1)/2+eml(1,1)/2;
m(1,14)=ml(1,13)/2+eml(1,1)/2;
for i=1:12
m(1,(i+1))=(ml(1,i)+ml(1,(i+1))+eml(1,i)+eml(1,(i+1)))/2;
end
jl=[0.0009 0 0 0 0 0 0 0 0 0 0 0 0];
j=[0,0,0,0,0,0,0,0,0,0,0,0,0,0];
j(1,1)=jl(1,1)/2;
j(1,14)=jl(1,13)/2;
for i=1:12
j(1,(i+1))=(jl(1,i)+jl(1,(i+1)))/2;
end
for i=1:13
I(1,i)=pi*d(1,i)^4/64;
end
E=200e9;
syms w;
T1=[1 0 0 0;0 1 0 0;0 -j(1,1)*w^2 1 0;m(1,1)*w^2 0 0 1];
T2=[1 0 0 0;0 1 0 0;0 -j(1,2)*w^2 1 0;m(1,2)*w^2 0 0 1]*[1 l(1,1) l(1,1)^2/(2*E*I(1,1)) l(1,1)^3/(6*E*I(1,1));0 1 l(1,1)/(E*I(1,1)) l(1,1)^2/(2*E*I(1,1));0 0 1 l(1,1);0 0 0 1];
T3=[1 0 0 0;0 1 0 0;0 -j(1,3)*w^2 1 0;m(1,3)*w^2 0 0 1]*[1 l(1,2) l(1,2)^2/(2*E*I(1,2)) l(1,2)^3/(6*E*I(1,2));0 1 l(1,2)/(E*I(1,2)) l(1,2)^2/(2*E*I(1,2));0 0 1 l(1,2);0 0 0 1];
T4=[1 0 0 0;0 1 0 0;0 -j(1,4)*w^2 1 0;m(1,4)*w^2 0 0 1]*[1 l(1,3) l(1,3)^2/(2*E*I(1,3)) l(1,3)^3/(6*E*I(1,3));0 1 l(1,3)/(E*I(1,3)) l(1,3)^2/(2*E*I(1,3));0 0 1 l(1,3);0 0 0 1];
T5=[1 0 0 0;0 1 0 0;0 -j(1,5)*w^2 1 0;m(1,5)*w^2-c5*w-k5 0 0 1]*[1 l(1,4) l(1,4)^2/(2*E*I(1,4)) l(1,4)^3/(6*E*I(1,4));0 1 l(1,4)/(E*I(1,4)) l(1,4)^2/(2*E*I(1,4));0 0 1 l(1,4);0 0 0 1];
T6=[1 0 0 0;0 1 0 0;0 -j(1,6)*w^2 1 0;m(1,6)*w^2 0 0 1]*[1 l(1,5) l(1,5)^2/(2*E*I(1,5)) l(1,5)^3/(6*E*I(1,5));0 1 l(1,5)/(E*I(1,5)) l(1,5)^2/(2*E*I(1,5));0 0 1 l(1,5);0 0 0 1];
T7=[1 0 0 0;0 1 0 0;0 -j(1,7)*w^2 1 0;m(1,7)*w^2 0 0 1]*[1 l(1,6) l(1,6)^2/(2*E*I(1,6)) l(1,6)^3/(6*E*I(1,6));0 1 l(1,6)/(E*I(1,6)) l(1,6)^2/(2*E*I(1,6));0 0 1 l(1,6);0 0 0 1];
T8=[1 0 0 0;0 1 0 0;0 -j(1,8)*w^2 1 0;m(1,8)*w^2 0 0 1]*[1 l(1,7) l(1,7)^2/(2*E*I(1,7)) l(1,7)^3/(6*E*I(1,7));0 1 l(1,7)/(E*I(1,7)) l(1,7)^2/(2*E*I(1,7));0 0 1 l(1,7);0 0 0 1];
T9=[1 0 0 0;0 1 0 0;0 -j(1,9)*w^2 1 0;m(1,9)*w^2 0 0 1]*[1 l(1,8) l(1,8)^2/(2*E*I(1,8)) l(1,8)^3/(6*E*I(1,8));0 1 l(1,8)/(E*I(1,8)) l(1,8)^2/(2*E*I(1,8));0 0 1 l(1,8);0 0 0 1];
T10=[1 0 0 0;0 1 0 0;0 -j(1,10)*w^2 1 0;m(1,10)*w^2-c10*w-k10 0 0 1]*[1 l(1,9) l(1,9)^2/(2*E*I(1,9)) l(1,9)^3/(6*E*I(1,9));0 1 l(1,9)/(E*I(1,9)) l(1,9)^2/(2*E*I(1,9));0 0 1 l(1,9);0 0 0 1];
T11=[1 0 0 0;0 1 0 0;0 -j(1,11)*w^2 1 0;m(1,11)*w^2 0 0 1]*[1 l(1,10) l(1,10)^2/(2*E*I(1,10)) l(1,10)^3/(6*E*I(1,10));0 1 l(1,10)/(E*I(1,10)) l(1,10)^2/(2*E*I(1,10));0 0 1 l(1,10);0 0 0 1];
T12=[1 0 0 0;0 1 0 0;0 -j(1,12)*w^2 1 0;m(1,12)*w^2 0 0 1]*[1 l(1,11) l(1,11)^2/(2*E*I(1,11)) l(1,11)^3/(6*E*I(1,11));0 1 l(1,11)/(E*I(1,11)) l(1,11)^2/(2*E*I(1,11));0 0 1 l(1,11);0 0 0 1];
T13=[1 0 0 0;0 1 0 0;0 -j(1,13)*w^2 1 0;m(1,13)*w^2 0 0 1]*[1 l(1,12) l(1,12)^2/(2*E*I(1,12)) l(1,12)^3/(6*E*I(1,12));0 1 l(1,12)/(E*I(1,12)) l(1,12)^2/(2*E*I(1,12));0 0 1 l(1,12);0 0 0 1];
T14=[1 0 0 0;0 1 0 0;0 -j(1,14)*w^2 1 0;m(1,14)*w^2 0 0 1]*[1 l(1,13) l(1,13)^2/(2*E*I(1,13)) l(1,13)^3/(6*E*I(1,13));0 1 l(1,13)/(E*I(1,13)) l(1,13)^2/(2*E*I(1,13));0 0 1 l(1,13);0 0 0 1];
T=T14*T13*T12*T11*T10*T9*T8*T7*T6*T5*T4*T3*T2*T1;
delta=T(3,1)*T(4,2)-T(3,2)*T(4,1);
solve(delta,w);

前面的好象没问题，后面的矩阵T1，T2。。。。。T14里只有一个未知数w，然后它们进行矩阵连乘求解w，但是求不出，初学者，还请高手指点一下。
谢谢了！:victory:

hector1982 · 发表于 2007-6-3 20:08

你的主轴是什么样的啊　？是简单的盘轴结构吗？
如果是，你可以参考一下
http://forum.vibunion.com/thread-44831-1-2.html

[ 本帖最后由 hector1982 于 2007-6-3 20:14 编辑 ]

gh688 · 发表于 2007-6-9 18:01

你的程序要用一下zeros命令来定义矩阵的大小就行了.表达太烦.
另外你对求固有频率没有理解,delta=T(3,1)*T(4,2)-T(3,2)*T(4,1);
代表的是剩余量,直接这样是计算不出来的.delta=T(3,1)*T(4,2)-T(3,2)*T(4,1)=0
所对应的横坐标的点是系统的固有频率.你可以参考高等转子动力学(闻邦椿)

wy558558558 · 发表于 2007-6-11 17:37

咋看不出来是几个轮盘的呢？好象轴和盘的直径差的太小了吧．

yejet · 发表于 2007-6-14 10:09

原帖由 wy558558558 于 2007-6-11 17:37 发表
咋看不出来是几个轮盘的呢？好象轴和盘的直径差的太小了吧．

在实际中，这种转子是普遍存在的
搂主的算例中全部认为是轴段，没有轮盘附加在质量

yejet · 发表于 2007-6-14 10:10

求w的时候改二分法吧，你这样求解类似于求解解析解，有点不可思议

laolang8311 · 发表于 2010-3-23 19:45

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=[ll1 ll2 ll3 ll4 ll5 ll6 ll7 ll8 ll9 ll10 ll11 ll12 ll13 ll14 ll15 ll16] ;
Lg=[lg1 lg2 lg3 lg4 lg5 lg6 lg7 lg8 lg9 lg10 lg11 lg12 lg13 lg14 lg15 ];
Lh=[lh1 lh2 lh3 lh4 lh5 lh6 lh7 lh8 lh9 lh10 lh11 lh12 lh13];
vl=[vl1 vl2 vl3 vl4 vl5 vl6 vl7 vl8 vl9 vl10 vl11 vl12 vl13 vl14 vl15 vl16];
vg=[vg1 vg2 vg3 vg4 vg5 vg6 vg7 vg8 vg9 vg10 vg11 vg12 vg13 vg14 vg15];
vh=[vh1 vh2 vh3 vh4 vh5 vh6 vh7 vh8 vh9  vh10 vh11 vh12 vh13];
Jdl=[Jdl1 Jdl2 Jdl3 Jdl4 Jdl5 Jdl6 Jdl7 Jdl8 Jdl9 Jdl10 Jdl11 Jdl12 Jdl13 Jdl14 Jdl15 Jdl16];
Jdh=[Jdh1 Jdh2 Jdh3 Jdh4 Jdh5 Jdh6 Jdh7 Jdh8 Jdh9 Jdh10 Jdh11 Jdh12 Jdh13 ];
Jdg=[Jdg1 Jdg2 Jdg3 Jdg4 Jdg5 Jdg6 Jdg7 Jdg8 Jdg9 Jdg10 Jdg11 Jdg12 Jdg13 Jdg14 Jdg15];
Jpl=[Jpl1 Jpl2 Jpl3 Jpl4 Jpl5 Jpl6 Jpl7 Jpl8 Jpl9 Jpl10 Jpl11 Jpl12 Jpl13 Jpl14 Jpl15 Jpl16];
Jph=[Jph1 Jph2 Jph3 Jph4 Jph5 Jph6 Jph7 Jph8 Jph9 Jph10 Jph11 Jph12 Jph13 ];
Jpg=[Jpg1 Jpg2 Jpg3 Jpg4 Jpg5 Jpg6 Jpg7 Jpg8 Jpg9 Jpg10 Jpg11 Jpg12 Jpg13 Jpg14 Jpg15];
Ml=[ml1 ml2 ml3 ml4 ml5 ml6 ml7 ml8 ml9 ml10 ml11 ml12 ml13 ml14 ml15 ml16];
Mg=[mg1 mg2 mg3 mg4 mg5 mg6 mg7 mg8 mg9 mg10 mg11 mg12 mg13 mg14 mg15 ];
Mh=[mh1  mh2 mh3 mh4  mh5 mh6 mh7  mh8 mh9 mh10  mh11 mh12 mh13];
Il=[Il1 Il2 Il3 Il4 Il5 Il6 Il7 Il8 Il9 Il10 Il11 Il12 Il13 Il14 Il15 Il16];
Ig=[Ig1 Ig2 Ig3 Ig4 Ig5 Ig6 Ig7 Ig8 Ig9 Ig10 Ig11 Ig12 Ig13 Ig14 Ig15];
Ih=[Ih1 Ih2 Ih3 Ih4 Ih5 Ih6 Ih7 Ih8 Ih9 Ih10 Ih11 Ih12 Ih13];
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=[0 0 0  Kx 0 0 0 0 0 0 0 Kx 0 0 0 0];
   Kyyl=[0 0 0  Ky 0 0 0 0 0 0 0 Ky 0 0 0 0];
   Kxxg=[0 Kx 0 0 0 0 0 Kx 0 0 0 0 0 0 0];
   Kyyg=[0 Ky 0 0 0 0 0 Ky 0 0 0 0 0 0 0];
   Kxxh=[0 0 0 Kx 0 0 0 0 0 0 0 Kx 0];
   Kyyh=[0 0 0 Ky 0 0 0 0 0 0 0 Ky 0];
   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 0  Lg(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) 0  0 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

心灯 · 发表于 2010-8-6 00:37

程序代码比较恐怖啊，呵呵，采用matlab，却没有充分发挥matlab向量和矩阵运算的优势，可读性不太好，很难检查错误，也不适用于多个算例。

wrcno1 · 发表于 2010-12-16 14:18

这个可以直接建矩阵嘛，没必要每个元素单个设置

幻化成雨 · 发表于 2010-12-16 14:43

wrcno1 发表于 2010-12-16 14:18

登录/注册后可看大图

这个可以直接建矩阵嘛，没必要每个元素单个设置

个人编程习惯问题

bicy2002 · 发表于 2013-2-22 15:44

请问楼主均布质量的变截面轴是否也适用这个程序？

茹果刘儿 · 发表于 2013-3-14 10:53

求计算转子临界转速时riccati传递矩阵法的编程方法，而且，由于有耦合作用，需要计算复频率。谢谢~

mwglikai · 发表于 2013-4-2 10:00

laolang8311 发表于 2010-3-23 19:45

登录/注册后可看大图

clc
clear
%等截面轴参数

您这轴承约束考虑为主刚度系数，没有考虑交叉刚度系数，如果考虑交叉刚度系数要怎么处理？谢谢！

alisbin · 发表于 2013-4-15 21:13

非常感谢楼主

账号		自动登录	找回密码
密码			我要加入

[转子动力学] 用matlab编写的传递矩阵法求主轴的临界转速

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分

本帖被以下淘专辑推荐:

评分

回复 7楼 laolang8311 的帖子