liewenjingque1.m

内燃机转子仿真

clear all;
clc;
close all;
format long;
n=1;
m=2;
c=0.1;
k=20;
dk=4;
f=6;
w=60*pi;
h=pi/400;
t=16;
st=rem(t,h);
tt=t-st;
nn=tt/h;
x10=0.1;
x20=0.2;
x30=(f-c*x20-k*x10)/m;
a0=1/h^2;
a1=1/2/h;
a2=2*a0;
a3=1/a2;
x0=x10-h*x20+a3*x30;
m1=a0*m+a1*c;
x=zeros(3*n+1,nn);
x(:,1)=[0;x10;x20;x30];
a=pi/6;
kn=0;
for i=2:nn
    t=(i-1)*h;
    oo=w*t/60;

if w>0

   
    if oo<=2*kn*pi+0.5*pi-a&oo>=2*kn*pi

        door=1;
    elseif oo<=2*kn*pi+0.5*pi+a&oo>2*kn*pi+0.5*pi-a

        door=0.5*(1+cos((0.5*(oo-0.5*pi+a)/a)*pi));
    elseif oo<=2*kn*pi+3*pi/2-a&oo>2*kn*pi+0.5*pi+a

        door=0;
    elseif oo<=2*kn*pi+3*pi/2+a&oo>2*kn*pi+3*pi/2-a

        door=0.5*(1-cos((0.5*(oo-3*pi/2+a)/a)*pi));
    %elseif oo>2*k*pi+3*pi/2+a&oo<=2*pi+2*k*pi
    elseif 2*kn*pi+3*pi/2+a<oo<=2*pi+2*kn*pi
         door=1;
    end
door1(:,i)=door;
else 
    door=0;
end
kt=k-door*dk;
    hf=f;
r=hf-(kt-a2*m)*x10-(a0*m-a1*c)*x0;
x11=inv(m1)*r;
x31=a0*(x0-2*x10+x11);
x21=a1*(-x0+x11);
xx10=x11;
xx20=x21;
xx30=x31;

for j=1:10
df=f-m*xx30-c*xx20-kt*xx10;
hf=df;
r=hf-(kt-a2*m)*x10-(a0*m-a1*c)*x0;
x11=inv(m1)*r;
x31=a0*(x0-2*x10+x11);
x21=a1*(-x0+x11);
xx10=xx10+x11;
xx20=xx20+x21;
xx30=xx30+x31;
end
x0=x10;
x10=xx10;
x20=xx20;
x30=xx30;
x(:,i)=[t;xx10;xx20;xx30];
end
save('zhongxinchafen','x');
plot(x(1,:),x(2,:),x(1,:),x(3,:),x(1,:),x(4,:));