tailei01.m

内燃机转子仿真

clear all
clc
close all
mm=[5.98000 1.08000 1.04000 2.91300 2.91300 2.91300 2.91300 2.91300 51.46300 0.60000 1.11500 3.94400];
cc=[0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000];
kk=[0.00000 8.20000 392.20000 150.00000 112.78000 112.78000 112.78000 112.78000 112.78000 169.66000 0.50000 0.50000 50.29000 0.00000];
FF=[0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000];
%n=1;
n=length(FF);
m=zeros(n);
c=zeros(n);
k=zeros(n);
for i=1:n
    m(i,i)=mm(i);
end
kk=kk*1.0e5;
for i=1:n
    k(i,i)=kk(i)+kk(i+1);
    c(i,i)=cc(i)+cc(i+1);
end
for i=1:n-1
            k(i,i+1)=-kk(i+1);
            c(i,i+1)=-cc(i+1);
            k(i+1,i)=-kk(i+1);
            c(i+1,i)=-cc(i+1);
end
F=FF';

S=inv(m)*k;
[px pn]=eig(S);

%w=0;
x=[0.500 0.500 1.000 1.200 0.500 0.500 1.000 1.200 0.500 0.500 1.000 1.200];
x=x';
y=zeros(n,1);
z0=inv(m)*F;

%w=1200;
%dt=60/w/360*2;
h=pn(1,1)/10;
t=6;
%t1=1;
dn=rem(t,h);
dn1=t-dn;
n1=dn1/h;

%k1=0.5;
%dn1=rem(t1,dt);
%dn2=t1-dn;
%n2=dn2/dt;
f=zeros(n,1);
df=zeros(n,1);
dx=zeros(n,1);
dy=zeros(n,1);
d2xx=zeros(n,1);
d2yy=zeros(n,1);
d3x=zeros(n,1);
d3y=zeros(n,1);
xyz1=zeros(2*n+1,n1);
xyz1(2:n+1,1)=x;
xyz1(n+2:2*n+1,1)=y;
for i=2:n1
    f=F;
    df=0;
    d2f=0;
    dx=y;
    dy=inv(m)*(f-c*y-k*x);
    d2x=dy;
    d2y=inv(m)*(df-c*dy-k*y);
    d3x=d2y;
    d3y=inv(m)*(d2f-c*d2y-k*dy);
    x=x+h*(dx+0.5*h*(d2x+h*d2x/3));
    y=y+h*(dy+0.5*h*(d2y+h*d3y/3));
    if (0.8<abs(y(1,1))<1.0|0.8<abs(y(2,1))<1.0)
        x=xyz1(2:n+1,i-1);
        y=xyz1(n+2:2*n+1,i-1);
        h1=h/30;
        %i=i-1;
        for j=1:30
            dx=y;
            dy=inv(m)*(f-c*y-k*x);
            d2x=dy;
            d2y=inv(m)*(df-c*dy-k*y);
            d3x=d2y;
            d3y=inv(m)*(d2f-c*d2y-k*dy);
            x=x+h1*(dx+0.5*h1*(d2x+h1*d2x/3));
            y=y+h1*(dy+0.5*h1*(d2y+h1*d3y/3));
        end
    elseif (0.6<abs(y(1,1))<=0.8|abs(0.6<y(2,1))<=0.8)
        x=xyz1(2:n+1,i-1);
        y=xyz1(n+2:2*n+1,i-1);
        h2=h/90;
        for j=1:90
            dx=y;
            dy=inv(m)*(f-c*y-k*x);
            d2x=dy;
            d2y=inv(m)*(df-c*dy-k*y);
            d3x=d2y;
            d3y=inv(m)*(d2f-c*d2y-k*dy);
            x=x+h2*(dx+0.5*h2*(d2x+h2*d2x/3));
            y=y+h2*(dy+0.5*h2*(d2y+h2*d3y/3));
        end
        elseif (0.4<abs(y(1,1))<=0.6|0.4<abs(y(2,1))<=0.6)
        x=xyz1(2:n+1,i-1);
        y=xyz1(n+2:2*n+1,i-1);
        h3=h/270;
        for j=1:270
            dx=y;
            dy=inv(m)*(f-c*y-k*x);
            d2x=dy;
            d2y=inv(m)*(df-c*dy-k*y);
            d3x=d2y;
            d3y=inv(m)*(d2f-c*d2y-k*dy);
            x=x+h3*(dx+0.5*h3*(d2x+h3*d2x/3));
            y=y+h3*(dy+0.5*h3*(d2y+h3*d3y/3));
        end
        elseif (0.2<abs(y(1,1))<=0.4|abs(0.2<y(2,1))<=0.4)
        x=xyz1(2:n+1,i-1);
        y=xyz1(n+2:2*n+1,i-1);
        h4=h/380;
        for j=1:380
            dx=y;
            dy=inv(m)*(f-c*y-k*x);
            d2x=dy;
            d2y=inv(m)*(df-c*dy-k*y);
            d3x=d2y;
            d3y=inv(m)*(d2f-c*d2y-k*dy);
            x=x+h1*(dx+0.5*h4*(d2x+h4*d2x/3));
            y=y+h1*(dy+0.5*h4*(d2y+h4*d3y/3));
        end
        elseif (0.1<abs(y(1,1))<=0.2|0.1<abs(y(2,1))<=0.4)
        x=xyz1(2:n+1,i-1);
        y=xyz1(n+2:2*n+1,i-1);
        h5=h/520;
        for j=1:520
            dx=y;
            dy=inv(m)*(f-c*y-k*x);
            d2x=dy;
            d2y=inv(m)*(df-c*dy-k*y);
            d3x=d2y;
            d3y=inv(m)*(d2f-c*d2y-k*dy);
            x=x+h5*(dx+0.5*h5*(d2x+h5*d2x/3));
            y=y+h5*(dy+0.5*h5*(d2y+h5*d3y/3));
        end
        elseif (abs(y(1,1))<=0.1|abs(y(2,1))<=0.1)
        x=xyz1(2:n+1,i-1);
        y=xyz1(n+2:2*n+1,i-1);
        h6=h/640;
        for j=1:640
            dx=y;
            dy=inv(m)*(f-c*y-k*x);
            d2x=dy;
            d2y=inv(m)*(df-c*dy-k*y);
            d3x=d2y;
            d3y=inv(m)*(d2f-c*d2y-k*dy);
            x=x+h6*(dx+0.5*h6*(d2x+h6*d2x/3));
            y=y+h6*(dy+0.5*h6*(d2y+h6*d3y/3));
        end
    end
    xyz1(1,i)=i;
    xyz1((2:n+1),i)=x;
    xyz1((n+2):(2*n+1),i)=y;
end
figure(1);
plot(xyz1(1,:)*h,xyz1(2,:));
figure(2);
plot(xyz1(1,:)*h,xyz1(3,:));