为什么激励频率和响应频率不一样

bruce001.

我算单盘对称转子的动力特性，计算出来的稳态响应的频率和激励频率无关，一直稳定在0.08Hz附近，。高手们支支招啊程序：
function a
clc
clear
T=2*pi/1000;
tt=0:T/50:100000*T;
hold on
[t,yy]=ode45('fun',tt,[0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001]);
subplot(411)   
plot(t,yy(:,1));
title('轮盘处x向时域波形')
grid on
subplot(412)
plot(t,yy(:,2));
title('轮盘处y向时域波形')
grid on
subplot(413)
plot(t,yy(:,3));
title('轴承处x向时域波形')
grid on
subplot(414)
plot(t,yy(:,4));
title('轴承处y向时域波形')
grid on


function ydot=fun(t,y)
c=y(1);d=y(2);e=y(3);f=y(4);
c1=y(5);d1=y(6);e1=y(7);f1=y(8);
% 短轴承非线性油膜力解析式
alpha=atan((f+2*e1)/(e-2*f1))-pi/2*sign((f+2*e1)/(e-2*f1))-pi/2*sign(f+2*e1);
S=(e*cos(alpha)+f*sin(alpha))/(1-(e*cos(alpha)+f*sin(alpha))^2);
G=2/(1-e^2-f^2)^0.5*(pi/2+atan((f*cos(alpha)-e*sin(alpha))/(1-e^2-f^2)^0.5));
V=(2+(f*cos(alpha)-e*sin(alpha))*G)/(1-e^2-f^2);
Z=-((e-2*f1)^2+(f+2*e1)^2)^0.5/(1-e^2-f^2);
fx=Z*(3*e*V-sin(alpha)*G-2*cos(alpha)*S);
fy=Z*(3*f*V+cos(alpha)*G-2*sin(alpha)*S);

miu=0.018;             % 滑油粘度
omega=1000;         % 转速
r=0.057;                 % 转子轴半径/圆盘内径
cc=0.0002;             % 轴承间隙
L=0.0285;               % 轴承长度。文章里还有个几何长度，意义？
C=3000;                  % 轴承处结构阻尼
m1=420;m2=50;       % 质量
rou=0.1;                  % 无量纲偏心
a=3e3/omega;
k=2.105e8/2/omega^2;
tau=omega*t;
G1=9.8/(cc*omega^2);
delta=miu*omega*r*L*(r/cc)^2*(L/(2*r))^2;
m22=m2*omega^2*cc/delta;

ydot=[c1;
    d1;
    e1;
    f1;
    -a/m1*c1-2*k/m1*(c-e)+rou*cos(tau);         % 轮盘
    -a/m1*d1-2*k/m1*(d-f)+rou*sin(tau)-G1;
    -k/m2*(e-c)+1/m22*fx;                       % 轴承处
    -k/m2*(f-d)+1/m22*fy-G1;
    ];


运动微分方程如上图

bruce001.

这是运动微分方程

bruce001.

上面被论坛签名挡住了

玉林

如果是持续激励的话，那么程序中应该有一项多余的激励，并且比真实激励要大。波形图中响应并没有衰减，那么就是说跟阻尼无关。

犟牛

从结果上看最大的可能还是某个或几个参数搞错了，自己检查吧

伤痕累累

tau=omega*t;
这是无量纲的时间吗？这样表示是不对的
从你的时域图上看周期很大，方程表示的周期应该是T=2π/omega，应该很小的周期啊

bruce001.

玉林 发表于 2012-11-2 10:39

                               
登录/注册后可看大图

如果是持续激励的话，那么程序中应该有一项多余的激励，并且比真实激励要大。波形图中响应并没有衰减，那么 ...

我现在发现上面的波形里还有频率为omega/2pi的成分的，只是振幅太小，画成这样看不出来了。

bruce001.

犟牛 发表于 2012-11-2 15:59

                               
登录/注册后可看大图

从结果上看最大的可能还是某个或几个参数搞错了，自己检查吧

我严格按照一篇文章上的参数写的，真是没有头绪啊。请问哪位有类似的数据，或者程序？我想了解一下流程。

bruce001.

伤痕累累 发表于 2012-11-4 16:59

                               
登录/注册后可看大图

tau=omega*t;
这是无量纲的时间吗？这样表示是不对的
从你的时域图上看周期很大，方程表示的周期应该是T= ...

是的，我的理解是，至少以一倍频为主，就是周期为T=2pi/omega，可以很明显看出这个频率成分才对。还有就是，我的计算时间取的是t=0:T/100:1000*T。我感觉这跟无量纲时间tau=omega*t没关系吧？总感觉时间无量纲化没有必要。请问正确的tau应该是怎样的呀？

伤痕累累

bruce001. 发表于 2012-11-5 10:37

                               
登录/注册后可看大图

是的，我的理解是，至少以一倍频为主，就是周期为T=2pi/omega，可以很明显看出这个频率成分才对。还有就是 ...

tao=1*t;

玉林

做一下fft看看有多少频率、

tiandong0213

你好！能否告知一下您是参考了哪篇文章的呢？

bruce001.

tiandong0213 发表于 2013-6-5 10:27

                               
登录/注册后可看大图

你好！能否告知一下您是参考了哪篇文章的呢？

荆建平等，《连续转子轴承系统的非线性动力学行为研究》，热能与动力工程，2002.11   17（102）：607~610.  我是用的里面的数据计算的。问题找到了，是在fun子函数里时间tau没有无量纲化。

tiandong0213

bruce001. 发表于 2013-6-5 15:26

                               
登录/注册后可看大图

荆建平等，《连续转子轴承系统的非线性动力学行为研究》，热能与动力工程，2002.11   17（102）：607~610 ...

谢谢您的回答，小弟最近在做滑动轴承轴心轨迹的分析，可否将您的最终程序发我一份呢？我的邮箱是234388662@qq.com。谢谢了~~