本帖最后由 肥振 于 2011-8-11 11:24 编辑
请问各位牛人，小弟把楼主的程序（红色部分）改了下，为什么出不了相似结果呢？是不是我在哪些地方出错了？
代码如下：
%主分析程序
clear all
h=pi/256;
w=2100;
tf=300000*2*pi/w;
tspan=0:h:tf;
y0=[0.05,0.5,0.5,0.05,0.05,0.5,0.5,0.05,0.05,0.5,0.5,0.05,0.05,0.5];
%options=odeset('RelTol',10^-6,'AbsTol',10^-6);
[t,y]=ode45('fun11',tspan,y0);
figure;
subplot(2,2,1);
plot(y(:,1),y(:,3))
title('未松动端轴心轨迹');xlabel('x1');ylabel('y1');
subplot(2,2,2);
plot(y(:,9),y(:,11))
title('松动端轴心轨迹');xlabel('x3');ylabel('y3');
subplot(2,2,3);
plot(y(:,5),y(:,7))
title('圆盘轴心轨迹');xlabel('x2');ylabel('y2');
fun11子函数：
function d=fun11(t,y)
d=zeros(14,1);
%N=length(y);
w=2100;
m1=4;%两端滑动轴承处等效集中质量
m2=32.1; %转子圆盘等效集中质量
m3=50.0;%轴承支座处等效集中质量
g=9.81;
e=0.00005; %偏心距
k=2.5e7;%弹性轴刚度
delta2=0.6e-3;%初始间隙
c=delta2;
c1=1050;%转子圆盘处阻尼系数
c2=2100;%转子在轴承处阻尼系数
k1=7.5e7;
k2=2.5e9;
cb1=350;
cb2=500;
ox1=y(1);%未松动端竖直方向位移x1
odx1=y(2);
oy1=y(3);%未松动端竖直方向位移y1
ody1=y(4);
ox2=y(5);%圆盘位移x2
odx2=y(6);
oy2=y(7);%圆盘位移y2
ody2=y(8);
ox3=y(9);%松动端轴心位移x3
odx3=y(10);
oy3=y(11);%松动端轴心位移y3
ody3=y(12);
oy4=y(13);%质量m3在竖直方向位移y4
ody4=y(14);
if oy4<0
cb=cb2;
kb=k2;
elseif (oy4>=0)&(oy4<=delta2)
cb=0;
kb=0;
else
cb=cb1;
kb=k1;
end
c1=c1/w;
c2=c2/w;
cb=cb/w;
k=k/w^2;
kb=kb/w^2;
m1=m1/c/w^2;
m2=m2/c/w^2;
m3=m3/c/w^2;
fx=oilx( ox1, oy1, odx1, ody1, w);
fy=oily( ox1, oy1, odx1, ody1, w);
fx1=oilx( ox3,oy3-oy4,odx3,ody3-ody4,w);
fy1=oily( ox3,oy3-oy4,odx3,ody3-ody4,w);
fx=fx/c/w^2;
fy=fy/c/w^2;
fx1=fx1/c/w^2;
fy1=fy1/c/w^2;
d(1)=odx1; %主要是这里改动
d(2)=-(c1/m1)*odx1-(k/m1)*(ox1-ox2)+(fx/m1);
d(3)=ody1;
d(4)=-(c1/m1)*ody1-(k/m1)*(oy1-oy2)+(fy/m1)-g;
d(5)=odx2;
d(6)=-(c2/m2)*odx2-(k/m2)*(ox2-ox1)-(k/m2)*(ox2-ox3)+e*w^2*cos(w*t);
d(7)=ody2;
d(8)=-(c2/m2)*ody2-(k/m2)*(oy2-oy1)-(k/m2)*(oy2-oy3)+e*w^2*sin(w*t)-g;
d(9)=odx3;
d(10)=-(c1/m1)*odx3-(k/m1)*(ox3-ox2)+(fx1/m1);
d(11)=ody3;
d(12)=-(c1/m1)*ody3-(k/m1)*(oy3-oy2)+(fy1/m1)-g;
d(13)=ody4;
d(14)=-(cb/m3)*ody4-(kb/m3)*oy4-(fy1/m3)-g;
t % 显示时间
end
油膜力子函数oilx，oily不变的
谢谢各位了！！！ |
[复制链接]
发表于 2010-10-28 21:30
显身卡
发表于 2010-10-29 11:33
发表于 2010-11-13 10:19
好东西啊！