求助：有关求四自由度的刚性转子固有频率和振型的程序

lnly

我现在要求一个４自由度的刚性转子的固有频率，需要用ＭＡＴＬＡＢ编程．我自己编了一个但是始终有问题，总是不对，算的固有频率总是零．现在我不知道是我的程序有问题还是运动方程有问题，急着写论文，希望大家帮我看看程序是不是有问题，谢谢啊！



clc
clear
%syms w m l a b j1 j2 k1 k2 k3 k4 k5 k6
w=36000;m=4.94;l=198.1;a=98;b=100.1;j1=1.2e+3;j2=600;k1=3.15;k2=0.018;k3=243.24;k4=204;k5=12.2;k6=10.2;
M=[m*b/l,m*a/l,0,0;0,0,m*b/l,m*a/l;-j2/l,j2/l,0,0;0,0,-j2/l,j2/l];
Ｃ=[-k3*k2,-k4*k2,0,0;0,0,-k3*k2,-k4*k2;-a*k3*k2,b*k4*k2,-w*j1/l,w*j1/l;-w*j1/l,w*j1/l,-a*k3*k2,b*k4*k2];
K=[-(k5+k3*k1),-(k6+k4*k1),0,0;0,0,-(k5+k3*k1),-(k6+k4*k1);-a*(k3*k1+k5),b*(k6+k4*k1),0,0;0,0,-a*(k3*k1+k5),b*(k6+k4*k1)];
I=eye(4);
mm=zeros(4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%转换为标准特征值问题
Ｇ=[mm,M;-inv(M)*K,-inv(M)*Ｃ];
cri=eig(Ｇ);
pp=imag((cri));
wn=(abs(imag(cri)));
firstwn=wn(1)     %%%第一阶固有频率
kkk=min(abs((real(cri))))
secondwn=wn(3)   %%%第二阶固有频率
kkk=min(abs((real(cri))));
thirdwn=wn(5)    %%%第三阶固有频率

[ 本帖最后由 xinyuxf 于 2007-5-30 11:44 编辑 ]

yejet

1. clc
   
2. clear
   
3. %syms w m l a b j1 j2 k1 k2 k3 k4 k5 k6
   
4. w=36000;m=4.94;l=198.1;a=98;b=100.1;j1=1.2e+3;j2=600;k1=3.15;k2=0.018;k3=243.24;k4=204;k5=12.2;k6=10.2;
   
5. M=[m*b/l,m*a/l,0,0;0,0,m*b/l,m*a/l;-j2/l,j2/l,0,0;0,0,-j2/l,j2/l];
   
6. C=[-k3*k2,-k4*k2,0,0;0,0,-k3*k2,-k4*k2;-a*k3*k2,b*k4*k2,-w*j1/l,w*j1/l;-w*j1/l,w*j1/l,-a*k3*k2,b*k4*k2];
   
7. K=[-(k5+k3*k1),-(k6+k4*k1),0,0;0,0,-(k5+k3*k1),-(k6+k4*k1);-a*(k3*k1+k5),b*(k6+k4*k1),0,0;0,0,-a*(k3*k1+k5),b*(k6+k4*k1)];
   
8. I=eye(4);
   
9. mm=zeros(4);
   
10. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%转换为标准特征值问题
    
11. G=[mm,M;-inv(M)*K,-inv(M)*C];
    
12. cri=eig(G);
    
13. 
    
14. wn=cri(imag(cri)==0);
    
15. wn=wn (wn > 0);
    
16. wn=sort(wn)

复制代码

lnly

谢谢你的指教！但我还有个问题想请教！
为什么固有频率取得是特征值为正实数的部分而不是复数特征值的虚数部分啊？
不是说如果刚度，阻尼矩阵（例如比例阻尼）是对称矩阵，算出来的特征值是实数，则特征值是系统的固有频率，单位是rad/s.如果系统的刚度阻尼矩阵不是对称矩阵，则特征值是复数。例如考虑了转子系统考虑了陀螺效应，引入了负对称矩阵，考虑了不对称的轴承油膜力。此时特征值实部表示系统的衰减特性，正负可以表征系统的稳定与否，虚部表示系统固有频率。如果原振动方程转化到状态空间中，方程维数扩为原2倍，特征值为一对对共轭复数，具有负虚部的不考虑，选择具有正虚部的来计算系统的稳定性和固有频率。
而我这个怎么算出来的特征值既有实数又有复数啊？