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);%无量纲gP=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参数设置gamma=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=gamma/delta-1;a5=dt*(gamma/(2*delta)-1);a6=dt*(1-gamma);a7=gamma*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);哪位高手帮忙写下每步什么意思
|
