newmark for rotor.m

用matlab法求解椭圆轴承支承的刚性jeffcott转子的振动方程。采用capone模型计算轴承油膜力。得到圆盘中心的轨迹图。

clc;
clear all;

m=1; %转子质量
r=57e-3;%轴颈半径 
c=0.2e-3;%轴承半径间隙 
yt=18e-3 ;%润滑油黏度 
L=28.5e-3;%轴瓦宽度
%  kp=1.052e8;c1=4000;c2=9000;
w=1000;%转速
e=0.06e-3;%质量偏心
g=9.8;%重力加速度
D=2*r;%轴颈直径
fai=c/r;%间隙比
pi=3.14159265;
%%%%%%%%%%%%%%无量纲化
G=g/(c*w^2);%无量纲g
P=e/c;%无量纲质量偏心率
MM=(m*c*w*fai^2)/(yt*L*r);%无量纲质量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%x=zeros(2,1);
%x1=zeros(2,1);
%x2=zeros(2,1);
%x0=[0.8,0.8,0.8,0.8]';
%dt=pi/100;

u=[-0.8638,0.2512]';%位移初始值
v=[-0.2039,-0.7229]';%速度初始值
a=[0,0]';%加速度初始值
%t(1)=0;               %时间
x(:,1)=u;             %位移
x1(:,1)=v;            %速度
x2(:,1)=a;            %加速度


 
%M=zeros(2);
K=zeros(2);
C=zeros(2);
%Q=zeros(2,1);
%F=zeros(2,1);
M=[MM,0;
    0,MM];


%newmark参数设置
gama=0.5;
dt=pi/100;
delta=0.25;
a0=1/(delta*dt^2);
a1=gama/(delta*dt);
a2=1/(delta*dt);
a3=1/(2*delta)-1;
a4=gama/delta-1;
a5=dt*(gama/(2*delta)-1);
a6=dt*(1-gama);
a7=gama*dt; 
%等效刚度矩阵
Ke=K+a0*M+a1*C;
t_max=1*pi;       %计算时间总长
i=1;
t(1)=0;
%q=zeros(2,1);
%fid1=fopen('weiyi.dat','wt');
while t(i)<t_max 
 %%%%%%广义力设定
%Qx=w^2*ww*ee*cos(t(i));
%Qy=w^2*ww*ee*sin(t(i));%此时不考虑重力
tao=w*t(i);%无量纲时间
Qx=MM*P*sin(tao);
Qy=MM*P*cos(tao)+MM*G;
Q=[Qx,Qy]';
%Q(i)=[MM*P*sin(tao(i)),MM*P*cos(tao(i))+MM*G]';
arf=atan((x1(1,i)+2*x(2,i))/(x(1,i)-2*x1(2,i)))-pi/2*sign((x1(1,i)+2*x(2,i))/(x(1,i)-2*x1(2,i)))-pi/2*sign(x1(1,i)+2*x(2,i));
Ga=2/sqrt(1-x(1,i)^2-x1(1,i)^2)*(pi/2+atan((x1(1,i)*cos(arf)-x(1,i)*sin(arf))/sqrt(1-x(1,i)^2-x1(1,i)^2)));
Va=(2+(x1(1,i)*cos(arf)-x(1,i)*sin(arf))*Ga)/(1-x(1,i)^2-x1(1,i)^2);
Sa=(x(1,i)*cos(arf)+x1(1,i)*sin(arf))/(1-(x(1,i)*cos(arf)+x1(1,i)*sin(arf))^2);
f11=-sqrt((x(1,i)-2*x1(2,i))^2+(x1(1,i)+2*x(2,i))^2)/(1-x(1,i)^2-x1(1,i)^2);
f12=3*x(1,i)*Va-Ga*sin(arf)-2*cos(arf)*Sa;
f22=3*x1(1,i)*Va+Ga*cos(arf)-2*sin(arf)*Sa;
 fx=f11*f12;
 fy=f11*f22;
Fx=-fx;
Fy=-fy;
F=[Fx,Fy]';
    %%%%%%%%%%%%%%%%%%%%%
    %Q=[0 0 fx 0 0 0 0 0 fy 0 0 0]';
      q(:,i+1)=F+Q+M*(a0*x(:,i)+a2*x1(:,i)+a3*x2(:,i))+C*(a1*x(:,i)+a4*x1(:,i)+a5*x2(:,i)); 
      x(:,i+1)=inv(Ke)*q(:,i+1);
      x2(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*x1(:,i)-a3*x2(:,i);
%       x1(i,i+1)=a1*(x(i,i+1)-x(i,i))-a4*x1(i,i)-a5*x2(i,i);
      x1(:,i+1)=x1(:,i)+a6*x2(:,i)+a7*x2(:,i+1);
%       fprintf(fid1,'%5d %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e\n',i,x(i,i+1));
  plot(t(i),x(1,i),'k*');
%       plot(t,x1(1,i),'ro');
%       hold on;
% plot(t,x2(1,i),'go');
% hold on;
      i=i+1;
      t(i)=t(i-1)+dt;
 %       plot(t(i),x(1,i));
  end
 %   figure(1) 
  %plot(t(:),x(1,:),'k*');   
  %  plot(t(i),x(1,i));
      
  %  figure(2)
  %  plot(t,x(7,i));
  %  figure(3)
  %  plot(x(1,i),x(7,i));
  %   hold on;
% fclose(fid1);