马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
这是我做的时间历程,相图和幅值谱,可是前面两个和频谱对不上,从相图可以看到theta的运动是4倍周期运动,可频谱出现了5种频率成分,这应该怎么解释呀,高手帮忙给解决一下。
下面是我的程序
function dy=moxing3(t,y,flag)
global c;
theta0=pi*0.1;
g=10;
p=1.29;
T=26000;
ml=1.7;h=10/1000;
a=0.4;
dd=0.028;
mi=dd*h*0.9;
D=dd+h*sin(theta0);
L=244;
U=17;
ky=4*T*3/L;
kt=864;
kx=4*T*3/L;
S=3*mi*dd/2*L*(2/pi)^2;
m=(2/pi)^2*L*3*(ml+mi);
J=(2/pi)^2*L*(ml+mi)*a^2;
xiy=0.04;
xix=0.04;
xit=0.06;
omegay=sqrt(ky/m);
omegax=sqrt(kx/m);
omegat=sqrt(kt/J);
a1=-4.5;
a2=1;
a3=2;
b1=-1;
b2=-1;
b3=1;
c1=-5.2;
c2=3;
c3=1;
theta1=y(5)+theta0-y(2)/U-(dd/2/U)*y(6);
cy=a1*theta1+a2*theta1^2+a3*theta1^3;
cx=b1*theta1+b2*theta1^2+b3*theta1^3;
cm=c1*theta1+c2*theta1^2+c3*theta1^3;
Fy=1/2*p*U^2*D*cy*L;
Fx=1/2*p*U^2*D*cx*L;
M=1/2*p*U^2*D^2*cm*L;
func1=y(2);
func2=-(-(2*m*xiy*omegay*y(2)+ky*y(1)-Fy)*S^2*sin(theta0)^2+S^2*cos(theta0)*sin(theta0)*(2*m*xix*omegax*y(4)+kx*y(3)-Fx)-S*cos(theta0)*(2*J*xit*omegat*y(6)+kt*y(5)-M)*m+(2*m*xiy*omegay*y(2)+ky*y(1)-Fy)*J*m)/(-S^2+J*m)/m;
func3=y(4);
func4=(-S^2*sin(theta0)*cos(theta0)*(2*m*xiy*omegay*y(2)+ky*y(1)-Fy)-S^2*sin(theta0)^2*(2*m*xix*omegax*y(4)+kx*y(3)-Fx)+S*sin(theta0)*(2*J*xit*omegat*y(6)+kt*y(5)-M)*m+(2*m*xix*omegax*y(4)+kx*y(3)-Fx)*S^2-(2*m*xix*omegax*y(4)+kx*y(3)-Fx)*J*m)/(-S^2+J*m)/m;
func5=y(6);
func6=-(-S*cos(theta0)*(2*m*xiy*omegay*y(2)+ky*y(1)-Fy)-S*sin(theta0)*(2*m*xix*omegax*y(4)+kx*y(3)-Fx)+(2*J*xit*omegat*y(6)+kt*y(5)-M)*m)/(-S^2+J*m);
dy=[func1;func2;func3;func4;func5;func6];
clear;
step=0.01;
tstop=1000;
y0=rand(1,6);
pt=500*1/step;
[t,y]=ode45('moxing3',[0:step:tstop] ,y0,[]);
set(gcf,'color','w')
subplot(2,1,1);
plot(t,y(:,5));
xlabel('t')
ylabel('theta')
axis([900 1000 -Inf Inf])
grid;
subplot(2,2,3);
plot(y(end-pt:end,5),y(end-pt:end,6));
xlabel('theta')
ylabel('dtheta')
grid;
subplot(2,2,4);
yy=fft(y(end-pt:end,5));
N=length(yy);
power=abs(yy);
freq=(1:N-1)*1/step/N;
plot(freq(1:N/2),power(1:N/2));
xlabel('f(theta)')
ylabel('')
axis([0 1.5 0 Inf]) |