pimsdo.m

内燃机转子仿真

%function []=pim1x(bf)
clear all;
clc;
close all;
format long;
n=1;
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];
m=2;
K=8;
dk1=0;
%c=0.2*(m*K)^0.5;
c=0.0;
F=4;
t=0;
a=pi/6;
fz=(K/m)^0.5;
fn=0.8;
w=1;
tf=10;
h=0.01;
%h=Tin/10;
N=20;
N1=20;
dt=h/2^N;
v=zeros(2*n+1,1);
v(1,1)=t;
vk=[0;0];
%vk=[-1;-3;-2;1;2;3]
o1=zeros(n,1);
F0=[0;F];
%[Hx Hn]=eig(H);
kn=0;
for i=1:tf/h
t(:,i)=(i-1)*h;
v(1,i)=t(:,i);
%door=(1+cos(w*t(:,i)))/2;
oo=w*t(:,i);

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*dk1;
kti(:,i)=kt;
Fi(:,i)=F0*sin(w*(i-1)*h);
H=[-c/2/m 1/m;0.25*c^2/m-kt -0.5*c/m];
[Hx Hn]=eig(H);
    I=eye(size(H));
 Ta=H*dt+(H*dt)^2*(I+(H*dt)/3+(H*dt)^2/12+(H*dt)^3/60)/2;%泰勒展开高阶项
 %Ta=H*dt+(H*dt)^2*(I+(H*dt)/3)/2;
  for j=1:N%求总Ta
    Ta=2*Ta+Ta*Ta;
  end
  T=I+Ta;
v(2:2*n+1,i)=vk(:);%给定初值
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*cos(w*(t(:,i)+0.5*h+0.5*h*xi(ii)));
end
vk=T*vk+0.5*h*vf;    
if oo>2*(kn+1)*pi&w>0
    kn=kn+1;
end
end
 % PIM end
 d1k=dk1(1);
 vx2=v;
 save('sjjl','v','m','fn','door1','K','dk1','h');
 figure(1);
plot(t,v(2,:));
title('一质量振动时域波形图');
xlabel('时间(s)');
ylabel('振幅');
if w>0
figure(2);
plot(t,door1,'k');
xlabel('时间(s)');
ylabel('幅值');
title('开关函数');
end
 figure(3);
plot(t,kti(:),'k');
ylim([0,22]);
xlabel('时间(s)');
ylabel('刚度');
title('刚度随时间的变化');
%figure(3);
%plot(t,v(4,:));