临界转速Matlab计算程序

用matlab编写的传递矩阵法求主轴的临界转速
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];
T1
2=[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，但是求不出，初学者，还请高手指点一下。




你的程序要用一下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
所对应的横坐标的点是系统的固有频率.你可以参考 高等转子动力学(闻邦椿)