请大家帮忙分析系统处于何种状态(附分岔、轨迹和poincare图)

chunshui2003

x1方向分岔图

以上分别为cz(x1，y1)的分岔图，cz=0.00100,   0.00230，0.00550时的poincare映射图和轨迹图。
这是去掉稳态后200个周期的结果。
我的问题是：
1.     cz=0.00100时，由于poincare映射图出现了一系列不规则的点，但是轨迹图又好像蛮规则的。结合分岔图是否可以认为此时系统处于混沌状态。
2.     cz=0.00230时，虽然poincare映射图出现了一系列不规则的点，但是如果考虑到小数点后第3位的话，是否可以认为映射图上实际为1点（但好像又同分岔图相互矛盾）。
3.     cz=0.00550时，其情况同cz=0.00230时类似。结合分岔图、轨迹图和映射图，在考虑精度后，此时系统是应该处于何种运动状态呢？周期运动、周期-n运动、还是混沌呢？

有些糊涂了，请大家帮忙分析一下。

chunshui2003

写帖子的时候已经调整好图片的大小了，但是不知道为什么上传后变成这么大。另外，前几张图片也不清楚为什么没有文件名，而最后四张我本来是没有粘贴进去的，但是也显示出来了并且还附带名称。晕了！

hsfy919

回复 2 # chunshui2003 的帖子

从你以上的图中应该说不能完全确定系统处于什么状态，我不知道你研究的是什么故障，但好像轨迹的范围变化不大，有些点虽然很多但坐标变化不大，有可能其实是一个点，也就是周期1，有可能是精度不够造成的，另外分岔图太过凌乱，如果愿意把你的程序贴上来吧，可以让大家帮你分析一下

chunshui2003

回复 3 # hsfy919 的帖子

今天我也考虑到精度上的问题了。之前计算时一直采用默认绝对和相对误差，今天把误差值变换后就得到了不同的结果。另外，分岔图凌乱指的是什么意思呢？点太多了还是其他的原因。程序很简单，贴出来麻烦给分析一下。

chunshui2003

我现在所学的非线性动力学知识几乎有80%都是来自于振动论坛，所以没有什么所谓保留。即使我知道了新的方法或者知识，也会同大家分享，因为我觉得只有这样才能够共同进步。麻烦大家看一下。

1. 
   
2. 
   
3. tic
   
4. clc
   
5. clear
   
6. global cz
   
7. period = 2*pi;               
   
8. tspan = (0:period/100:500*period);
   
9. y0 = [0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001];
   
10. 
    
11.                   %%%%%%%%%%%%%%%%         参数变化范围          %%%%%%%%%%%%%%%
    
12. %                                                                          
    
13.                             cz_sub = 0.20*10^-3 : 0.005*10^-3 : 6.0*10^-3;                                             
    
14.                            
    
15. %
    
16. %
    
17. row = length(cz_sub);                                       
    
18. column = length(40100:100:50001);                                               
    
19. 
    
20. U_x1 = zeros(row,column);                                   %  转子轴心x1的ode计算结果存储
    
21. U_y1 = zeros(row,column);                                   %  转子轴心y1的ode计算结果存储
    
22. U_plus_cz_x1 = zeros(row,column+1);                         %  x1和cz的混合结果存储
    
23. U_plus_cz_y1 = zeros(row,column+1);                         %  y1和cz的混合结果存储
    
24. %
    
25.                  %%%%%%%%%%%%%%%%%      数值计算      %%%%%%%%%%%%%%%%%%%
    
26. 
    
27. for z = 1:length(cz_sub)
    
28.    cz = cz_sub(z);
    
29.    disp(cz)
    
30.    [T,u] = ode15s(@equ_elastic_without_UMP,tspan,y0);
    
31.    y0 = u(end,:);                                          %将最后的计算结果作为下一次的初值
    
32.    
    
33.    U_x1(z,:) = u(40100:100:end,1)';                                             
    
34.    U_y1(z,:) = u(40100:100:end,3)';                                             
    
35. end
    
36. 
    
37.                  %%%%%%%%%%%%%%%%%      计算结果存储      %%%%%%%%%%%%%%%%%%%
    
38. 
    
39. %
    
40. U_plus_cz_x1(:,1) = cz_sub';
    
41. U_plus_cz_x1(:,2:end) = U_x1;
    
42. U_plus_cz_y1(:,1) = cz_sub';
    
43. U_plus_cz_y1(:,2:end) = U_y1;
    
44. %
    
45. m = row;
    
46. n = column + 1;
    
47. 
    
48.                 %%%%%%%%%%%%%%%%       计算结果写入       %%%%%%%%%%%%%%%%
    
49.                
    
50.                
    
51. %      x1的分岔结果写入
    
52. 
    
53. fid = fopen('长径比=0.20,e0=0.0006,忽略UMP时cz_x1分岔结果_Situ_E.dat','w+');                        
    
54. 
    
55. for i = 1:m
    
56.     for j= 1:n
    
57.         fprintf(fid,'%6.4e        ',U_plus_cz_x1(i,j));
    
58.     end
    
59.     fprintf(fid,'\n');
    
60. end
    
61. fclose(fid)
    
62. %      y1的分岔结果写入
    
63. 
    
64. fid = fopen('长径比=0.20,e0=0.0006,忽略UMP时cz_y1分岔结果_Situ_E.dat','w+');                     
    
65. 
    
66. for i = 1:m
    
67.     for j= 1:n
    
68.         fprintf(fid,'%6.4e        ',U_plus_cz_y1(i,j));
    
69.     end
    
70.     fprintf(fid,'\n');
    
71. end
    
72. fclose(fid)
    
73. 
    
74. save cz_change_Situ_E                                                                           
    
75. 
    
76. toc
    
77. 
    

复制代码

%%%%%%%%%%%%微分方程%%%%%%%%%%%%%

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

复制代码

chunshui2003

然后绘制分岔图就是hsfy919 之前所告诉的方法，在origin中画，这个不好叙述，不过很简单，把计算结果导入到origin即可。

hsfy919

回复 6 # chunshui2003 的帖子

我大概看了一下，你的“参数变化范围”那部分程序是什么意思，另外你把积分方法换成ode45试试，不好意思，我现在要出去一下，一会回来仔细看看

chunshui2003

回复 7 # hsfy919 的帖子

参数变化范围指的就是以cz为控制参数。cz_sub是cz的替代品，因为采用的是全局变量，所以将cz_sub的值附给cz即可，效果是一样的。
ode45算的时间非常久，10分钟后也没有结果，我的方程大概是刚性的，所以采用ode15s计算，7秒左右可以算一次。

gghhjj

从你给的程序和计算结果来看，初步可以判断计算没有收敛
1. 从庞加莱图看，大部分状态下基本上是一条直线（你可以调整一下坐标的尺度就能看出来），现在纵横坐标的尺度不一致
2. 从计算程序所给的参数看，个人初步感觉非线性现象未必存在

chunshui2003

回复 9 # gghhjj 的帖子

学长你好。非常感谢你的回答，其中有几个问题不太清楚。
1 你说计算结果没有收敛是从哪里判断出的呢？分岔图还是计算周期？如果是从分岔图角度讲，如何判断，其达到收敛的标志是什么？
2“ 映射图几乎是直线”这个我能理解，如果将纵坐标的精度控制一下就能达到。不过你说的横纵坐标尺度不一指的是数量级不同吗？
3 我起初计算的时候是采用初值不变的原则，只绘制出系统的轴心轨迹，发现非常混乱（稳态后），按照一般按概念应该可能会出现拟周期或者混沌现象。学长所说的从参数看，系统的非线性现象未必存在的判据是否是根据经验得来？能否告之一下。

hsfy919

回复 10 # chunshui2003 的帖子

在非线性油膜力作用下，系统会表现出自激振动特征，从你的分叉图上看不出任何明显特征，建议你从原始方程入手，查找建模的合理性、程序输入的准确性，参数选取的是否合适等等。我有一次好像也出现了解不收敛的情况，弄了好几天发现是方程输错了，郁闷

chunshui2003

回复 11 # hsfy919 的帖子

多谢，我在查看一下，也许是方程的问题。

chunshui2003

回复 11 # hsfy919 的帖子

听了你的建议之后我又重新推导了一次方程，发现同之前的一样，所以方程没有错误；然后又逐个核对程序的符号、参数，发现也没有问题。至于建立模型的合理性，如果你有时间的话可以看一下“成玫，孟光，荆建平. 非线性“转子-轴承-密封”系统动力分析”.振动工程学报，2006，19(4)：519-524上面的模型。我的模型和她的相似，只不过我的属于立式机组。
那么只剩下参数的选取了，我觉得我选择的可能不合适。
你说非线性油膜力会使系统表现出自激振动特征。我想问一下，如果在参数选取合适的情况下“自激振动”效应在分岔图上应该表现出如何的特征呢，能描述一下吗。这样我就能判断参数的选取得到的结果是否合适了。谢谢。

gghhjj

油膜力自激振动在分岔图上一般是随着转速升高先发生二倍周期分岔（对应半频涡动，频率略低于1/2转速频率，0.45-0.48左右），到一定转速后发生hopf分岔（对应油膜振荡）
对于立式转子而言，经常会出现直接发生hopf分岔的现象

chunshui2003

回复 14 # gghhjj 的帖子

感谢gghhjj学长的解答。我今天就按照您说的再试一下改变转速再做分岔图。
另外，学长，再向您请教一个问题：在选取分岔控制参数的时候，是否一般选取同激励项有关的参数。看了很多相关的文献，主要是偏心、转速，如果涉及到轴系不对中时，还有平行不对中量和偏角，根据学长的经验，如果将轴承间隙cz作为控制参数（油膜力项，也属于外激励），能否出现非线性现象呢。
我的无量纲化是对delta_0进行的，delta_0=4.5*10^-3m。之前经过hsfy919的提醒，提到轴心轨迹的范围太小，正负0.02之前，如果换算成单位制为0.02*4.5*10^-3m，确实很小。看到文献中的无量纲结果大部分都在1附近，那是否说明结构参数的选择不合理亦或是无量纲操作出现了问题（对于x1,y1，无量纲操作为delta_0；对于x2,y2，无量纲操作为cz。是否这方面也有影响）。
学长，一下提了这么多问题，不好意思，真是迫切想解决问题，麻烦你了。