转子-轴承油膜力程序问题，求指导

tiandong0213

    小弟最近在弄滑动轴承轴颈中心的轴心轨迹，根据论坛上的一个帖子和一篇英文文章，但是所得到的的结果不收敛，一般情况下经过无量纲话之后轴颈中心的xy方向的位移都是【-1,1】之间的，但是为什么总是发散的呢？我把那个帖子的网址贴出来，希望大家知道的能指导一下小弟哪里有问题，谢谢了~~~
由于不能贴网址，我把具体的内容移过来，希望原作者见谅，我也是想集思广益，解决问题。
function book_oilfilm_force
clear all
clc
y0=[0.01;0.01;0.0;0.0];
ff=2*pi;
fs=256*3;
tspan=[0:ff/fs:1000*ff];
options=odeset('rel',1e-6,'abs','1e-6');
[t,y]=ode45(@oilfilm,tspan,y0,options);
plot(t,y(:,1));
title('水平方向时域图')
xlabel('时间/t')
ylabel('x')
figure
plot(y(:,1),y(:,3))
title('phase figure')
xlabel('x')
ylabel('dx')
function yy=oilfilm(t,y)
%油膜力表达式
%%%%%%%%%%%%%%%%
ppp1=(y(2)+2.0*y(3))/(y(1)-2.0*y(4));
sign1=sign(ppp1);      
ppp2=y(2)+2.0*y(3);
sign2=sign(ppp2);
alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y(2)*cos(alpha)-y(1)*sin(alpha))/sqrt(abs(1.0-abs(y(1)*y(1))-abs(y(2)*y(2)))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(y(1)*y(1))-abs(y(2)*y(2))));
fv=(2.0+(y(2)*cos(alpha)-y(1)*sin(alpha))*fg)/(1.0-abs(y(1)*y(1))-abs(y(2)*y(2)));
fs=(y(1)*cos(alpha)+y(2)*sin(alpha))/(1.0-abs((y(1)*cos(alpha)+y(2)*sin(alpha))*(y(1)*cos(alpha)+y(2)*sin(alpha))));
f1=sqrt(abs(abs((y(1)-2.0*y(4))*(y(1)-2.0*y(4)))+abs((y(2)+2.0*y(3))*(y(2)+2.0*y(3)))))/(1.0-abs(y(1)*y(1))-abs(y(2)*y(2)));
fx=-1.0*f1*(3.0*y(1)*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y(2)*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);
%%%%%%%%%%%%%%%%%%%
%参数
segma=0.1;%文中的 σ
p=0.3;%小写的平，为 偏心率ρ
M=70.0;%文中的m
P=170/2.0;%转子重力，大写的P
Fx=fx/(segma*P);%无量纲化的油膜力
Fy=fy/(segma*P);
G=1/(segma*M); %无量纲化后的重力      
yy=[y(3);
    y(4);
    p*sin(t)+Fx/M;
    p*cos(t)-G+Fy/M;
    ];

运动方程

论文中结果

编程结果1

file:///d:/users/administrator/appdata/roaming/360se6/User Data/Temp/forum.php?mod=attachment&aid=NTUxNzR8NGNiNzczNzl8MTM3MDU2ODk5MnwxODM3NDR8MTExMDQ5&noupdate=yes
file:///d:/users/administrator/appdata/roaming/360se6/User Data/Temp/forum.php?mod=attachment&aid=NTUxNzR8NGNiNzczNzl8MTM3MDU2ODk5MnwxODM3NDR8MTExMDQ5&noupdate=yes

tiandong0213

自己顶一个，各位大神路过的求指点一下，小弟初搞这一块，有点迷茫，还请大家不吝赐教~~~

tiandong0213

说明一下，上面图片中前两个是论文中的运动方程和仿真结果，下面两个图是上面那个程序运行之后得到的结果

萨布江胡

这个和国内张文的《一种非稳态油膜力模型下刚性转子的分岔与混沌特性》这篇文章是一样的吧

mwglikai

请问这个程序调出来没，我改变时间间隔调后才是没有封闭

墨香

学长你最终是怎么决绝的这个发散的问题

墨香

我在做碰磨转子，闻邦春书上的算例，就是发散的

lmn789

你好！我最近也在弄滑动轴承轴颈中心的轴心轨迹程序，可否交流一下

NASA

转子轴承系统计算过程中出现发散问题是因为轴承模型中存在这根号，根号中的量在迭代过程中有可能出现负数，从而导致计算失败

相关的问题论坛前就有人讨论，而且还给出了一些处理的经验办法，大家可以找找老帖，应该会有所启发

NASA

比如：
http://forum.vibunion.com/thread-30740-1-1.html
http://forum.vibunion.com/thread-118612-1-1.html

turkeyhazel

楼主这个问题你最后解决了么 ，我现在也在做这个

怪咖先生

重复贴？？？