|
楼主 |
发表于 2006-10-26 09:44
|
显示全部楼层
%主程序:
clear all
global m1 m2 K1 K2 C1 C2 w Me F0 T M N Tn ms
m1=1640;
m2=3100;
ms=0.1*m2;
K1=2.08e6;
K2=5.6e7;
w1=(((m2*K1+(m1+ms)*K2+(m1+ms)*K1)-((m2*K1+(m1+ms)*K2+(m1+ms)*K1)^2-4*(m1+ms)*m2*K1*K2)^0.5)/(2*(m1+ms)*m2))^0.5;
f1=w1/(2*pi)
w2=(((m2*K1+(m1+ms)*K2+(m1+ms)*K1)+((m2*K1+(m1+ms)*K2+(m1+ms)*K1)^2-4*(m1+ms)*m2*K1*K2)^0.5)/(2*(m1+ms)*m2))^0.5;
f2=w2/(2*pi)
%静平衡位移
x2balance=(m1+m2)*9.8/K2;
x1balance=m1*9.8/K1+x2balance;
C1=1625;
C2=3.66e4;
w=2000*2*pi/60;
f=w/(2*pi);
T=1/f;
Me=2.1;
F0=Me*w^2;
N=2000;%积分时间上限
M=20;%T/M步长
Tn=1;%稳态的起始周期数
t0=[0:T/M:N*T];
y0=[0,0,0,0];
y=ode4('odex',[0:T/M:N*T],y0);
%[px,py]=poincare([Tn*T:T/M:N*T],y(Tn*M:N*M,:));
figure(1)
plot([Tn*T:T/M:N*T],y(Tn*M:N*M,3),'k')
xlabel('t/s')
ylabel('振动轮位移/m')
set(gca,'box','off')
%odex函数
function ydot= odex(t,y)
global m1 m2 K1 K2 C1 C2 w Me F0 T M N Tn ms
Qc=-(m1+m2)*9.8;
Xc=Qc/K2;
ydot=zeros(4,1);
ydot(1)=y(2);
ydot(2)=(C1*y(4)+K1*y(3)-C1*y(2)-K1*y(1))/m1;
ydot(3)=y(4);
if y(3)<=Xc
ydot(4)=(F0*sin(w*t)+K1*(y(1)-y(3))+C1*(y(2)-y(4))+m2*9.8)/m2;
else
ydot(4)=(F0*sin(w*t)+K1*(y(1)-y(3))+C1*(y(2)-y(4))-(C2*y(4)+K2*y(3)))/(m2+ms);
end |
评分
-
1
查看全部评分
-
|