slpim2x.m

内燃机转子仿真

%m1=0.39,m2=0.066,e1=89.4*1e-7,e2=658*1e7,e3=345.68*1e7,m1+m2=0.456,wn=402.
%58,vn=3839r/min,767.8,255.93
%function [xdx,h]=pim2x(w,iu,II);
clear all;
clc;
close all;
format long;
n=2;
we=[0.2369269 0.4786287 0.5688889 0.4786287 0.2369269];
xi=[-0.9061798 -0.538469 0 0.538469 0.9061798];
%m=[1 1 0;1 2 0;0 0 3];
%c=[1 0 0;0 0 1;0 1 0];
%k1=[2 1 0;1 3 0;0 0 4];
%F=[2;4;6];
%f=1:2:60
%w=2*pi*60;
c1=10;
%dc=0.2;
dc=0.6;
M=[0.39 0;0 0.066];
KK1=9.148461228821310e+003;

%C=[c+0.1 -c;-c c+0.1];
%dk=120;
dk=60;
pw=89*10.7878;
w=2*pi*600/60;
%w  (rad/min);
f1=(1000*60/2/pi)*pw/w;
f2=-(1000*60/2/pi)*pw/w;


t=0;
a=pi/40;


tf=20;

h=0.001;
ddt=rem(tf,h);
dzt=tf-ddt;
n1=dzt/h;

N=20;
N1=20;
dt=h/2^N;
v=zeros(2*n,n1);

vk=zeros(2*n,1);


o1=zeros(n,1);
F=[f1;f2];
F0=[o1;F];
%[Hx Hn]=eig(H);
kn=0;

for i=1:n1
t(:,i)=(i-1)*h;

%door=(1+cos(w*t(:,i)))/2;
oo=w*t(:,i);



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

        door(:,i)=1;
    elseif oo<=2*kn*pi+0.5*pi+a&oo>2*kn*pi+0.5*pi-a

        door(:,i)=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(:,i)=0;
    elseif oo<=2*kn*pi+3*pi/2+a&oo>2*kn*pi+3*pi/2-a

        door(:,i)=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(:,i)=1;
    end


KK=KK1-door(:,i)*dk;
K=[KK+0.1 -KK;-KK KK+0.1];

 a1=pi/20;
     a2=pi/20;
 if oo<=2*kn*pi+0.5*pi-a1&oo>=2*kn*pi

       
         ffc(:,i)=0;
    elseif oo<=2*kn*pi+0.5*pi+a2&oo>2*kn*pi+0.5*pi-a2

        ffc(:,i)=0.5*(1+cos((0.5*(oo-0.5*pi+a2)/a2)*pi));
        
    elseif oo<=2*kn*pi+3*pi/2-a1&oo>2*kn*pi+0.5*pi+a2

       
         ffc(:,i)=0;
    elseif oo<=2*kn*pi+3*pi/2+a1&oo>2*kn*pi+3*pi/2-a1

        ffc(:,i)=0.5*(1-cos((0.5*(oo-3*pi/2+a1)/a1)*pi));
        
    %elseif oo>2*k*pi+3*pi/2+a&oo<=2*pi+2*k*pi
    elseif 2*kn*pi+3*pi/2+a2<oo<=2*pi+2*kn*pi
     
         ffc(:,i)=0;
 end   
 c=c1+dc*ffc(:,i) ;
C=[c+0.01 -c;-c c+0.01]; 
H=[-inv(M)*C/2 inv(M);0.25*C*inv(M)*C-K -0.5*C*inv(M)];
    I=eye(size(H));
 Ta=H*dt+(H*dt)^2*(I+(H*dt)/3+(H*dt)^2/12+(H*dt)^3/60)/2;%泰勒展开高阶项
  for j=1:N%求总Ta
    Ta=2*Ta+Ta*Ta;
  end
  T=I+Ta;
%给定初值
vf=0;
    for ii=1:5
        h1=h*(1-xi(ii))/2;
        dt1=h1/2^N1;
       Ta1=H*dt1+(H*dt1)^2*(I+(H*dt1)/3+(H*dt1)^2/12+(H*dt1)^3/60)/2;
       for j=1:N1%求总Ta1
    Ta1=2*Ta1+Ta1*Ta1;
  end
  T1=I+Ta1;
  vf=vf+we(ii)*T1*F0;
%vf=vf+we(ii)*T1*F0*sin(w*(t(:,i)+0.5*h+0.5*h*xi(ii)));
end
vk=T*vk+0.5*h*vf;  
v(:,i)=vk(:);
xdx1(1,i)=vk(1,1)-vk(2,1);
fjnj(1,i)=KK*xdx1(:,i);
if oo>2*(kn+1)*pi
    kn=kn+1;
end
end
save('nb','xdx1','h','tf','t');
 % PIM end
 x=(xdx1(1,:)-mean(xdx1(1,:)))*180/pi;
 step=h;

sf=1/step;
%t=0:step:80;
%y1=sin(2*pi*10*t);
%y2=sin(2*pi*5*t);
ng=length(x);
%nfft=2^nextpow2(ng);
nfft=2^14;
yy=fft(x,nfft);
pyy=yy.*conj(yy)/nfft;
ff=(0:sf/nfft:sf/2-sf/nfft);%%%%%%%%%%%%%%%%%%%%%%%%
nn=1:length(ff);
plot(ff(nn),pyy(nn),'k');