一般力学和振动理论 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

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

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

ml=;
jl=;
eml=;
d=;
l=;
I=;
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=/100;
l=/100;
for i=1:13
ml(1,i)=ro*pi*d(1,i)*d(1,i)*l(1,i)/4;
end
m=;
eml=;
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=;
j=;
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=;
T2=*;
T3=*;
T4=*;
T5=*;
T6=*;
T7=*;
T8=*;
T9=*;
T10=*;
T11=*;
T12=*;
T13=*;
T14=*;
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= ;
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

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

回复 7楼 laolang8311 的帖子

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

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

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

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

wrcno1 发表于 2010-12-16 14:18 static/image/common/back.gif
这个可以直接建矩阵嘛，没必要每个元素单个设置

个人编程习惯问题

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

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

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

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

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

laolang8311 发表于 2010-3-23 19:45 static/image/common/back.gif
clc
clear
%等截面轴参数

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

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

非常感谢楼主

页: [1]

声振论坛's Archiver

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

回复 7楼 laolang8311 的帖子