振动方程无量纲之后的疑问

伤痕累累

方程如附件所示，搞不懂的就是无量纲之后的时间τ，在龙格库塔法计算时，直接用仿真设置的时间t来表示吗。还有就是，利用里面的方程和参数计算之后，结果很不一样啊。我检查不出是什么原因来。
我给程序也发上来吧，大家共同探讨一下。
function ccyy
clear all
clc
y0=zeros(4,1);
t0=[0:0.01:1000];
options=odeset('rel',1e-6);
[t,y]=ode45(@mfun,t0,y0,options);
plot(y(end-15000:end-5000,1),y(end-15000:end-5000,2))
title('轴心轨迹')
function ydot=mfun(t,x)
global R
q=3;r=2;%p为刚度比β，b为间隙δ
p=4;b=1.6e-4;
%u表示偏心量，a表示阻尼比ξ,w0为固有频率，f为摩擦系数
u=8.45e-5;f=0.0025;a=0.0025;w0=25;G=9.81/w0^2;
R=0.75;%取中间的一个参数R=0.75，测试，发现结果不是那样的，而应该是二倍周期的，根据分叉图来看。
e=sqrt(x(1)^2+x(2)^2);
if e>=b   
    fx=((e-b)^(q/r))*(x(1)-f*x(2))/e;
    fy=((e-b)^(q/r))*(f*x(1)+x(2))/e;
else fx=0;fy=0;
end
ydot=[x(3);
      x(4);
     -2*a*x(3)-x(1)+p*fx+u*R^2*cos(R*t);
     -2*a*x(4)-x(2)+p*fy+u*R^2*sin(R*t)-G;
    ];

伤痕累累

求大虾们帮忙分析一下吧。。。。

伤痕累累

求指导。。。。。。。。。

伤痕累累

顶起。。。偌大一个论坛难道没有人能回答我的问题吗？

无水1324

你的意思是，用原文中的参数不能得到论文中的类似的图？

伤痕累累

回复 5 # 无水1324 的帖子

是的。感觉他的参数是错误的。也有可能是无量纲之后的余弦、正弦项那个无量纲时间没处理好，迷惑中

无水1324

回复 6 # 伤痕累累 的帖子

那你把它全部推导一遍，如果存在问题，你就可以跟元作者联系。

伤痕累累

回复 7 # 无水1324 的帖子

我和他得疑问是一样的。不知道该怎么处理了。
请大虾帮忙
http://forum.vibunion.com/forum. ... id=94302&page=1

paopaotaiqiang

无量纲化的时间，就是转子绕着轴承转动的角度啦，计算的时候用角度就解决啦，看过你的连接了，他上面的最后求解时间长单位是时间

伤痕累累

回复 9 # paopaotaiqiang 的帖子

貌似你说的可不对哦。。。亲

paopaotaiqiang

我看了你的程序，是转自碰磨的吧，调试

n=8000

出来了吗？我也在编着方面的程序，可以交流下,这是我计算的一个短轴承，轴心，可就是稳定后的稳定圆有点小，看着不好，而且轴心和转的方向有点不符

伤痕累累

回复 11 # paopaotaiqiang 的帖子

如果有油膜力，很显然这个图是不符的、

paopaotaiqiang

我想问下，那个地方不符合呀，我刚研究，有些地方还不太清楚

伤痕累累

回复 13 # paopaotaiqiang 的帖子

你看你的轨迹图是个单位圆。说明，在非线性油膜力作用下，发生的是线型振动，显然跟实际是不符合的。不过，有可能是油膜力作用太弱吧，显示不出来。我只知道，在非线性力作用下，体现的应该是非线性的特征

paopaotaiqiang

恩，对，产生的应该是非线性特征，这是3000转的时候轴心放大的样子，中心的那个那个椭圆就是稳定时候的轨迹。我想问下，你的程序调试好了吗？让我看看你的轴心轨迹把