gneiranjincrackok.m

内燃机转子仿真

%
%function [YLmax YLfeng Tfeng mTi mxdi mixdi]=gpa6l280m10(zhuansu) 
%function []=gneiranjincrackok(kic) 
clear all;
clc;



close all;
format long;
[xi we]=grule(15);
%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;

if9=[360 600 120 480 240 0]/3;

kic=5;
%bcn=1;
exn=0;
tu1=2;
Ts=8;                                
d=60.0/zhuansu/360*3.0;  
h=d/5;
TD=rem(Ts,h);
TT=Ts-TD;
nn=TT/h;                            
         

n=10;
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=[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]; %stiffness(n*m/rad)
%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;
if kic==0;
bgd=[0 0 0 0 0 0 0 0 0];
elseif kic==1;
    bgd=[0 0 0 0 0 6.7155e+005 0 0 0];%0.1
elseif kic==2;
    bgd=[0 0 0 0 0 3.2438e+006 0 0 0];%0.2
elseif kic==3;
    bgd=[0 0 0 0 0 7.0588e+006 0 0 0];%0.3
elseif kic==4;
    bgd=[0 0 0 0 0 1.0825e+007 0 0 0];%0.4
 elseif kic==5;
    bgd=[0 0 0 0 0 1.2495e+007 0 0 0]; %  0.45
end
kcn=5;
adamp=[0 0 230 230 230 230 230 230 0 0];
rdamp=[2700 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
%%%%%%%%%%%%%me
%for i=1:n
   % ac(i,i)=adamp(i);
%end
%%%%%%%%%%%%%%%%%me
M=dm;
C=dc;



w=2*pi*zhuansu/60;


mj=73.5;

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


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


alph=alph*pi/180;



 for i=1:240
     %t=(i-1)*d;
 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))+5*(lamd)^3*sin(5*alph(ii))/32+3*(lamd)^4*sin(6*alph(ii))/32);
  
   Mg8(i)=(Mg81(i)+Mg82(i)); 
 end
Mg8(241)=0; 
xx1=(0:d:240*d);
plot(xx1,Mg8);
jmh=d/50;
xx2=(0:jmh:240*d);

%Mg8s=spline(xx1,Mg8,'*',xx2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%&&&&me
hham=48;
 tq=(0:240)*4*pi/241;
 tq1=(1:241)*4*pi/241;
aa0=sum(Mg8)/241 ;
xxx=w*(0:d:240*d);
aa00=trapz(xxx,Mg8)/4/pi;%%%im
aa1=zeros(24,1);
for iw=0.5:0.5:hham
    yyy=Mg8.*cos(iw*xxx);
aa11(2*iw,1)=trapz(xxx,yyy)/2/pi;%%%im
end
aa2=zeros(24,1);
for iw=0.5:0.5:hham
    yyy=Mg8.*sin(iw*xxx);
aa22(2*iw,1)=trapz(xxx,yyy)/2/pi;%%%im
end
for i=1:hham
    aa1(i,1)=2*sum(Mg8.*cos(i*tq/2))/241;
    aa2(i,1)=2*sum(Mg8.*sin(i*tq/2))/241;
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=aa00; 
  ccct=aa0;
for j=0.5:0.5:24
    sums1=sums1+aa11(j*2,1)*cos(j*w*tt(i))+aa22(j*2,1)*sin(j*w*tt(i));%1000r/500T
    ccct=ccct+cc(2*j)*sin(j*w*tt(i)+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((0:d:240*d),Mg8);

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

  kn=0;
  www(1:n,1)=zhuansu/60;
 vk=zeros(2*n,1);
% vk(n+1:2*n,1)=M*2*pi*www+0.5*C*vk(1:n,1);
 xyz(:,1)=vk(1:n,1);
xyz1(:,1)=vk(1:n,1);
tt0=if9*d;
kkn=0; 
vv0=0;
    
%a=pi/50;
%a1=pi/2;
%a2=pi/20;
k=0;

for i=1:nn
    
        adamp=[0 0 230 230 230 230 230 230 0 0];
    rdamp=[2700 0 0 0 0 0 0 0 0]; 
 
    if kic==0
    brdamp=[0 0 0 0 0 0 0 0 0];
    a=0;
    elseif kic==1
    brdamp=[0 0 0 0 0 120 0 0 0];
    a=pi/40;
    elseif kic==2
    brdamp=[0 0 0 0 0 120 0 0 0];
    a=pi/30;
    elseif kic==3
    brdamp=[0 0 0 0 0 120 0 0 0];
    elseif kic==4
        a=pi/20;
    brdamp=[0 0 0 0 0 120 0 0 0];
    elseif kic==5
    brdamp=[0 0 0 0 0 120 0 0 0];
    a=pi/10;
    end
    
    
    t(:,i)=(i-1)*h;
    %%%%%%%%%%%%%%%%me
     oo=w*t(:,i);
   
    if oo<=2*k*pi+0.5*pi-a&oo>=2*k*pi

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

        ffk(:,i)=0.5*(1+cos((0.5*(oo-0.5*pi+a)/a)*pi));
      
    elseif oo<=2*k*pi+3*pi/2-a&oo>2*k*pi+0.5*pi+a

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

        ffk(:,i)=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*k*pi+3*pi/2+a<oo<=2*pi+2*k*pi
        ffk(:,i)=1;
        
    end

    
 if kic==0
     a1=pi/50;
     a2=pi/50;
 elseif  kic==1
     a1=pi/50;
     a2=pi/50;
 elseif kic==2
     a1=pi/50;
     a2=pi/50;
 elseif kic==3
     a1=pi/50;
     a2=pi/50;
 elseif kic==4
     a1=pi/50;
     a2=pi/50;
 elseif kic==5
     a1=pi/50;
     a2=pi/50;
 end
 if oo<=2*k*pi+0.5*pi-a1&oo>=2*k*pi

       
         ffc(:,i)=0;
    elseif oo<=2*k*pi+0.5*pi+a2&oo>2*k*pi+0.5*pi-a2

        ffc(:,i)=0.5*(1+cos((0.5*(oo-0.5*pi+a2)/a2)*pi));
        
    elseif oo<=2*k*pi+3*pi/2-a1&oo>2*k*pi+0.5*pi+a2

       
         ffc(:,i)=0;
    elseif oo<=2*k*pi+3*pi/2+a1&oo>2*k*pi+3*pi/2-a1

        ffc(:,i)=0.5*(1-cos((0.5*(oo-3*pi/2+a1)/a1)*pi));
        
    %elseif oo>2*k*pi+3*pi/2+a&oo<=2*pi+2*k*pi
    elseif 2*k*pi+3*pi/2+a2<oo<=2*pi+2*k*pi
     
         ffc(:,i)=0;
    end   
    
    %%%%%%%%%%%%%%%%%%%me
  
    kt(2:n)=gangdu(2:n)-bgd*ffk(:,i);
      rdamp(1:n-1)=rdamp(1:n-1)+brdamp*ffc(:,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

for ii=2:n-1
   dc(ii,ii)=rdamp(ii-1)+rdamp(ii);
end
dc(1,1)=rdamp(1);
dc(n,n)=rdamp(n-1);
for ii=1:n-1;
   dc(ii,ii+1)=-rdamp(ii);
   dc(ii+1,ii)=-rdamp(ii);
end
for ii=1:n
   dc(ii,ii)=dc(ii,ii)+adamp(ii);
end


C=dc;


   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*w*t(:,i)+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:15
        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