|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
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)
这是本论坛某位朋友根据闻梆春院士《高等转子动力学》中P11算例3.2.3编写的程序，当时他有几个问题没有解决。今天再提出来是因为本人最近一直在看这方面的内容，理论已经看了好几次，但总有几点搞不清楚，所以希望大家能给解答一下。
1试算频率
步长取多大为宜，上面红色的 for n=1864:1864是否有问题，如果是试算的话，那么这个循环等于没有做，只是计算了频率n=1864时的情况，不知道我理解是否正确
2剩余量方程或者频率方程
这里是我最糊涂的地方。首先，剩余量是什么意思；其次，如果求解这个剩余量方程，里面除了频率n是未知以外，其余都是已知的，为什么还要进行试算，取不同的频率值n。这个剩余量方程究竟要求的是什么呢。还有，D=det(S); 这是要求解特征值吗？
3书上说“以一定的步长选取试算频率，用式（1.3-22）递推，算得剩余量|S|N+1值，就可以画出1.3-6那样剩余量S|N+1随频率变化的曲线，曲线和横坐标的各个交点就是所求的各阶临界角速度”
这里我还是不明白，算得|S|N+1的值是指式（1.3-22）中的Si+1吗，否则按照式（1.3-24），|S|N+1是等于0啊？
4上面的程序在D=det(S) 之后不清楚是什么意思，是求交点吗？请解释一下。
另外，这个程序的作者已说过程序存在问题，我在matlab中运行后得到结果如“Riccati结果.jpg”所示。
我希望论坛里有朋友能够帮忙解答一下，我不是一个懒惰或者不劳而获的人，我是真心想学东西，但苦于自己的理解能力有限，所以渴望有人能够给指引一下，也许大家的一句话就能够让我少走很多弯路。当我掌握后也会无私地将经验同大家分享的，希望版主和各位朋友帮助，谢谢！
|
-
Riccati递推公式
-
系统频率方程
-
剩余量曲线
-
根据上面程序得到的结果
|
发表于 2010-5-14 14:30
提升卡
变色卡
显身卡
楼主