pim12.m

内燃机转子仿真

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];
%k=[2 1 0;1 3 0;0 0 4];
%F=[2;4;6];
m=[2 0;0 1];
c=[1.2 -0.4;-0.4 0.6];
k=[6 -2;-2 4];
F=[0;10];
%[M K C]=arryf(m,k,c,n);
t=0;
w=pi;
dk=zeros(n);
ci=2;
dk(ci,ci)=0;
tf=1;
h=0.05;
%h=pi/40;
N=20;
N1=20;
dt=h/2^N;
v=zeros(2*n+1,1);
v(1,1)=t;
%vk=diag(eye(2*n))*0;
vk=[0;0;0;0];
door=dk*(1+cos(w*t))/2;
o1=zeros(n,1);
F0=[o1;F];
H=[-inv(m)*c/2 inv(m);0.25*c*inv(m)*c-k -0.5*c*inv(m)];
[Hx Hn]=eig(H);
for i=1:tf/h
t(:,i)=(i-1)*h;
v(1,i)=t(:,i);
%door=(1+cos(w*t(:,i)))/2;
    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;
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%求总Ta
    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;    
end
 % PIM end
 save('pimdata','v');
 figure(1);
plot(t,v(2,:),t,v(3,:));