palls.m

内燃机转子仿真

%PIM methold 

close all;
clear;
close all;
format long;
we=[0.2369269 0.4786287 0.5688889 0.4786287 0.2369269];
xi=[-0.9061798 -0.538469 0 0.538469 0.9061798];
zhuansu=1000;
f0=zeros(16,1);
if9=[0 20 200 220 40 60 160 180 80 100 120 140];

tu1=1;
Ts=0.2;                                
d=60.0/zhuansu/360*3.0;               
TD=rem(Ts,d);
TT=Ts-TD;
nn=TT/d;                            

CT1=1;                                    
CTD1=rem(CT1,d);
CTT1=CT1-CTD1;
Cnn1=CTT1/d;                    

CT2=5.5;                                     
CTD2=rem(CT2,d);
CTT2=CT2-CTD2;
Cnn2=CTT2/d;                     

CT3=6.5;                                         
CTD3=rem(CT3,d);
CTT3=CT3-CTD3;
Cnn3=CTT3/d;                     

n=16;
f=zeros(n,1);
o1=zeros(n,480);
o2=zeros(n,1);
xyz=zeros(n,nn);
xyz1=zeros(n,nn);
Txyz=zeros(n,nn);
xy=zeros(n,nn);

countx2=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]';
x=countx2;
temp=diag(eye(2*n));

guanliang=[9.8790 2.6700 8.1140 5.9840 8.1140 8.1140 5.9840 8.1140 64.6100 470.8800]; 
gangdu=[0.1 0.679*1e6 0.2096*1e9 0.3413*1e8 0.3413*1e8 0.3413*1e8 0.3413*1e8 0.3413*1e8 0.66*1e8 0.981*1e8 0.1];
bgd=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
%adamp=[0 0 0 5.1 5.1 5.1 5.1 5.1 5.1 0 0 0 0 0 0 0 ];
%rdamp=[4700 0 0 3078  3078 3078 3078 3078 0 0 0 0 0 0 0 ];
adamp=[0 0 230 230 230 230 230 230 0 0]; 
rdamp=[2700 0 0 0 0 0 0 0 0 0]; 


dm=zeros(n,n);
dc=zeros(n,n);
for i=1:n
   dm(i,i)=guanliang(i);
end
for i=2:n-1
   dc(i,i)=rdamp(i-1)+rdamp(i);
end
dc(1,1)=rdamp(1);
dc(n,n)=rdamp(n-1);
for i=1:n-1;
   dc(i,i+1)=-rdamp(i);
   dc(i+1,i)=-rdamp(i);
end
for i=1:n
   dc(i,i)=dc(i,i)+adamp(i);
end


M=dm;
C=dc;



w=zhuansu*2*pi/60;


mj=73.5;

lamd=0.2544;
stroke=0.29;  
R=stroke/2; 
D=0.28;    

load data11
load PA6.dat
alph=PA6(:,1);
  
pressure1=PA6(:,2);


alph=alph*pi/180;



 for i=1:480
     
 if i>240
         ii=rem(i,240)+1;
  else
       ii=i;
   end 
   
   beta(i)=asin(lamd*sin(alph(ii)));
   
   Mg81(i)=(pi*(D^2)*R*pressure1(ii)/4)*sin(alph(ii)+beta(i))/cos(beta(ii));
 
  Mg82(i)=mj*(R^2)*(w^2)*(0.25*lamd*sin(alph(ii))-0.5*sin(2*alph(ii))-0.75*lamd*sin(3*alph(ii))-0.25*(lamd)^2*sin(4*alph(ii)));
  
   Mg8(i)=Mg81(i)+Mg82(i); 
 end
 %Mgx8(1:480)=Mg8(480:719);



 for i=1:480
     
   for ii=1:12
      ik=i+if9(ii);
      if(ik>240)
         ik=rem(ik,240)+1;
      end
      Mg(ii)=Mg8(ik);
        
   end
   for ii=1:6
       f(ii+3)=Mg(ii*2-1)+Mg(ii*2);
   end
   f(1)=0;
   f(2)=0;
   f(3)=0;
   f(15)=0;
   f(16)=0;
   for jj=10:14
       f(jj)=0;
   end
 
      
af(:,i)=f;
 end
%######
%xx=0:0.2:nn;
%yy=spline(1:nn,Mg8(:),xx);



 %********************击振力***************




  t(1,1)=0;
  F1=zeros(32,1);
  N=20;
  N1=20;
  a=pi/6;
  h=d;
  dt=h/2^N;
  w=2*pi*zhuansu/60;
  kn=0;
 vk=zeros(32,1);
 xyz(:,1)=vk(1:n,1);
xyz1(:,1)=vk(1:n,1);
AF=zeros(32,480);
AF=[o1;af];

   for i=1:nn 
 if rem(i,480)==0
     ij=1;
 else ij=rem(i,480);
 end
       
t(:,i)=(i-1)*h;

    oo=w*t(:,i);
   
    if oo<=2*kn*pi+0.5*pi-a&oo>=2*kn*pi

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

        ff(:,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

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

        ff(:,i)=0.5*(1-cos((0.5*(oo-3*pi/2+a)/a)*pi));
    
    elseif 2*kn*pi+3*pi/2+a<oo<=2*pi+2*kn*pi
         ff(:,i)=1;
    end
   
    kt(2:16)=gangdu(2:16)-bgd*ff(:,i);
    kt(1)=0.1;
    kt(n+1)=0.1;
    for ii=1:n
   K(ii,ii)=kt(ii)+kt(ii+1);
    end

%K(1,1)=kt(1);


%K(n,n)=kt(n-1);
for ii=2:n;

  K(ii-1,ii)=-kt(ii);
   K(ii,ii-1)=-kt(ii);
end

   H=[-inv(M)*C/2 inv(M);0.25*C*inv(M)*C-K -0.5*C*inv(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;
   for jj=1:N
    Ta=2*Ta+Ta*Ta;
  end
  T=I+Ta;
 
 if i<=Cnn1
       AF(31,ij)=42975*i*d;
       AF(32,ij)=42975*i*d;
   elseif i<=Cnn2
       AF(31,ij)=42975;
       AF(32,ij)=42975;
   elseif i<=Cnn3
       AF(31,ij)=42975*(6.5-i*d);
       AF(32,ij)=42975*(6.5-i*d);
   else
      AF(31,ij)=0;
      AF(32,ij)=0;
      end
  
ww(1:n,1)=w;
vv(1:n,1)=w*i*d;
%AF(:,ij)=AF(:,ij)-[o2;C*ww];

%if ij==1
    %F1=af(:,240);
%elseif ij==0
  %  F2==af(:,1);
%end 
%F0=[o1;f];


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=2*Ta1+Ta1*Ta1;
  end
  T1=I+Ta1;
if ij==1
    if ii==1|ii==2
  vf=vf+we(ii)*T1*(AF(:,1)+xi(ii)*(AF(:,1)-AF(:,479)));
elseif ii==4|ii==5
  vf=vf+we(ii)*T1*(AF(:,ij)+xi(ii)*(AF(:,ij+1)-AF(:,ij)));
else
    vf=vf+we(ii)*T1*AF(:,ij);
    end
else
   if ii==1|ii==2
  vf=vf+we(ii)*T1*(AF(:,ij)+xi(ii)*(AF(:,ij)-AF(:,ij-1)));
elseif ii==4|ii==5
  vf=vf+we(ii)*T1*(AF(:,ij)+xi(ii)*(AF(:,ij+1)-AF(:,ij)));
else
    vf=vf+we(ii)*T1*AF(:,ij);
    end
    end
  %vf=vf+we(ii)*T1*F0*sin(w*(t(:,i)+0.5*h+0.5*h*xi(ii)));

  %vf=vf+we(ii)*T1*(F1+(F0-F1)*(0.5+0.5*x(ii)));
 
end
vk=T*vk+0.5*h*vf; 
  % if (i>1)%在 这里是表示在每一步计算的时候都把角度的值清零，但是还考虑到速度和加速度的影响 
     % countx2=x;
     %end
    temp=vk;
     x=vk(1:n,1);
     %vk=vk-w*i/2400;
     Txyz(1:n,i)=x-w*i/2400;
     
     xy(1:n,i)=x;
     zhuanzhix=x';
     
          for ll=1:(n-1)
         xiangduix(ll)=zhuanzhix(ll)-zhuanzhix(ll+1);
         xdnj(ll)=xiangduix(ll)*gangdu(ll);
     end
    
     
   
     xyz(2:n,i)=xdnj';
    xyz1(2:n,i)=xiangduix;
if oo>2*(kn+1)*pi
        kn=kn+1;
end
  

   end     
save('data11','ff0','ff1');
save('neiranji','xyz','t','nn','Ts','h');
save('neiranji1','xyz1','t','nn','Ts','h');
save('neiranji2','Txyz','t','nn','Ts','h');
save('neiranji3','xy','t','nn','Ts','h');
%%%%%%%%%%%%%%%%%%%%%%%% 15个轴段应力  %%%%%%%%%%%%%%%%%%%%%%%%%%%
wn12=pi*(0.1)^3/16;
wn23=pi*(0.21)^3/16;
wn34=pi*(0.21)^3/16;
wn45=pi*(0.21)^3/16;
wn56=pi*(0.21)^3/16;
wn67=pi*(0.21)^3/16;
wn78=pi*(0.21)^3/16;
wn89=pi*(0.3)^3/16;
wn910=pi*(0.22)^3/16;

xyz12=xyz(2,:);
yingli12=xyz12/wn12/1e6;
maxyl12=max(yingli12);
minyl12=min(yingli12);
feng12=maxyl12-minyl12;

xyz23=xyz(3,:);
yingli23=xyz23/wn23/1e6;
maxyl23=max(yingli23);
minyl23=min(yingli23);
feng23=maxyl23-minyl23;

xyz34=xyz(4,:);
yingli34=xyz34/wn34/1e6;
maxyl34=max(yingli34);
minyl34=min(yingli34);
feng34=maxyl34-minyl34;

xyz45=xyz(5,:);
yingli45=xyz45/wn45/1e6;
maxyl45=max(yingli45);
minyl45=min(yingli45);
feng45=maxyl45-minyl45;

xyz56=xyz(6,:);
yingli56=xyz56/wn56/1e6;
maxyl56=max(yingli56);
minyl56=min(yingli56);
feng56=maxyl56-minyl56;

xyz67=xyz(7,:);
yingli67=xyz67/wn67/1e6;
maxyl67=max(yingli67);
minyl67=min(yingli67);
feng67=maxyl67-minyl67;

xyz78=xyz(8,:);
yingli78=xyz78/wn78/1e6;
maxyl78=max(yingli78);
minyl78=min(yingli78);
feng78=maxyl78-minyl78;

xyz89=xyz(9,:);
yingli89=xyz89/wn89/1e6;
maxyl89=max(yingli89);
minyl89=min(yingli89);
feng89=maxyl89-minyl89;

xyz910=xyz(10,:);
yingli910=xyz910/wn910/1e6;
maxyl910=max(yingli910);
minyl910=min(yingli910);
feng910=maxyl910-minyl910;



%%%%%%%%%%%%%%%%%%%%%%曲线%%%%%%%%%%%%%

if tu1==0
figure(1)
plot(t,yingli12,'g');title('第12轴段变化曲线');


figure(2)
plot(t,yingli23,'g');title('第23轴段变化曲线');


figure(3)
plot(t,yingli34,'g');title('第34轴段变化曲线');
figure(4)
plot(t,yingli45,'r');title('第45轴段变化曲线');
figure(5)
plot(t,yingli56,'g');title('第56轴段变化曲线');
figure(6)
plot(t,yingli67,'r');title('第67轴段变化曲线');
figure(7)
plot(t,yingli78,'g');title('第78轴段变化曲线');
figure(8)
plot(t,yingli89,'r');title('第89轴段变化曲线');

elseif tu1==1
figure(1)
plot(t,xyz(2,:),'g');title('第12轴段变化曲线');


figure(2)
plot(t,xyz(3,:),'g');title('第23轴段变化曲线');


figure(3)
plot(t,xyz(4,:),'g');title('第34轴段变化曲线');
figure(4)
plot(t,xyz(5,:),'r');title('第45轴段变化曲线');
figure(5)
plot(t,xyz(6,:),'g');title('第56轴段变化曲线');
figure(6)
plot(t,xyz(7,:),'r');title('第67轴段变化曲线');
figure(7)
plot(t,xyz(8,:),'g');title('第78轴段变化曲线');
figure(8)
plot(t,xyz(9,:),'r');title('第89轴段变化曲线');

elseif tu1==2
figure(1)
plot(t,xy(1,:),'g');title('第12轴段变化曲线');


figure(2)
plot(3*t/d,xy(2,:),'g');title('第23轴段变化曲线');


figure(3)
plot(3*t/d,xy(3,:),'g');title('第34轴段变化曲线');
figure(4)
plot(3*t/d,xy(4,:),'r');title('第45轴段变化曲线');
figure(5)
plot(3*t/d,xy(5,:),'g');title('第56轴段变化曲线');
figure(6)
plot(3*t/d,xy(6,:),'r');title('第67轴段变化曲线');
figure(7)
plot(3*t/d,xy(7,:),'g');title('第78轴段变化曲线');
figure(8)
plot(3*t/d,xy(8,:),'r');title('第89轴段变化曲线');
figure(9)
plot(3*t/d,xy(9,:),'r');title('第910轴段变化曲线');

elseif tu1==3
    figure(1)
plot(t,xyz1(2,:),'g');title('第12轴段变化曲线');


figure(2)
plot(t,xyz1(3,:),'g');title('第23轴段变化曲线');


figure(3)
plot(t,xyz1(4,:),'g');title('第34轴段变化曲线');
figure(4)
plot(t,xyz1(5,:),'r');title('第45轴段变化曲线');
figure(5)
plot(t,xyz1(6,:),'g');title('第56轴段变化曲线');
figure(6)
plot(t,xyz1(7,:),'r');title('第67轴段变化曲线');
figure(7)
plot(t,xyz1(8,:),'g');title('第78轴段变化曲线');
figure(8)
plot(t,xyz1(9,:),'r');title('第89轴段变化曲线');


end