带油膜力碰摩转子非线性特性分析 程序出不了分岔图 望指导

boom

function zdot=pengmo(tt,z,flag,w)
zdot=zeros(8,1);


m1=4.0;%转子在轴承处集中质量kg
m2=32.1;%转子在圆盘处等效的集中质量kg
b1=0.00005;%偏心量m
R=0.025;%轴承半径m
L=0.012;%轴承长度m
c=0.00011;%平均油膜厚度m
u=0.018;%润滑油黏度N*s/m
%w=900;%转子转动角速度
f=0.1;
c1=1050;%转子在轴承出阻尼系数N*S/m
c2=2100;%转子圆盘轴承系数N*s/m
k=2.5*10^7;%静子刚度N/m
g=9.8;%N/kg
%delta=0.005;%m
delta=0.000016;
k1= 3.5*10^6;
G1=g/(c*w^2);
b=b1/c;

tau=w*tt;
P=m2*g/2;%转子圆盘重量的一半
s=(u*w*R*L*((R/c)^2)*(L/(2*R))^2)/P; % s为Sommerfeld修正系数，c平均油膜厚度
e=sqrt(z(5)^2+z(7)^2);

if e>=delta
Px=-c*((e-delta)/e)*k1*(z(5)-f*z(7));
Py=-c*((e-delta)/e)*k1*(z(5)*f+z(7));

else
Px=0;
Py=0;
end
   

alpha=atan((z(3)+2*z(2))/(z(1)-2*z(4)))-pi/2*sin((z(3)+2*z(2))/(z(1)-2*z(4)))-pi/2*sin(z(3)+2*z(2));
S=(z(1)*cos(alpha)+z(3)*sin(alpha))/(1-abs((z(1)*cos(alpha)+z(3)*sin(alpha))^2));
G=(2/abs(sqrt(1-abs(z(1)^2)-abs(z(3)^2))))*(pi/2+atan((z(3)*cos(alpha)-z(1)*sin(alpha))/abs(sqrt(1-abs(z(1)^2)-abs(z(3)^2)))));
V=(2+(z(3)*cos(alpha)-z(1)*sin(alpha))*G)/abs((1-abs(z(1)^2)-abs(z(3)^2)));

fx=-abs((abs((z(1)-2*z(4))^2)+abs((z(3)+2*z(2))^2))^(1/2))/(1-abs(z(1)^2)-abs(z(3)^2))*(3*z(1)*V-G*sin(alpha)-2*S*cos(alpha));
fy=-abs((abs((z(1)-2*z(4))^2)+abs((z(3)+2*z(2))^2))^(1/2))/(1-abs(z(1)^2)-abs(z(3)^2))*(3*z(3)*V-G*cos(alpha)-2*S*sin(alpha));%油膜力

zdot(1)=z(2);
zdot(2)=-c1*z(2)/(w*m1)-k*(z(1)-z(5))/(w^2*m1)+s*P*fx/(w^2*m1*c)+b*cos(tau);
zdot(3)=z(4);
zdot(4)=-c1*z(4)/(w*m1)-k*(z(3)-z(7))/(w^2*m1)+s*P*fy/(w^2*m1*c)+b*sin(tau)-G1;
zdot(5)=z(6);
zdot(6)=-c2*z(6)/(w*m2)-2*k*(z(5)-z(1))/(w^2*m2)+Px/(c*m2*w^2);
zdot(7)=z(8);
zdot(8)=-c2*z(8)/(w*m2)-2*k*(z(7)-z(3))/(w^2*m2)+Py/(c*m2*w^2)-G1;

clc
clear
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=0:25:1000;           %转子的转速是以rad/s为单位的转速


%位移、速度、转角、转角速度的初值
z1=1e-005;   z2=1e-005;   z3=-1e-5;     z4=3e-005;      z5=2e-005;   z6=-3e-5;    z7=-1e-005;   z8=1e-005;     

for  n=1:length(w);
      T=2*pi;
      ts=[0:T/200:T*100];                          % 时间区间/w(n)
      z0=[z1;z2;z3;z4;z5;z6;z7;z8;];       % xo包含的三十个向量           
      [tt,z]=ode45('pengmo',ts,z0,[],w(n));
   
      figure(1)
      plot(w(n),z(100:100:end,1),'*');
      xlabel('\fontsize{18}\omega');
      ylabel('\fontsize{18}x');grid
      hold on
   end
toc

wqpower

@boom，你好。貌似你的程序中没有对系统运动方程进行无量纲化处理，有可能是这个原因。

杨一枭

你的这个问题弄出来吗？