方程如下:
(无量纲化后的):p分别取0.1,0.2,0.3,或0.4(无具体顺序)得到以上图形
function dxdt=force(t,x,flag,w)
r=0.21;
L=0.2016; c=0.000699;m=10700;p=0.068;u=0.0178;g=9.8;
% w=ω为转速,u=η为润滑油动力粘度,p为偏心率
a1=c/r; M=m*c*w*a1^2/(L*r*u);
%a1为间隙比,L为轴瓦宽度,c为轴承半径间隙,r为轴承半径
G1=g/c/w^2;%T=w*t ,T为无量纲时间
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign(x(3)+2*x(2));
G=2/sqrt(1-x(1)^2-x(3)^2)*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/sqrt(1-x(1)^2-x(3)^2)));
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cos(a)+x(3)*sin(a))^2); % x(1)=x,x(2)=x',x(3)=y,x(4)=y'
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G)/(1-x(1)^2-x(3)^2);
fx=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(1)*V-sin(a)*G-2*cos(a)*S);
fy=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(3)*V+cos(a)*G-2*sin(a)*S);
dxdt=zeros(4,1);
dxdt=[ x(2);
-fx/M+p*sin(t);
x(4);
-fy/M+p*cos(t)+G1]; |