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

chunshui2003 · 发表于 2010-12-14 16:12

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

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

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

function uu = equ_elastic_without_UMP(T,u,e0,cz,Lb)
m1 = 60; %转子质量 kg
m2 = 25; %轴颈质量 kg
c1 = 4000; %转子阻尼 N.s/m
c2 = 1200; %轴承阻尼 N.s/m
Ke = 6.2*10^6; %转轴刚度 N/m
% e0 = 0.6*10^-3; %转子质量偏心 m
% Rr = 60*10^-3; %转子半径 m
% Lr = 150*10^-3; %转子长 m
delta_0 = 4.5*10^-3; %均匀气隙大小 m
% cz = 0.2*10^-3; %轴颈间隙 m
Rb = 50*10^-3; %轴承半径 m
% Lb = 20*10^-3; %轴承长 m
mu = 18*10^-3; %绝对润滑油粘度无单位
omega = 13; %大轴旋转角速度
% global cz
sigma = mu * omega * Rb * Lb * (Rb/cz)^2 * (Lb/2/Rb)^2; %Sommerfeld修正系数
A1 = u(7) + 2*u(6);
A2 = u(5) - 2*u(8);
alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1);
E = sqrt(u(5).^2 +u(7).^2); %轴颈轴心轨迹
E_minus = sqrt(1 -E.^2);
B1 = u(7)*cos(alpha) - u(5)*sin(alpha);
B2 = u(5)*cos(alpha) + u(7)*sin(alpha);
G = 2./E_minus * ( pi/2 + atan(B1./ E_minus) );
V = (2 + B1 * G) ./ E_minus.^2;
S = B2 ./ ( 1 - B2.^2 );
F_oil_same = -sqrt( A2.^2 + A1.^2 ) ./ E_minus.^2;
f_x = F_oil_same * ( 3* u(5) *V - sin(alpha) * G - 2 *cos(alpha)*S );
f_y = F_oil_same * ( 3* u(7) *V + cos(alpha) * G - 2 *sin(alpha)*S ); %油膜力的无量纲表达式
fx = sigma * f_x;
fy = sigma * f_y; %非线性油膜力
uu = zeros(8,1);
uu(1) = u(2);
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 ;
uu(3) = u(4);
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 ;
uu(5) = u(6);
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);
uu(7) = u(8);
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);

复制代码

频闪法：

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

复制代码

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

function linshi2
tic
period = 2*pi; %%%因为对时间t采用了无量纲化，cos(oemga*t)变成了cos(T)，所以周期变为2*pi/1，暂时这样
e0 = 0.6*10^-3; Lb = 100*10^-3; %Situ_H_imp
cz = 0.0002:0.0001:0.0060;
options = odeset('Reltol',1e-3,'Abstol',1e-6);
k = 0;
step = period/100;
for ww = 1:length(cz)
y0 = [0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001];
disp(cz(ww))
k = k+1;
tspan = 0:step:800*period;
[t,y] = ode15s(@equ_elastic_without_UMP,tspan,y0,options,e0,cz(ww),Lb);
y0 = y(end,:);
j = 1;
for i = 800:1000
tspan = i*period:step:(i+1)*period;
[t,y] = ode15s(@equ_elastic_without_UMP,tspan,y0,options,e0,cz(ww),Lb);
YY1(k,j)=y(end,1);
j = j+1;
y0 = y(end,:);
end
end
bifdata = YY1(:,end-20:end);
plot(cz,bifdata,'r.','markersize',5);
toc

复制代码

频闪法

liliangbiao法

chunshui2003 · 发表于 2010-12-14 16:28

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

liliangbiao法poincare

频闪法轨迹

频闪法poincare

chunshui2003 · 发表于 2010-12-14 16:32

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

octopussheng · 发表于 2010-12-15 12:42

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

chunshui2003 · 发表于 2010-12-15 14:44

回复 4 # octopussheng 的帖子

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

octopussheng · 发表于 2010-12-15 15:09

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

chunshui2003 · 发表于 2010-12-15 18:45

回复 6 # octopussheng 的帖子

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

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

复制代码

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

hsfy919 · 发表于 2010-12-15 20:35

回复 7 # chunshui2003 的帖子

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

chunshui2003 · 发表于 2010-12-16 08:43

回复 8 # hsfy919 的帖子

转子-轴承系统

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

jgwang · 发表于 2010-12-16 09:16

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

chunshui2003 · 发表于 2010-12-16 15:00

回复 10 # jgwang 的帖子

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

hsfy919 · 发表于 2010-12-16 15:29

回复 9 # chunshui2003 的帖子

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

chunshui2003 · 发表于 2010-12-16 15:44

回复 12 # hsfy919 的帖子

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

gghhjj · 发表于 2010-12-17 09:18

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

chunshui2003 · 发表于 2010-12-18 07:44

回复 14 # gghhjj 的帖子

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

账号		自动登录	找回密码
密码			我要加入

[稳定性与分岔] 请大家看一下，比较哪一种分岔程序更优

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。