liewenjingque.m

内燃机转子仿真

%function dy=liewenjingque(t,y)
clear all;
clc;
close all;
format long;
n=1;
c=0.0;
m=2;
k=20;
dk=4;
q0=0;
q1=1;
f=0;
%dy=zeros(2,1);
a=pi/6;
h=0.01;
t=10;
w=2*pi*30;
bt=rem(t,h);
bbt=t-bt;
n1=bbt/h;
v=zeros(n+2,n1);
v(:,1)=[0;q0;q1];
    kn=0;
for i=2:n1
    t(:,i)=(i-1)*h;
    oo=w*t(:,i)/60;  
    if oo<=2*kn*pi+0.5*pi-a&oo>=2*kn*pi
         ff=1;
        dff=0;
    elseif oo<=2*kn*pi+0.5*pi+a&oo>2*kn*pi+0.5*pi-a
         ff=0.5*(1+cos((0.5*(oo-0.5*pi+a)/a)*pi));
        dff=(-w*pi/240/a)*sin((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
         ff=0;
        dff=0;
    elseif oo<=2*kn*pi+3*pi/2+a&oo>2*kn*pi+3*pi/2-a
         ff=0.5*(1-cos((0.5*(oo-3*pi/2+a)/a)*pi));
        dff=-(w*pi/240/a)*sin((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
        ff=1;
         dff=0;
    end
     if oo>2*(kn+1)*pi
        kn=kn+1;
     end
  wn=(k-dk*ff)^0.5/m^0.5;
 dwn=dk*dff/2/(m*(k-dk*ff))^0.5;
 y=q0*cos(wn*t(:,i))+q1*sin(wn*t(:,i))/wn;
dy=-q0*sin(wn*t(:,i))*(dwn*t(:,i)+wn)+(q1*cos(wn*t(:,i))*(dwn*t(:,i)+wn)*wn-q1*sin(wn*t(:,i))*dwn)/wn^2 ;
v(1,i)=t(:,i);
v(2:3,i)=[y;dy];
q0=y;
q1=dy;
end  
figure(1);
plot(v(1,:),v(2,:));   
figure(2);
plot(v(1,:),v(3,:));