%非线性刚度不平衡转子动力学行为研究

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%系统方程（圆盘不平衡）
%    x'' + 2ξx' + x +γ( x.^2 + y.^2) x = e*ω_bar.^2cosτ
%    y'' + 2ξy' + y +γ( x.^2 + y.^2) y = e*ω_bar.^2sinτ - G
%    G =g/(ω.^2*δ) ;τ=ω*t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%初始参数
%    2ξ = 0. 24
%    ω_bar = 2. 1
%      γ = 0. 5
%偏心量e :[0. 08 ～0. 18 ]
%e = 0. 1175、0. 1384 和0. 1750
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
keci=0.12;
omega_bar=2.1;
gama=0.5;
delta=0.1;
omega=2.0;
G=9.8/(omega.^2*delta);
% e=0.08:0.01:0.18;
t=0:500;

tao=omega*t;
initial=[0,0,0,0];
global coef;
cf=[];
ccf=[];

% e=[ 0.1175, 0.1384 ,0.1750];
e=0.1750
% for i=1:length(e)
%     if e(i)==0.1175 || e(i)==0.1384||e(i)==0. 1750
        coef=[2*keci,1,gama,e*omega_bar.^2,-G];       
%        for t=0:0.1:500;%500个周期（omega）
            [t,x]=ode45('right',tao,initial,[],coef);
            cf=[cf;x];
%        ccf=[ccf;cf];
%   end
% end

last_300=length(cf(:,1));
last_300=last_300-300;
range=last_300:length(cf(:,1));
figure

plot(cf(range,2),cf(range,1),'b',cf(range,4),cf(range,2),'g')
title('PointCare 图');
xlabel('位移')
ylabel('速度')

figure
plot(cf(range,1),cf(range,3),'r')
title('轴心轨迹图');
xlabel('x')
ylabel('y')

figure
xfft=fft(cf(range,1))
xfft_Amplit=abs(xfft);
subplot(211)
plot(xfft_Amplit)
title('x方向图频谱');
xlabel('f')
ylabel('Ampli')

subplot(212)
plot(cf(range,1))
title('x方向时域图形');
xlabel('t')
ylabel('x')

 

 

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f=right(t,x)
global coef;

f=[x(2);
   -x(2)*coef(1)-coef(2)*x(1)-coef(3)*(x(1).^2+x(3).^2)*x(1)+coef(4)*cos(t);
   x(4);
   -x(4)*coef(1)-coef(2)*x(3)-coef(3)*(x(1).^2+x(3).^2)*x(1)+coef(4)*sin(t)-coef(5);  
]
%