Riccati传递矩阵

奔腾6

     各位师兄师姐：大家好，很抱歉打扰你们，我最近开始做课题，是转子动力学方面的，用matlab编程求临界转速和振型，是按照闻邦椿《高等转子动力学》上面的Riccati传递矩阵公式编的，算例也是用上面的数据，但是得出的剩余量曲线和书上的相差很大，反复检查好长时间，就是不知哪里的问题，很头疼啊，只能向前辈请教了，请帮忙看一看，或者有正确的程序给我发一个作为参考，多谢，祝身体健康。
clear;
l=[1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,0];      %13个结点
m=[2940,5880,5880,5880,5880,5880,5880,5880,5880,5880,5880,5880,2940];
Jp=[0,0,0,0,0,0,0,0,0,0,0,0,0];
Jd=[0,0,0,0,0,0,0,0,0,0,0,0,0];       %不计转动惯量和陀螺力矩
I=[100,100,100,100,100,100,100,100,100,100,100,100,100];
E=4393;
v=[0,0,0,0,0,0,0,0,0,0,0,0,0];     %不计剪切影响
k=[1.96*10^9,0,0,1.96*10^9,0,0,1.96*10^9,0,0,1.96*10^9,0,0,1.96*10^9];    %油膜刚度
kb=[2.7048*10^9,0,0,2.7048*10^9,0,0,2.7048*10^9,0,0,2.7048*10^9,0,0,2.7048*10^9];       %轴承座刚度
mb=[3577,3577,3577,3577,3577,3577,3577,3577,3577,3577,3577,3577,3577];     %参振质量
S=[0,0;0,0];    %Riccati第一矩阵
s=1;
x=[];     %用于记录剩余量
for n=1864:1864       %试算频率
for i=1:13
K=k(i)*(kb(i)-mb(i)*n^2)/(k(i)+kb(i)-mb(i)*n^2)      %总刚度
u11=[1,l(i);0,1];
u12=[l(i)*(m(i)*n^2-K),(Jp(i)-Jd(i))*n^2;m(i)*n^2-K,0];
u21=(l(i)/(E*I(i))).*[l(i)/2,l(i)^2*(1-v(i))/6;1,l(i)/2];
u22=[1+l(i)^3*(1-v(i))*(m(i)*n^2-K)/(6*E*I(i)),l(i)+l(i)^2*(Jp(i)-Jd(i))*n^2/(2*E*I(i));l(i)^2*(m(i)*n^2-K)/(2*E*I(i)),1+l(i)*(Jp(i)-Jd(i))*n^2/(E*I(i))];
S=[u11*S+u12]*inv([u21*S+u22])
end
D=det(S);
x=[x,D];        %记录剩余量
end
n=1864:1:1864;     %产生曲线的横坐标
grid on
plot(n,x)

sssssxxxxx921

什么年代了怎么还用这种方法做啊  找九十年代的研究生问吧

奔腾6

那用什么方法？？？

hector1982

看＜转子动力学＞吧，
用传递矩阵法吧
参考
http://forum.vibunion.com/forum/ ... 6amp%3Btypeid%3D157

yangyongwhu@163

Riccati传递矩阵
比普通的传递矩阵法更为优化啊，你们不知道的吗？Riccati传递矩阵
法就是在普通传递矩阵法的基础上修改而来的

lhy

Riccati传递矩阵法及其改进算法在进行振动分析时，只需对一些阶次很低的传递矩阵进行连续的矩阵乘法运算，并且和系统自由度没有直接关系，大大节省计算工作量，特别是在只要求计算若干阶低阶的固有频率和振型时更加方便。Riccati传递矩阵法的基本思想是：把系统分割成一系列具有简单动力学特性的单元（两端，三端或者多端），振动时系统的状态可以用该单元端点的状态矢量来表示，相邻单元之间状态矢量的关系，用传递矩阵表示。

yejet

sssssxxxxx921 发表于 2007-7-26 19:34

                               
登录/注册后可看大图

什么年代了怎么还用这种方法做啊  找九十年代的研究生问吧

看来你已经很高级了，不知你用什么方法呢

yejet

一个例子，仅作参考，比保证程序的正确性

1. clc
   
2. clear
   
3. format short e
   
4. m1=2940;m13=2940;m2=5880;
   
5. M=[m1 m2 m2 m2 m2 m2 m2 m2 m2 m2 m2 m2 m13];
   
6. l=1.3;
   
7. L=[l l l l l l l l l l l l 0];
   
8. a=2.9592e-9;
   
9. A=[a a a a a a a a a a a a a];
   
10. k1=1.960000000e9;
    
11. k2=2.7048e9;
    
12. mb=3577;
    
13. k=0;
    
14. S=[0 0;0 0];
    
15. for w=0:0.01:5000
    
16.     for i=1:13
    
17.        k3=1.0e9;
    
18.        K=[k3 0 0 k3 0 0 k3 0 0 k3 0 0 k3];  
    
19.     U11=[1 L(i);0 1];
    
20.     U12=[L(i)*(M(i)*w^2-K(i)) 0;(M(i)*w^2-K(i)) 0];
    
21.     U21=[L(i)*A(i)/2 L(i)^2*A(i)/6;A(i) L(i)*A(i)/2];
    
22.     U22=[1+L(i)^2*A(i)*(M(i)*w^2-K(i))/6 L(i);L(i)*A(i)*(M(i)*w^2-K(i))/2 1];
    
23.     S=(U11*S+U12)/(U21*S+U22);
    
24.    end
    
25.      F=S(1,1)*S(2,2)-S(1,2)*S(2,1);
    
26. if F*(-1)^k< 0 %求解临界转速
    
27.        k=k+1;
    
28.       wi(k)=w;
    
29.        w=wi(k)
    
30.        ni(k)=wi(k)*30/pi
    
31.     end
    
32.   
    
33. end

复制代码

huangziyuan212

huangziyuan212

单元长度为0是什么意思啊？？？

璀璨星辰

huangziyuan212

dollfish000

二楼挺搞笑，我们要相互学习，取长补短！