szbjf1.m

内燃机转子仿真

%function [x1 x2 x3]=zbjf(K M C f b)
clear all
clc
close all
%fni=input('the name of the film:','s');
%fid=fopen(fni,'r');
%n=fscanf(fid,'%d',1);
%wind=fscanf(fid,'%d',1);
%TW=fscanf(fid,'%f',1);
%cn=fscanf(fid,'%d',[2,1]);
%KK=fscanf(fid,'%f',1);
%t=fscanf(fid,'%f',1);%仿真时间
%t1=fscanf(fid,'%f',1);
%sn=fscanf(fid,'%d',1);
%b=fscanf(fid,'%f',1);
%K=fscanf(fid,'%f',[n,n]);
%M=fscanf(fid,'%f',[n,n]);
%C=fscanf(fid,'%f',[n,n]);
%fclose(fid);
M=[2 0;0 1];
C=[0 0;0 0];
K=[6 -2;-2 4];
F=[0;10];
n=length(F);
x10=zeros(n,1);
%dK=zeros(1,180);
x20=zeros(n,1);
f0=zeros(n,1);
x30=inv(M)*(f0-C*x20-K*x10); 

f=zeros(n,180);
dt=60.0/sn/360*2.0;%每步两度时间
td=rem(t,dt);
tdd=t-td;
tn=tdd/dt;
t1d=rem(t1,dt);
t1dd=t1-t1d;
t1n=t1dd/dt;
xx=zeros(3,tn);
yy=zeros(3,tn);
zz=zeros(3,tn);

K1=zeros(n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a0=6/b^2/h^2;
a1=3/b/h;
a2=2*a1;
a3=b*h/2;
a4=a0/b;
a5=-a2/b;
a6=1-3/b;
a7=h/2;
a8=h^2/6;
fni=input('the name of resul:','s');
fid=fopen(fni,'wt');
%开始循环
%for jj=1:10
    for i=1:tn;%2
        i
         if i>180
         in=rem(i,180)+1;
         else
          in=i;
         end
     dK(in)=KK*sin(pi*in*2/180);%裂纹类型不一样，裂纹所引起的刚度变化就不一样
     K(cn(1,1),cn(2,1))=K(cn(1,1),cn(2,1))+dK(in);%非线性刚度

         if i<=t1n
          f(:,i)=TW*i*dt;
          f1=f(:,i)+b*TW*dt;
         else
         f(:,i)=TW;
          f1=(1-b)*f(:,i)+b*TW;
          end


%对K1作ldlt分解
     K1=K+a2*M+a8*C;
     K2=K1;
     k1=zeros(n-1,1);
     k=zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          for j=1:n;
               for i=1:j-1 
               k1(i)=K1(j,i)*K1(i,i);
               end
          k(j)=K1(j,j)-K1(j,1:j-1)*k1(1:j-1);
          K1(j,j)=k(j);
          K1(j+1:n,j)=(K1(j+1:n,j)-K1(j+1:n,1:j-1)*k(1:j-1))/k(j);
          end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     l=tril(K1,0)-diag(diag(K1))*diag(ones(n,1))+diag(ones(n,1));
     D=diag(diag(K1))*diag(ones(n,1));
     fa=f1+M*(a6*x10+a7*x20+2*x30)+C*(a8*x10+2.*x20+a9*x30);
     x11=inv(l')*inv(D)*inv(l)*fa;
     x3=a1*(x11-x10)-a2*x20+a3*x30;
     x2=x20+a4*(x3+x30);
     x1=x10+dt*x20+(x3+2*x30)*a5;
     temp1=x10(1,1);
     x10=x1;
     temp2=x1(1,1);
     x20=x2;
     x30=x3;
     pp=temp2-temp1
     while (norm(pp,2)/norm(temp2,2))>1e-10
      gg=1;
      x30=inv(M)*(f(:,i)-C*x20-K*x10);   
      fa=f1+M*(a6*x10+a7*x20+2*x30)+C*(a8*x10+2*x20+a9*x30);
      x11=inv(l')*inv(D)*inv(l)*fa;
      x3=a1*(x11-x10)-a2*x20+a3*x30;
      x2=x20+a4*(x3+x30);
      x1=x10+dt*x20+(x3+2*x30)*a5;
      temp1=x10(1,1);
      x10=x1;
      temp2=x1(1,1);
      x20=x2;
     %x30=x3;
      pp=temp2-temp1
      gg=gg+1;
      end
ci=i-180;%周期 当i大于240时进入下一个周期
      if (ci<1)
      ci=1;
       else 
       ci=ci;
       end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xx(:,ci)=x1;%xx是一个16*（tn-180）阵  第一个周期只保存当前值，以后每步计算值都要保存，有足够的值保证图形的精度
yy(:,ci)=x2;
zz(:,ci)=x3;
       for j=1:3
             if j==1
             fprintf(fid,'%3.2f   %4.3f    ',(i-1)*dt,x1(j,1));
             else
             fprintf(fid,'%4.3f    ',x1(j,1)); 
             end
       end
       for j=1:3
       fprintf(fid,'%4.3f     ',x2(j,1));
       end
       for j=1:3
             if j<3
             fprintf(fid,'%5.4f      ',x3(j,1));
             else
             fprintf(fid,'%5.4f\n',x3(j,1));   
             end
        end
    end
%end
%对仿真数据进行分析
%时域分析