24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

chunshui2003

木虫 (著名写手)

应助: 0 (幼儿园)
金币: 15830.9
散金: 1992
红花: 14
帖子: 1323
在线: 168小时
虫号: 1000957
注册: 2010-04-19
专业: 结构工程

[交流] 【求助】Riccati传递矩阵法求解临界转速（有部分程序）已有3人参与

我的专业是水轮发电机组振动，最近在研究Riccati传递矩阵法。但没有人交流，自己的思路总是有限，可能到了某一个节点就绕不过去了。下面的一段程序是根据闻邦椿院士编写的《高等转子动力学》P11的一个算例而来，涉及到Riccati传递矩阵法，希望在这里能有朋友帮忙看一下是哪里出了问题。毕竟我的这个专业有些偏，在算法和编程上力学的朋友更强。

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)

也许我发贴的版块可能不太对，但我觉得有一分希望永远比没有好。

[ Last edited by nono2009 on 2010-5-14 at 13:01 ]

回复此楼