lj.m

内燃机转子仿真

clear all;
clc;
format long;
%在400——1100rpm转速下,从0.5－10谐次的 振动幅值，扭矩，应力 及合成值
close all
we=[0.2369269 0.4786287 0.5688889 0.4786287 0.2369269];
xi=[-0.9061798 -0.538469 0 0.538469 0.9061798];
tu1=0;
ts=1;
n=10;
f1=zeros(n,1);
f2=zeros(n,1);
aa0=0;
o2=zeros(n,1);
jj=zeros(n,n);
kk=zeros(n,n);
ca=zeros(n,n);
cr=zeros(n,n);
Mgg=zeros(n,20);
Mgi=zeros(n,20);
Mg=zeros(n,20);
E=[9.8790 2.6700 8.1140 5.9840 8.1140 8.1140 5.9840 8.1140 64.6100 470.8800];       
gd=[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];
ABS_DAMP=[0 0 230 230 230 230 230 230 0 0];    
REIA_DAMP=[2700 0 0 0 0 0 0 0 0 0];            
DIAM=[0.1 210 210 210 210 210 210 300 220];   
lamd=0.2544;
llc=[0.25*lamd -0.5 -0.75 -0.25*lamd^2];
stroke=29;   
r=stroke/2; 
d=28;     
speed=1000;
w=2*pi*speed/60;
jiaodu=[360 600 120 480 240 0]*pi/180; 
delt=jiaodu/w;
h=3*60/360/speed;
remt1=rem(ts,h);
remt2=ts-remt1;
tn=remt2/h;

for i=1:n
    jj(i,i)=E(i);
end



for i=1:n
    ca(i,i)=ABS_DAMP(i);            
end

for i=2:n
    cr(i,i)=REIA_DAMP(i-1)+REIA_DAMP(i);
end
for i=1:n-1
    cr(i,i+1)=-REIA_DAMP(i);
    cr(i+1,i)=-REIA_DAMP(i);
end
cr(1,1)=REIA_DAMP(1);
c=ca+cr;  


M=jj;
C=c;
    MEP=9.80665*1e4*20.153*(speed*60/2/pi/1000)^0;        

load('qPA6.dat');

for i=1:20            
     re(i)=(qPA6(2*i-1,1)*MEP^3+qPA6(2*i-1,2)*MEP^2+qPA6(2*i-1,3)*MEP+qPA6(2*i-1,4))*1e-4;   
     im(i)=(qPA6(2*i,1)*MEP^3+qPA6(2*i,2)*MEP^2+qPA6(2*i,3)*MEP+qPA6(2*i,4))*1e-4;      
end

cci=sqrt(re.^2+im.^2);
atanf=atan(re./im);



  F=zeros(2*n,1);
  N=20;
  N1=20;
  a=pi/6;
  dt=h/2^N;
 f=zeros(n,1);
  kn=0;
 vk=zeros(2*n,1);

 xyz(:,1)=vk(1:n,1);
xyz1(:,1)=vk(1:n,1);

for i=1:tn
     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

        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:n)=gd(2:n)-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;
%%%%%%%%%%%%%%%%%%%%%%% &&&me 
 sumf7=20.153;
for j=0.5:0.5:10
    sumf7=sumf7+cci(2*j)*sin(j*200*pi*t(:,i)/6+atanf(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:2
      f1(ix,1)=0;
  end
      
      sumf1=aa0;
      sumf2=aa0;
      sumf3=aa0;
      sumf4=aa0;
      sumf5=aa0;
      sumf6=aa0;
     
      for iw=0.5:0.5:10
    
sumf1=sumf1+cci(2*iw)*sin(iw*200*pi/6*(t(:,i)+delt(1)+0.5*h+0.5*h*xi(ii))+atanf(2*iw));
sumf2=sumf2+cci(2*iw)*sin(iw*200*pi/6*(t(:,i)+delt(2)+0.5*h+0.5*h*xi(ii))+atanf(2*iw));
sumf3=sumf3+cci(2*iw)*sin(iw*200*pi/6*(t(:,i)+delt(3)+0.5*h+0.5*h*xi(ii))+atanf(2*iw));
sumf4=sumf4+cci(2*iw)*sin(iw*200*pi/6*(t(:,i)+delt(4)+0.5*h+0.5*h*xi(ii))+atanf(2*iw));
sumf5=sumf5+cci(2*iw)*sin(iw*200*pi/6*(t(:,i)+delt(5)+0.5*h+0.5*h*xi(ii))+atanf(2*iw));
sumf6=sumf6+cci(2*iw)*sin(iw*200*pi/6*(t(:,i)+delt(6)+0.5*h+0.5*h*xi(ii))+atanf(2*iw));
       

end
f1(3,1)=sumf1;
f1(4,1)=sumf2;
f1(5,1)=sumf3;
f1(6,1)=sumf4;
f1(7,1)=sumf5;
f1(8,1)=sumf6;    

for ix=9:10;
    f1(ix)=0;
end

%vf=vf+we(ii)*T1*(F0*sin(w*(t(:,i)+0.5*h+0.5*h*xi(ii)));
F1(:,1)=[o2;f1]*pi*d^2/4;
for ix=3:8
    f2(ix,1)=0;
for iw=1:4
    f2(ix,1)=f2(ix)+73.5*r^2*w^2*llc(iw)*sin(iw*w*(t(:,i)+delt(ix-2)+0.5*h+0.5*h*xi(ii)))/10000;
end
end
F2(:,1)=[o2;f2]    ;
F(:,1)=F1+F2;

Ffff(:,i)=F1;    
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)*gd(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','tn','ts','h');
save('neiranji1','xyz1','t','tn','ts','h');
save('neiranji2','Txyz','t','tn','ts','h');
 
 
%wn12=pi*(0.0001)^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=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;



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

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轴段变化曲线');


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轴段变化曲线');

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轴段变化曲线');


end