请教各位大侠转子-滑动轴承的碰摩问题

yangxiaokang

分岔图还是没做好，轨迹图和映射图也没做出来，到底怎么回事啊:'( ，可能还是前面方程求解有问题，还请各位多多指点

无水1324

恩，是不是那个积分进行不下去，你要仔细检查一下程序，或者调整一下参数。

yangxiaokang

是的，在做轨迹图和时程图是就发现警告说积分最小误差不能满足，程序应该没问题，估计是里面的参数有问题，参数不知道怎么调好，从去年折腾到今年

hsfy919

w=400:5:600;
for i=1:length(w)
disp(w(i));
T=2*pi/w(i);            %既然你是非自治系统，里面是sin（t）、cos（t），所以T=2*pi
tspan=[0:T/100:100*T];
x0=[0.02;0.1;0.01;0.1;0.02;0.1;0.01;0.1];  %初值
[t,y]=ode45(@fangcheng,tspan,x0,[ ],w(i));
%y0=x(end,:);  %这里不明白                   %应该是x0=y（end,:），是指每计算一个周期就把该周期的最后一次迭代数据作为初值带入下一个周期，你再试试
plot(w(i),y(4000:100:end,1),'markersize',2);%这里应该是频率随位移变化的图形           %我计算了一下，用龙格库塔法解得y值只有几十个，这一步是画不出来的，估计你函数有问题，你再检查一下
hold on;
xlabel('频率');
ylabel('位移');
end

yangxiaokang

谢谢fsfy919的指点，不过问题还是多多，图形出不来，以下是出现的问题，还请多多赐教
   400

Warning: Failure at t=1.760105e+001.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.684342e-014) at time t.
> In ode45 at 371
  In myxiugai0316 at 14
??? Error using ==> plot
Vectors must be the same lengths.

w=400:5:600;
for i=1:length(w)

disp(w(i));
T=2*pi;
tspan=[0:T/100:100*T];
x0=[0.02;0.1;0.01;0.1;0.02;0.1;0.01;0.1];  
[t,y]=ode45(@fangcheng,tspan,x0,[ ],w(i));
x0=y(end,:);  
plot(w(i),y(4000:100:end,1),'markersize',2);
hold on;
xlabel('频率');
ylabel('位移');
end

hsfy919

帮你调试了一下，程序通了，但不知道是不是你要的
%%%%%%%%%%%%%%%%%%%函数文件
function dx=fangcheng_fun(t,x)
%参数
%global w
w=760;
kc=3.5e6;
m1=4;%轴承处集中质量
m2=32.1;%圆盘处集中质量
R=0.025;%轴承半径
L=0.012;%轴承长度
c=0.00011;%油膜厚度
mu=0.018;%油膜黏度
f=0.1;%摩擦系数
c1=1050;%轴承阻尼
c2=2100;%圆盘阻尼
k=2.5e7;%弹性轴刚度
r=0.00005;%圆盘质量偏心距
delta=0.000016;%转子与定子间隙
g=9.8;
d=r/c;
G1=g/(c*w^2);
G2=g/(c*w^2);
e=sqrt(x(5)^2+x(7)^2); %圆盘中心的径向位移
%判断是否发生碰摩
H=1/2*(sign(abs(e-delta))+sign(e-delta)); %判断是否发生碰摩，sign正数返回1，负数返回-1
%碰摩力
Px=-c*(e-delta)*kc/e*(x(5)-f*x(7))*H;
Py=-c*(e-delta)*kc/e*(f*x(5)+x(7))*H;
%非线性油膜力
s=mu*w*R*L*(R/c)^2*(L/(2*R))^2;  %短轴承油膜力的sommerfeld数
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign(x(3)+2*x(2));
G=2/sqrt(1-x(1)^2-x(3)^2)*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/sqrt(1-x(1)^2-x(3)^2)));
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cos(a)+x(3)*sin(a))^2);
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G)/(1-x(1)^2-x(3)^2);
fx=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(1)*V-sin(a)*G-2*cos(a)*S);
fy=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(3)*V+cos(a)*G-2*sin(a)*S);
%方程组
dx=[x(2);
    -c1/(w*m1)*x(2)-k/(w^2*m1)*(x(1)-x(5))+s/(c*w^2*m1)*fx;
    x(4);
    -c1/(w*m1)*x(4)-k/(w^2*m1)*(x(3)-x(7))+s/(c*w^2*m1)*fy-G1;
    x(6);
    -c2/(w*m2)*x(6)-2*k/(w^2*m2)*(x(5)-x(1))+Px/(c*w^2*m2)+d*cos(t);
    x(8);
    -c2/(w*m2)*x(8)-2*k/(w^2*m2)*(x(7)-x(3))+Py/(c*w^2*m2)+d*sin(t)-G2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%下面是主程序
clear
clc
hold on;
global w f
f=400:5:600;
x0=[0.02;0.1;0.01;0.1;0.02;0.1;0.01;0.1];
for i=1:length(f)
    disp(f(i));
    w=f(i);
T=2*pi;
tspan=0:T/100:100*T;
options = odeset('RelTol', 1e-1,'AbsTol',1e-1);
[t,x]=ode45('fangcheng_fun',tspan,x0,options,w);
x0=x(end,:);  
plot(f(i),x(8000:100:end,1),'markersize',2);
xlabel('频率');
ylabel('位移');
end
%%%%%%%%%%顺便问一下，你是搞转子那部分的

hsfy919

哦，对不起，我又试了一下，主程序中options项还可以设置更精确些
'RelTol', 1e-2,'AbsTol',1e-6
当然结果越精确越好了，呵呵
另外plot(f(i),x(8000:100:end,1),'markersize',2);
改为plot(f(i),x(4000:100:end,1),'markersize',2);    %%%%%%尽量点取多点

无水1324

能否说明一下，程序那里做了修改哈哈！

yangxiaokang

非常感谢hsfy919的热心帮助，小弟在此拜谢了。我是做转子碰摩的，我们这边就我一个在做，没人商量探讨，一个人搞了很长时间也没结果，几乎灰心了。

yangxiaokang

另外还要感谢无水大哥一如既往的关注和帮助，指点了很多，谢谢

yangxiaokang

我照着hsfy919修改的程序，运算了几遍，警告说x虚部被忽略，得到的图形也有问题。分叉图计算有些难，可不可以先做时程图和轴心轨迹图，还请各位帮忙:handshake

hsfy919

那些警告不用理会，意思是在计算中得到了共轭复根，但在画图中用实数轴不能表示，故忽略，在实数域不影响绘图，况且复数域并没有研究价值。
%%%%%%%%%%%%%%%%%这是时间响应曲线
clear
clc
tspan=[0,1500];         %%%%%%%%%%%%%时间范围可以自己改变
x0=[0.2;0.1;0.01;0.1;0.02;0.1;0.01;0.1];      
options = odeset('RelTol', 1e-2,'AbsTol',1e-6);
[t,x]=ode45('fangcheng_fun',tspan,x0,options);
plot(t(10000:end),x(10000:end,1))           %%%%%%%%%%%%%%我没见你的图，也不清楚你的x1表示什么，所以先取了x1，你可以更改
%%%%%%%%%%%%%%%%%%%%%%%%%函数文件见21楼，不变
我做的图如下

yangxiaokang

按你说的可以计算时程图，跟你计算出来图形一样，不过时间范围变大如tspan=[0,8000];  图形就有问题了
%%%%%%%%%%%%%%%%%这是时间响应曲线
clear
clc
tspan=[0,1500];         %%%%%%%%%%%%%时间范围可以自己改变
x0=[0.2;0.1;0.01;0.1;0.02;0.1;0.01;0.1];      
options = odeset('RelTol', 1e-2,'AbsTol',1e-6);
[t,x]=ode45('fangcheng_fun',tspan,x0,options);
plot(t(10000:end),x(10000:end,1))           %%%%%%%%%%%%%%我没见你的图，也不清楚你的x1表示什么，所以先取了x1，你可以更改。  我这里的x1表示左轴颈中心横坐标的位移，x3是纵坐标位移；x5是圆盘中心横坐标的位移，x7是纵坐标位移。可以看看我前面的模型。


%%%%%%%%%%%%以下是我的计算程序，方程程序没变
clear
clc
tspan=[0,8000];         
x0=[0.2;0.1;0.01;0.1;0.02;0.1;0.01;0.1];      
options = odeset('RelTol', 1e-2,'AbsTol',1e-6);
[t,x]=ode45('fangcheng_fun',tspan,x0,options);

figure(1)%左轴颈x1的时程图
plot(t(4000:end),x(4000:end,1));title('时程图'); xlabel('t');ylabel('x1');  
figure(2)%左轴颈中心轨迹图
plot(x(4000:end,1),x(4000:end,3)); title('轨迹图');xlabel('x1');ylabel('y1');
figure(3)%圆盘中心轨迹图
plot(x(4000:end,5),x(4000:end,7)); title('轨迹图');xlabel('x2');ylabel('y2');

hsfy919

呵呵，我有检查了一下程序，觉得计算程序没有问题，因为同样的程序对上田振子的庞加莱截面都是正确的，请仔细检查模型，以及方程的输入。
当然还有一种情况，就是方程本身有某些特性，比如不稳定性，病态等，那样的话就不能用matlab自带的解法了，要自己编程

zdeming

不知道有没有其他收藏的方法
我只好回帖了