请大家看一下，比较哪一种分岔程序更优

chunshui2003

本来以为已经很接近了，但算着算着就开始不断怀疑，陆陆续续弄了好久也对不上相图和分岔图的结果。换了不同的计算方法，更改精度、扩大计算周期、改变初值......但仍然无法得到满意的结果。把之前liliangbiao前辈的分岔程序和频闪法的分岔程序计算同一系统，希望大家看一下哪一种更优、更合理，谢谢了。

系统运动微分方程：因为在2009a平台运算，所以采用内嵌函数

1. 
   
2. function  uu = equ_elastic_without_UMP(T,u,e0,cz,Lb)
   
3. m1 = 60;               %转子质量 kg
   
4.   m2 = 25;               %轴颈质量 kg
   
5.   c1 = 4000;             %转子阻尼 N.s/m
   
6.   c2 = 1200;             %轴承阻尼 N.s/m
   
7.   Ke = 6.2*10^6;          %转轴刚度 N/m
   
8. %   e0 = 0.6*10^-3;        %转子质量偏心 m
   
9. %   Rr = 60*10^-3;         %转子半径 m
   
10. %   Lr = 150*10^-3;        %转子长 m
    
11.   delta_0 = 4.5*10^-3;   %均匀气隙大小 m
    
12. %   cz = 0.2*10^-3;        %轴颈间隙 m
    
13.   Rb = 50*10^-3;         %轴承半径 m
    
14. %   Lb = 20*10^-3;        %轴承长 m
    
15.   mu = 18*10^-3;         %绝对润滑油粘度 无单位
    
16.   
    
17.   omega = 13;            %大轴旋转角速度
    
18.   
    
19. %   global cz
    
20.   
    
21.   
    
22.   sigma = mu * omega * Rb * Lb * (Rb/cz)^2 * (Lb/2/Rb)^2;        %Sommerfeld修正系数
    
23.   
    
24.   A1 = u(7) + 2*u(6);
    
25.   A2 = u(5) - 2*u(8);
    
26.   
    
27.   alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1);
    
28.   
    
29. E = sqrt(u(5).^2 +u(7).^2);                     %轴颈轴心轨迹
    
30. E_minus = sqrt(1 -E.^2);
    
31. 
    
32. B1 = u(7)*cos(alpha) - u(5)*sin(alpha);
    
33. B2 = u(5)*cos(alpha) + u(7)*sin(alpha);
    
34.   
    
35.   G = 2./E_minus * ( pi/2 + atan(B1./ E_minus) );
    
36.   V = (2 + B1 * G) ./  E_minus.^2;
    
37.   S = B2 ./ ( 1 - B2.^2 );
    
38.   
    
39.   F_oil_same = -sqrt( A2.^2 + A1.^2 ) ./ E_minus.^2;
    
40.   
    
41.   f_x = F_oil_same * ( 3* u(5) *V - sin(alpha) * G - 2 *cos(alpha)*S );
    
42.   f_y = F_oil_same * ( 3* u(7) *V + cos(alpha) * G - 2 *sin(alpha)*S );     %油膜力的无量纲表达式
    
43.   
    
44.   fx = sigma * f_x;
    
45.   fy = sigma * f_y;              %非线性油膜力
    
46.   
    
47.   uu = zeros(8,1);
    
48.   uu(1) = u(2);
    
49.   uu(2) = -c1/(m1*omega) * u(2) - Ke/(m1*omega^2) * u(1) + Ke * cz/(m1*omega^2*delta_0)*u(5) + e0 * cos(T)/delta_0 ;
    
50.   uu(3) = u(4);
    
51.   uu(4) = -c1/(m1*omega) * u(4) - Ke/(m1*omega^2) * u(3) + Ke * cz/(m1*omega^2*delta_0)*u(7) + e0 * sin(T)/delta_0 ;
    
52.   uu(5) = u(6);
    
53.   uu(6) = -c2/(m2*omega) * u(6) + Ke * delta_0/(2*m2*omega^2*cz) * u(1) - Ke/(2*m2*omega^2) * u(5) + fx/(m2*omega^2 *cz);
    
54.   uu(7) = u(8);
    
55.   uu(8) = -c2/(m2*omega) * u(8) + Ke * delta_0/(2*m2*omega^2*cz) * u(3) - Ke/(2*m2*omega^2) * u(7) + fy/(m2*omega^2 *cz);
    

复制代码

频闪法：

1. 
   
2. function  [time,y] = jie_without_UMP
   
3. tic
   
4. period = 2*pi;               %%%因为对时间t采用了无量纲化，cos(oemga*t)变成了cos(T)，所以周期变为2*pi/1，暂时这样
   
5. tspan = (0:period/100:1000*period);
   
6. y0 = [0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001];
   
7. e0 = 0.6*10^-3; Lb = 100*10^-3;         
   
8.       
   
9. cz = 0.0002:0.0001:0.00600;
   
10. options = odeset('Reltol',1e-3,'Abstol',1e-6);
    
11. 
    
12. row = length(cz);
    
13. column = length(80100:100:100001);
    
14. U = zeros(row,column);
    
15. 
    
16. for i = 1:length(cz)
    
17.     disp(cz(i))
    
18. [time,y] = ode15s(@equ_elastic_without_UMP,tspan,y0,options,e0,cz(i),Lb);
    
19. y0 = y(end,:);
    
20. U(i,:) = y(80100:100:end,1)';
    
21. end
    
22. plot(cz,U,'r.','markersize',5);
    
23. saveas(gcf,'分岔图.bmp','bmp');
    
24. toc
    

复制代码

liliangbiao前辈的单个周期计算法（我是这么理解的，所以这样叫）：

1. 
   
2. function  linshi2
   
3. 
   
4. 
   
5. tic
   
6. period = 2*pi;               %%%因为对时间t采用了无量纲化，cos(oemga*t)变成了cos(T)，所以周期变为2*pi/1，暂时这样
   
7. e0 = 0.6*10^-3; Lb = 100*10^-3;         %Situ_H_imp
   
8. cz = 0.0002:0.0001:0.0060;
   
9. 
   
10. options = odeset('Reltol',1e-3,'Abstol',1e-6);
    
11. 
    
12. k = 0;
    
13. step = period/100;
    
14. 
    
15. for ww = 1:length(cz)
    
16.     y0 = [0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001];
    
17.     disp(cz(ww))
    
18.     k = k+1;
    
19.     tspan = 0:step:800*period;
    
20.     [t,y] = ode15s(@equ_elastic_without_UMP,tspan,y0,options,e0,cz(ww),Lb);
    
21.     y0 = y(end,:);
    
22.     j = 1;
    
23.     for i = 800:1000
    
24.         tspan = i*period:step:(i+1)*period;
    
25.         [t,y] = ode15s(@equ_elastic_without_UMP,tspan,y0,options,e0,cz(ww),Lb);
    
26.         YY1(k,j)=y(end,1);   
    
27.         j = j+1;
    
28.        y0 = y(end,:);
    
29.     end
    
30. 
    
31.    
    
32. end
    
33. bifdata = YY1(:,end-20:end);
    
34. plot(cz,bifdata,'r.','markersize',5);
    
35. toc
    

复制代码

频闪法

liliangbiao法

chunshui2003

以下分别是cz=0.0030时，采用liliangbiao法和频闪法得到的poincare图和轨迹图，请大家看一下。

liliangbiao法poincare

频闪法轨迹

频闪法poincare

chunshui2003

图片的上传方法搞不懂了。我采用的是图片上传，加入了描述，但有的图片就有名称，有的就没有，另外也没有按照我的顺序显示。请问该如何操作？

octopussheng

截面中，点的坐标都很接近，应该就是一个点，是周期运动

chunshui2003

回复 4 # octopussheng 的帖子

感谢你，学长，这方面的问题最近一直都是你在提供给我帮助，让我觉得上论坛是一件很值得期待的事。再向你问几个问题：
上面轨迹和映射图中，第一幅和第四幅对应着liliangbiao前辈的程序计算而来；第二幅和第三幅对应着频闪法而来。
学长，按你说的截面中点的坐标都很接近，是否指的是liliangbiao法的映射图。在频闪法的映射图中-0.0154——-0.0156之间的点距大概为0.2*10^-3量级，根据横坐标的量级为10^-3，这样的差距可以理解为点很接近吗？
另外，学长，从上面四幅图看，是不是liliangbiao前辈的方法更加好一些呢？（对于我这个微分方程计算而言）

octopussheng

具体哪个好我不敢说，但对于你这个系统，应该说分岔的图的效果还没有出来。
如果你把ode程序的计算误差设的小一点的话，计算可能会有问题。

chunshui2003

回复 6 # octopussheng 的帖子

学长你好，分岔图的控制参数为cz = 0.0002:0.0001:0.0060，大概画了200个周期的点。这是我为了方便在matlab画图，所以将控制参数的间隔取的比较大，实际上控制参数

1. 
   
2. cz = 0.2*10^-3:0.005*10^-3:6*10^-3;
   

复制代码

需要计算1161次，耗时较长。
如果分岔的效果是因为间隔太大而没有出来的话，这个可以解决；如果是周期取少了的话也可以解决。
学长的意思是指的哪一方面呢？或者是没有出现由周期运动迈向周期二运动、到周期N（或者拟周期），最终出现混沌吗？
请学长明示。

hsfy919

回复 7 # chunshui2003 的帖子

两种程序的思路应该都没错，你试着调高积分精度试试。在非线性油膜力作用下，系统的轴承轨迹为什么差距那么小，控制参数cz代表的是什么，建议你用频谱分析一下，看看频率构成

chunshui2003

回复 8 # hsfy919 的帖子

转子-轴承系统

我画的是转子轴心的轨迹，你说的差距小我不太明白是什么意思。因为只是在cz=0.003时采用两种方法得到的结果，并没有涉及到参数的改变。
cz是轴承与轴颈之间的间隙。
另外，我尝试过将abstol改为1e-7，这样精度提高了，但是计算速度也下降了，并且好像没有对轨迹产生太明显的影响。
我怀疑可能是系统其他的部分参数设置不合理，准备再改变一下。

jgwang

方程似乎没问题,主要还是参数的选择上如ke,delta,omega

chunshui2003

回复 10 # jgwang 的帖子

恩，说的有道理，之前一直没考虑这些问题，现在重新弄一下，谢谢提醒。

hsfy919

回复 9 # chunshui2003 的帖子

我说的差距小是指轴心轨迹的坐标变化范围很小，说明系统很稳定，应该是参数设置的问题，你再调调

chunshui2003

回复 12 # hsfy919 的帖子

谢谢你的提醒。经你这么一说，我再看看其他相关的文献，发现这个轴心轨迹的变化范围确实是有点小了，可能就是参数设置的问题。感谢！

gghhjj

这个模型很熟悉，好像是上交的吧

chunshui2003

回复 14 # gghhjj 的帖子

恩，是的，成玫老师在振动工程学报上发表的。我的模型和它的不一样，只是某些地方有些相似，拿过来说明一下情况。