关于滚动轴承转子系统的分叉图问题

冬夜 · 发表于 2012-4-24 23:56

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

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

x

大家好，我看了有关滚动轴承支承下转子系统动力学研究的文章，自己画得轴心轨迹对着，但分叉图老是不对，不知道什么原因，大家能否帮忙看下。结合http://forum.vibunion.com/thread-24346-1-1.html里面的程序，大家能看的清晰点（但是这里面没有进行无量纲化）下面是我的程序及画出的分叉图：
子程序：

function dy=jeff_fun118(t,y,flag,W)
%相关参
%m=0.6;NB=9;d0=45*10^(-6);
%R=0.038/2;r=0.028/2;
%M为圆盘质量，NB为轴承滚珠数目，d0为游隙，w为角速度
% gear bearing
% syms sita(),deltasita(),j,s1,s2,E(3),ma,l(e),k(e),kp,k(x),k(y);
%sita是每一个滚珠经历t后的角度，deltasita是第j个滚珠的弹性接触变形.
%Y=y/d;
%dy/dt=w*dx/dT;
%d(dy/dt)/dt=w^2*d(dy/dT)/dT;
%T=w*t;
%W=w*(d0/g)^(1/2);
%R1=R/d;r1=r/d;cc=c/(m*(g/d0)^(1/2));
m=0.6;ff=6;C=340;NB=9;a=0;R1=0.0265/2;r1=0.0173/2;d0=45*10^(-6);g=9.8;S1=0;S2=0;E=0.005;
cc=(C/m)*sqrt(d0/g);Kp=7.055*(10^9)*d0^(3/2)/(m*g);fff=ff/(m*g);ee=E/d0;
for j=1:1:NB;
Sita(j)=2*pi*(j-1)/NB+t*r1/(r1+R1);
Deltasita(j)=(y(1)*cos(Sita(j))+y(3)*(sin(Sita(j))))*cos(a)-1; %deltasita(j)=(y(1)*cos(sita(j))+y(3)*(sin(sita(j))))*cos(a)-d;
if Deltasita(j)<=0
Deltasita(j)=0;
end
Sa=((Deltasita(j)).^(3/2))*cos(Sita(j)); %sa=(deltasita(j)).^(3/2)*cos(sita(j)); sa=d^(3/2)*Sa
Sb=((Deltasita(j)).^(3/2))*sin(Sita(j)); %sb=(deltasita(j)).^(3/2)*sin(sita(j)); sb=d^(3/2)*Sb
S1=Sa+S1;
S2=Sb+S2;
end
Qx=Kp*S1;
Qy=Kp*S2;
t(end,1)
dy=[y(2);-cc*y(2)/W+ee*cos(t)-Qx./W^2+fff/W^2;y(4);-cc*y(4)/W+ee*sin(t)-Qy./W^2];
主程序：
clc
clear
d0=45*10^(-6);g=9.8;
f=(100*2*pi/60)*(d0/g)^(1/2):0.05:(20000*2*pi/60)*(d0/g)^(1/2); % 转速参数
k=0;
for i=1:length(f)
disp(f(i))
period=2*pi;
b=f(i); %转速的参数传递transmagic
k=k+1;
step=period/100;
tspan=0:step:500*period;
options=odeset('RelTol',1.0e-9);
y0=[0.1;0.2;0.1;0.1];
[t,y]=ode45('jeff_fun118',tspan,y0,options,b);
plot(f(i),y(45000:100:end,1),'r.','markersize',1) %画分叉图
pause(1);
hold on
end

复制代码

伤痕累累 · 发表于 2012-4-25 09:06

你这是那篇论文的啊。

冬夜 · 发表于 2012-4-25 10:16

回复 2 # 伤痕累累的帖子

滚动轴承－转子系统的非线性动力学研究.西北工业大学硕士论文

冬夜 · 发表于 2012-4-25 12:54

没人研究这个吗？

伤痕累累 · 发表于 2012-4-25 14:20

回复 4 # 冬夜的帖子

你的轴心轨迹图是对的？？？不对吧？我试了一下。跟原文差距比较大

冬夜 · 发表于 2012-4-26 17:13

回复 5 # 伤痕累累的帖子

里面有些东西我不会处理，滚动轴承系统其实是一种非线性参数和强迫两者联合的振动，程序里面的周期我是根据强迫力的周期进行计算，我画出来的轴心轨迹是一个椭圆呀，（她的看上去很乱）上面连接里的那个程序人家是用非线性参数算的，如果我的程序改成他那样的话，周期就不是2pi了，还跟转速有关了，这样画出来的分叉图也很怪，我不知道哪里出问题了,还需大家多多帮忙！

伤痕累累 · 发表于 2012-4-27 15:10

回复 6 # 冬夜的帖子

慢慢来。还得多看看书和文献。

伤痕累累 · 发表于 2012-4-27 15:11

回复 6 # 冬夜的帖子

可以先进行无量纲处理。让周期变成2*pi,然后用频闪法作出poincare截面。如果正确，下一步就是分叉图。

冬夜 · 发表于 2012-4-27 15:48

回复 8 # 伤痕累累的帖子

恩好吧我再多看看书

dongyuanzhixing · 发表于 2012-5-8 17:15

看着就蒙了……

ME! · 发表于 2013-7-30 20:23

伤痕累累发表于 2012-4-27 15:11

登录/注册后可看大图

回复 6 # 冬夜的帖子

可以先进行无量纲处理。让周期变成2*pi,然后用频闪法作出poincare截面。如果正确，下 ...

请问一下楼主，有没有用matlab画转子的轴心轨迹啊，我想请教一下
我对转子施加了力F(t)=F0sinwt，F(t)=F0coswt，然后根据徐斌的《matlab有限元结构动力学》求响应再画节点的轨迹，画出来是个圆，请问这对吗？无阻尼的

q0=zeros(length(K),1);
dq0=zeros(length(K),1); % 初始化位移和速度
eta0=Vnorm'*M*q0; deta0=Vnorm'*M*dq0; %2.2-53% 初始条件模态坐标的位移和速度
eta=zeros(nstep,sdof); % %初始化位移2.2-52
phase0=omega0*tt; %w*t omega0为激振力频率
for i=1:sdof % t（i）时刻的响应
gama=omega0/omega(i); %omega0为激振力的圆频率，omega(i)=w,系统的无阻尼固有频率
omegad=omega(i)*sqrt(1-zeta(i)^2);
phase=omegad*tt; %wd*t其中(omegad为有阻尼固有频率)
Exx=exp(-zeta(i)*omega(i)*tt); %中间变量
C1=eta0(i); %初始位移
C2=(deta0(i)+eta0(i)*zeta(i)*omega(i))/omegad; %中间变量
X0=sqrt((1-gama^2)^2+(2*zeta(i)*gama)^2);
XX=Fnorm(i)/(omega(i)^2*X0);
XP=atan((2*zeta(i)*gama)/(1-gama^2));
D1=(zeta(i)*omega(i)*cos(XP)+omega0*sin(XP))/omegad;
D2=cos(XP);
eta(:,i)=C1*Exx.*cos(phase)+C2*Exx.*sin(phase)...%系统对初始条件的响应
-XX*Exx.*(D1*sin(phase)+D2*cos(phase))...%伴随自由振动
+XX*cos(phase0); % 稳态响应
%eta(:,i)=XX*cos(phase0-XP); %稳态强迫振动
end
%------------------------------------------------------
eta=eta'; %模态坐标
y=Vnorm*eta; %物理坐标

复制代码

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

[分形与混沌] 关于滚动轴承转子系统的分叉图问题

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