shiyuxieci.m

内燃机转子仿真

%
%function 

clear;
clc;
close all;
format long;
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];
AA=42975;
aa=63;
zhuansu=1000;
f0=zeros(16,1);
if9=[0 20 200 220 40 60 160 180 80 100 120 140];

tu1=0;
Ts=0.6;                                
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);


ctt=(1:2:91);
ctt1=ctt';
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=[19.627 7.059 3.049 12.566  9.61  12.566 12.566 9.61 12.566 42.11 49.09 177.6 140.5 11.759585 308.155203 310.113323];
gangdu=[0.000001 0.651 29.86 32.75 32.75 32.75 32.75 32.75 32.75 325 0.316 205 9.7768 11.75483 3.11507 11.10479 0.000001]*1.0e6;
%E=[19.627 9.599 0.509 13.018  9.908  13.018 13.018 9.908 13.018 42.11 49.09 177.6 14007.5 11.759585 308.155203 310.113323];
%F=[0.651 39.86 32.75 32.75 32.75 32.75 32.75 32.75 325 0.316 205 9.7768 11.75483 3.11507 11.10479]*1.0e6;
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 ];


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.07;

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
 for i=1:240
     for j=1:12
    Mg(j,i)=Mg8(i+if9(j)) ;
     end
     for ii=1:6
         Mgg(ii,i)=Mg(2*ii-1,i)+Mg(2*ii,i);
     end
 end
%Mg8(241)=0; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%&&&&me
 tq=(0:239)*4*pi/240;
 for i=1:6
aa0(i,1)=sum(Mgg(i,:))/240 ;
 end
aa1=zeros(24,1);

aa2=zeros(24,1);
for i=1:24
    aa1(i,1)=2*sum(Mgg.*cos(i*tq/2))/240;
    aa2(i,1)=2*sum(Mgg.*sin(i*tq/2))/240;
end
%aa1()=sum(aa1);
%aa2=sum(aa2);
cc=sqrt(aa1.^2+aa2.^2);
tanff=tan(aa1./aa2);
tanf=atan(tanff);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%tt(1)=0;
tt(1:1001)=(0:d:1000*d);
ffm1=zeros(1,length(tt));

for i=1:1001
  sums1=aa0; 
  ccct=aa0;
for j=0.5:0.5:12
    sums1=sums1+aa1(j*2,1)*cos(j*200*pi*tt(i)/6)+aa2(j*2,1)*sin(j*200*pi*tt(i)/6);%1000r/500T
    ccct=ccct+cc(2*j)*sin(j*200*pi*tt(i)/6+tanf(2*j));
   
    
end
ffm1(1,i)=sums1;
ffm2(1,i)=ccct;
end
%ffm1(1,1)=0;
%ffm2(1,1)=0;
cci=(0.5:0.5:12);
cci=cci';
%ccct=aa0+sum(cc*sin(cci*200*pi*tt(i)/6+tanf));
%for ii=1:240
%f4(:,ii)=aa0+sum(cc.*sin(cci*200*pi/6*tt(ii)+tanf));
%end
ctt=(1:2:46)*0.48332193495134;

figure(1);
plot(tt(1:1001),ffm1(1:1001));

figure(2);
plot(tt(1:1001),ffm2(1:1001));

%figure(3);
%plot(tt,f4);
%kkk(1,:)=Mg8;
%kkk(2,:)=ffm1(1:240);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%&&&me
  w=2*pi*zhuansu/60;
t(1,1)=0;
  F1=zeros(32,1);
  N=20;
  N1=20;
  a=pi/6;
  h=d;
  dt=h/2^N;

  kn=0;
 vk=zeros(2*n,1);

 xyz(:,1)=vk(1:n,1);
xyz1(:,1)=vk(1:n,1);
tt0=if9*3*pi/180;
 
    
 


for i=1:nn
    
    t(:,i)=(i-1)*h;
     
    oo=w*t(:,i);
  % vk(1:n)=vk(1:n)+oo;

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

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

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

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

        ffi(:,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
         ffi(:,i)=1;
    end
   
    kt(2:n)=gangdu(2:n)-bgd*ffi(:,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;
%%%%%%%%%%%%%%%%%%%%%%% &&&me 
 sumf7=aa0;
for j=0.5:0.5:12
    sumf7=sumf7+cc(2*j)*sin(j*200*pi*t(:,i)/6+tanf(2*j));
    %ccct=sums1+cc(2*j)*sin(j*200*pi*tt(i)/6+tanf(2*j));
end
Ffff2(:,i)=sumf7;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%&&&me

%v(:,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;
  for ix=1:3
      f(ix)=0;
  end
      
      sumf1=0;
      sumf2=0;
      sumf3=0;
      sumf4=0;
      sumf5=0;
      sumf6=0;
      sumf7=0;
      sumf8=0;
      sumf9=0;
      sumf10=0;
      sumf11=0;
      sumf12=0;
      for iw=0.5:0.5:12
sumf1=sumf1+cc(2*iw)*sin(iw*200*pi/6*(t(:,i)+0.5*h+0.5*h*xi(ii))+tanf(2*iw));
sumf2=sumf2+cc(2*iw)*sin(iw*200*pi/6*(t(:,i)+0.5*h+0.5*h*xi(ii))+tanf(2*iw));
sumf3=sumf3+cc(2*iw)*sin(iw*200*pi/6*(t(:,i)+0.5*h+0.5*h*xi(ii))+tanf(2*iw));
sumf4=sumf4+cc(2*iw)*sin(iw*200*pi/6*(t(:,i)+0.5*h+0.5*h*xi(ii))+tanf(2*iw));
sumf5=sumf5+cc(2*iw)*sin(iw*200*pi/6*(t(:,i)+0.5*h+0.5*h*xi(ii))+tanf(2*iw));
sumf6=sumf6+cc(2*iw)*sin(iw*200*pi/6*(t(:,i)+0.5*h+0.5*h*xi(ii))+tanf(2*iw));


end
ff(1,1)=sumf1;
ff(2,1)=sumf2;
ff(3,1)=sumf3;
ff(4,1)=sumf4;
ff(5,1)=sumf5;
ff(6,1)=sumf6;   

for ix=10:14;
    f(ix)=0;
end
if i<=Cnn3
f(15)=(4*AA/(0.5*pi-aa*pi/180)/pi)*sum(sin(ctt1*(0.5*pi-aa*pi/180)).*sin(ctt1*0.48332193495134*(t(i)+0.5*h+0.5*h*xi(ii)))./ctt1.^2);
f(16)=(4*AA/(0.5*pi-aa*pi/180)/pi)*sum(sin(ctt1*(0.5*pi-aa*pi/180)).*sin(ctt1*0.48332193495134*(t(i)+0.5*h+0.5*h*xi(ii)))./ctt1.^2);
end
%vf=vf+we(ii)*T1*(F0*sin(w*(t(:,i)+0.5*h+0.5*h*xi(ii)));
F=[o2;f];

vf=vf+we(ii)*T1*F;
end
vk(:,1)=T*vk+0.5*h*vf;    

x=vk(1:n,1);
 zhuanzhix=x';

xyz(:,i)=vk(1:n);
v(:,i)=vk;
vv(1:n,i)=inv(M)*(vk(n+1:2*n)-0.5*C*vk(1:n));
for ll=1:(n-1)
         xiangduix(ll)=zhuanzhix(ll)-zhuanzhix(ll+1);
         xdnj(ll)=xiangduix(ll)*gangdu(ll+1);
end
Txyz(2:n,i)=xdnj';
xyz1(2:n,i)=xiangduix*180/pi;    

Ffff1(:,i)=F;

if oo>2*(kn+1)*pi&w>0
    kn=kn+1;
end;

end;




save('neiranji','xyz','t','nn','Ts','h');
save('neiranji1','xyz1','t','nn','Ts','h');
save('neiranji2','Txyz','t','nn','Ts','h');
 
 
wn12=pi*(0.21)^3/16;
wn23=pi*(0.21)^3/16;
wn34=pi*(0.23)^3/16;
wn45=pi*(0.23)^3/16;
wn56=pi*(0.23)^3/16;
wn67=pi*(0.23)^3/16;
wn78=pi*(0.23)^3/16;
wn89=pi*(0.23)^3/16;
wn910=pi*(0.23)^3/16;
wn1011=pi*(0.23)^3/16;%联轴器内外之间的连接轴
wn1112=pi*(0.23)^3/16;%离合器主从之间的连接轴
wn1213=pi*(0.314^4-0.18^4)/0.314/16;%一号传动轴
wn1314=pi*(0.314^4-0.18^4)/0.314/16;%2号传动轴
wn1415=pi*(0.25^4-0.15^4)/0.25/16;%电轴1
wn1516=pi*(0.38^4-0.064^4)/0.38/16;%两个电枢之间

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

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

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

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

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

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

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

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

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

xyz1011=Txyz(11,:);
yingli1011=xyz1011/wn1011/1e6;
maxyl1011=max(yingli1011);
minyl1011=min(yingli910);
feng1011=maxyl910-minyl910;

xyz1112=Txyz(12,:);
yingli1112=xyz1112/wn1112/1e6;
maxyl1112=max(yingli1112);
minyl1112=min(yingli910);
feng91112=maxyl910-minyl910;


xyz1213=Txyz(13,:);
yingli1213=xyz1213/wn1213/1e6;
maxyl1213=max(yingli1213);
minyl1213=min(yingli1213);
feng1213=maxyl1213-minyl1213;

xyz1314=Txyz(14,:);
yingli1314=xyz1314/wn1314/1e6;
maxyl1314=max(yingli1314);
minyl1314=min(yingli1314);
feng1314=maxyl1314-minyl1314;

xyz1415=Txyz(15,:);
yingli1415=xyz1415/wn1415/1e6;
maxyl1415=max(yingli1415);
minyl1415=min(yingli1415);
feng1415=maxyl1415-minyl1415;;

xyz1516=Txyz(16,:);
yingli1516=xyz1516/wn1516/1e6;
maxyl1516=max(yingli1516);
minyl1516=min(yingli1516);
feng1516=maxyl1516-minyl1516;

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

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轴段变化曲线');
figure(9)
plot(t,yingli910,'r');title('第910轴段变化曲线');
figure(10)
plot(t,yingli1011,'r');title('第1011轴段变化曲线');
figure(11)
plot(t,yingli1112,'r');title('第1112轴段变化曲线');

figure(12)
plot(t,yingli1213,'g');title('第1213轴段变化曲线');
figure(13)
plot(t,yingli1314,'r');title('第1314轴段变化曲线');
figure(14)
plot(t,yingli1415,'r');title('第1415轴段变化曲线');
figure(15)
plot(t,yingli1516,'g');title('第1516轴段变化曲线');

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轴段变化曲线');
figure(9)
plot(t,xyz(10,:),'r');title('第910轴段变化曲线');
figure(10)
plot(t,xyz(11,:),'r');title('第1011轴段变化曲线');
figure(11)
plot(t,xyz(12,:),'r');title('第1112轴段变化曲线');    
figure(12)
plot(t,xyz(13,:),'g');title('第1213轴段变化曲线');
figure(13)
plot(t,xyz(14,:),'r');title('第1314轴段变化曲线');
figure(14)
plot(t,xyz(15,:),'r');title('第1415轴段变化曲线');
figure(15)
plot(t,xyz(16,:),'g');title('第1516轴段变化曲线');

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轴段变化曲线');
figure(9)
plot(t,xyz1(10,:),'r');title('第910轴段变化曲线');
figure(10)
plot(t,xyz1(11,:),'r');title('第1011轴段变化曲线');
figure(11)
plot(t,xyz1(12,:),'r');title('第1112轴段变化曲线');  
figure(12)
plot(t,xyz1(13,:),'r');title('第1112轴段变化曲线'); 
figure(13)
plot(t,xyz1(14,:),'r');title('第1112轴段变化曲线'); 
figure(14)
plot(t,xyz1(15,:),'r');title('第1112轴段变化曲线'); 
figure(15)
plot(t,xyz1(16,:),'r');title('第1112轴段变化曲线'); 


end