xiaoci1205sun60.m

内燃机转子仿真

%newmark methold 

close all;
clear;

zhuansu=1000;%柴油机转速

T=8;                                 %计算时间
d=60.0/zhuansu/360*3.0;               
TD=rem(T,d);
TT=T-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;%有16个惯性质量快
x=zeros(n,1);
y=zeros(n,1);
z=zeros(n,1);
f=zeros(n,1);
x1=zeros(n,1);
y1=zeros(n,1);
z1=zeros(n,1);
fa=zeros(n,1);
xx=zeros(n,nn-240);
yy=zeros(n,nn-240);
zz=zeros(n,nn-240);
ppxx=zeros(n,nn);
XDNJ=zeros(n-1,nn-240);
countx2=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';
x=countx2;

%*****质量矩阵DM（10，10）***********
%输入惯量和刚度
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 147.5 11.759585 308.155203 310.113323];%inertia(kg*m**2)
gangdu=[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;;%转动惯量矩阵stiffness(n*m/rad)   %刚度有五个值相同，所以这个为六缸发动机
adamp=[0 0 0 5.1 5.1 5.1 5.1 5.1 5.1 0 0 0 0 0 0 0 ];%绝对阻尼 N*M/RAD/S=====================================
rdamp=[4700 0 0 3078  3078 3078 3078 3078 0 0 0 0 0 0 0 ];% %相对阻尼  N*M/RAD/S=============================================
%形成刚度矩阵和惯量矩阵
dm=zeros(n,n);
dk=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);
   dk(i,i)=gangdu(i-1)+gangdu(i);
end
dc(1,1)=rdamp(1);
dk(1,1)=gangdu(1);
dc(n,n)=rdamp(n-1);
dk(n,n)=gangdu(n-1);
for i=1:n-1;
   dc(i,i+1)=-rdamp(i);
   dc(i+1,i)=-rdamp(i);
   dk(i,i+1)=-gangdu(i);
   dk(i+1,i)=-gangdu(i);
end
for i=1:n
   dc(i,i)=dc(i,i)+adamp(i);
end

%*****质量刚度阻尼矩阵***********
M=dm;
C=dc;
K=dk;

%*****常数赋值***********
w=zhuansu*2*pi/60;% rad/s,弧度/秒 ，转速下的w
pi=3.1415926;

mj=73.07;%惯性 质量 ,kg

lamd=0.2544;%曲柄连杆比
stroke=0.29;   %m 冲程 
R=stroke/2;  %曲柄半径
D=0.28;     %气缸直径m

%*****气缸压力转化成曲柄转矩***********
load  PA6.dat%=========================================================================================change
alph=PA6(:,1);%曲柄转角%=======================================================================================change
  
pressure1=PA6(:,2);%气缸内气体压力,并且单位由原来的bar转换成pa%================================================change

%求解beta角度

alph=alph*pi/180;%，曲柄转角

 for i=1:nn%等同于时间
     
 if i>240
         ii=rem(i,240)+1;%rem(x,Y)表示if Y ~= 0, returns X - n.*Y where n = fix(X./Y)因为i小于240，n=0,ii=i,表明了以240为周期(3度一步
         %四冲程为两周720度，720/3=240度),每步时间步长为d，nn又为总时间与d之商，所以i=240代表转两周的时间
  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));%，为切向力矩阵,单位：N.M======套用公式为 教授 149页 6－3－2,6
  %Mg821(i)=mj*R*sin(alph(ii)+beta(ii))/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

%以上是周期力，只求一个周期就可

 %********************击振力***************只是为了求初加速度
if9=[0 60 600 660 120 180 480 540 240 300 360 420]/3;%发火间隔角，2度一个齿   顺序为：1-4-2-6-3-5 为四冲程双列内燃机
for i=1:12%除以3是转换成在那一时间步发火
  inc=if9(i);
  Mg9(i)=Mg8(inc+1);% 将单自由度转换成不同初相位的多自由度可见没有哪一个质量超过一个周期240
end

for i=1:6%循环控制量可以相同互不影响
    
  Mg(i)=Mg9(i*2-1)+Mg9(i*2);% [1 1 2 2 3 3 4 4 5 5 6 6]
end
for i=4:9%初始加速度 4-9为汽缸
   z(i)=Mg(i-3)/M(i,i);%Mg(1)表示4-4缸的合力
end

%***************************
ga=0.5;
bi=0.25;
a0=1.0/bi/d/d;
a1=ga/bi/d;
a2=1.0/bi/d;
a3=0.5/bi-1.0;
a4=ga/bi-1.0;
a5=d/2.0*(ga/bi-2.0);
a6=d*(1.0-ga);
a7=ga*d;
%形成等效刚度矩阵
ka=K+a0.*M+a1.*C;
%带入初值

x=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]';%初始位移
y=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]';%初始速度
z;%%初始加速度
  t=0;
 %***************开始进行计算******
%形成激励向量

for jj=1:10%整个用的大循环来迭代用以提高精度%可否用判断来代替
   for i=1:nn %每一部进行逐步积分计算

   ca=0;
   ma=0;
    f=0;         
   ca=C*(a1.*x+a4.*y+a5.*z);
   ma=M*(a0.*x+a2.*y+a3.*z);
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%将时间布与相位吻合
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对气缸爆发压力相位周期循环处理1
   for ii=1:12%ii示所应的气缸标号，从一到六
      ik=i+if9(ii);%if9(ii)示发火间隔角矩阵，加上i后表示各个气缸这个在该个积分步所对应的曲柄转角
      if(ik>240)
         ik=rem(ik,240)+1;%rem(x,Y)表示if Y ~= 0, returns X - n.*Y where n = fix(X./Y)，表示进入下一周期
      end
      Mg(ii)=Mg8(ik);%当有的质量相位超过720，将循环进入下一周期
   end
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对应气缸的合力1
   for ii=1:6
       f(ii+3)=Mg(ii*2-1)+Mg(ii*2);%对应汽缸对轴的合力
   end
   f(1)=0;
   f(2)=0;
   f(3)=0;
   for jj=10:14
       f(jj)=0;
   end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%电机上作用的冲击脉冲载荷的加载%%%%%%%%%%%%%%% 1
      if i<=Cnn1
       f(15)=42975*i*d;%表示时间步上的脉冲，程斜坡
       f(16)=42975*i*d;
   elseif i<=Cnn2
       f(15)=42975;%矩形脉冲
       f(16)=42975;
   elseif i<=Cnn3
       f(15)=42975*(6.5-i*d);%此与上两段脉冲成一梯形
       f(16)=42975*(6.5-i*d);
   else
      f(15)=0;
      f(16)=0;
  end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
fff(i,:)=f(:);
   f=f';  
   fff(i,:)=f(:);
   fa=f+ma+ca;  
   x1=ka\fa;
     z1=a0.*(x1-x)-a2.*y-a3.*z;
     y1=y+a6*z+a7*z1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     1
 
    if i>1%在 这里是表示在每一步计算的时候都把角度的值清零，但是还考虑到速度和加速度的影响 
      countx2=x;
     end
    
     x=x1;
     y=y1;
     z=z1;
     
     zhuanzhix=x';
     
          for ll=1:(n-1)
         xiangduix(ll)=zhuanzhix(ll)-zhuanzhix(ll+1);%相对扭转角
         xdnj(ll)=xiangduix(ll)*gangdu(ll);%弹性力
     end
    
     
     xiangduix=xiangduix';
     %xdnj=xdnj';
     
   ci=i-240;%周期 当i大于240时进入下一个周期
   if (ci<1)
      ci=1;
       else 
       ci=ci;
   end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
    xx(:,ci)=x;%xx是一个16*（nn-240）阵  第一个周期只保存当前值，以后每步计算值都要保存，有足够的值保证图形的精度
    yy(:,ci)=y;
    zz(:,ci)=z;
   
   
     XDNJ(:,ci)=xdnj';
    f4(:,ci)=f(4);% 
    f16(:,ci)=f(16);

     
   end 
    
    
end
%%%%%%%%%%%%%%%%%%%%%%%% 15个轴段应力  %%%%%%%%%%%%%%%%%%%%%%%%%%%
wn12=pi*(0.21)^3/16;
wn23=pi*(0.21)^3/16;
wn34=pi*(0.23)^3/16;
wn45=pi*(0.23)^3/16;%所以4到9为汽缸
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;%两个电枢之间

%ff4=ff(4,:);
%ff16=ff(16,:);


XDNJ12=XDNJ(1,:);
yingli12=XDNJ12/wn12/1e6;
maxyl12=max(yingli12);
minyl12=min(yingli12);
feng12=maxyl12-minyl12;%峰峰植

XDNJ23=XDNJ(2,:);
yingli23=XDNJ23/wn23/1e6;
maxyl23=max(yingli23);
minyl23=min(yingli23);
feng23=maxyl23-minyl23;

XDNJ34=XDNJ(3,:);
yingli34=XDNJ34/wn34/1e6;
maxyl34=max(yingli34);
minyl34=min(yingli34);
feng34=maxyl34-minyl34;

XDNJ45=XDNJ(4,:);
yingli45=XDNJ45/wn45/1e6;
maxyl45=max(yingli45);
minyl45=min(yingli45);
feng45=maxyl45-minyl45;

XDNJ56=XDNJ(5,:);
yingli56=XDNJ56/wn56/1e6;
maxyl56=max(yingli56);
minyl56=min(yingli56);
feng56=maxyl56-minyl56;

XDNJ67=XDNJ(6,:);
yingli67=XDNJ67/wn67/1e6;
maxyl67=max(yingli67);
minyl67=min(yingli67);
feng67=maxyl67-minyl67;

XDNJ78=XDNJ(7,:);
yingli78=XDNJ78/wn78/1e6;
maxyl78=max(yingli78);
minyl78=min(yingli78);
feng78=maxyl78-minyl78;

XDNJ89=XDNJ(8,:);
yingli89=XDNJ89/wn89/1e6;
maxyl89=max(yingli89);
minyl89=min(yingli89);
feng89=maxyl89-minyl89;

XDNJ910=XDNJ(9,:);
yingli910=XDNJ910/wn910/1e6;
maxyl910=max(yingli910);
minyl910=min(yingli910);
feng910=maxyl910-minyl910;

XDNJ1011=XDNJ(10,:);
yingli1011=XDNJ1011/wn1011/1e6;
maxyl1011=max(yingli1011);

XDNJ1112=XDNJ(11,:);
yingli1112=XDNJ1112/wn1112/1e6;
maxyl1112=max(yingli1112);


XDNJ1213=XDNJ(12,:);
yingli1213=XDNJ1213/wn1213/1e6;
maxyl1213=max(yingli1213);
minyl1213=min(yingli1213);
feng1213=maxyl1213-minyl1213;

XDNJ1314=XDNJ(13,:);
yingli1314=XDNJ1314/wn1314/1e6;
maxyl1314=max(yingli1314);
minyl1314=min(yingli1314);
feng1314=maxyl1314-minyl1314;

XDNJ1415=XDNJ(14,:);
yingli1415=XDNJ1415/wn1415/1e6;
maxyl1415=max(yingli1415)
minyl1415=min(yingli1415);
feng1415=maxyl1415-minyl1415;;

XDNJ1516=XDNJ(15,:);
yingli1516=XDNJ1516/wn1516/1e6;
maxyl1516=max(yingli1516);
minyl1516=min(yingli1516);
feng1516=maxyl1516-minyl1516;

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

i=1:(nn-240);%241到nn
%t=d*ij;


figure(1)
plot(d*i,yingli12,'g');title('第34轴段变化曲线');


figure(2)
plot(d*i,yingli23,'g');title('第34轴段变化曲线');


figure(3)
plot(d*i,yingli34,'g');title('第34轴段变化曲线');
figure(4)
plot(d*i,yingli45,'r');title('第45轴段变化曲线');
figure(5)
plot(d*i,yingli56,'g');title('第56轴段变化曲线');
figure(6)
plot(d*i,yingli67,'r');title('第67轴段变化曲线');
figure(7)
plot(d*i,yingli78,'g');title('第78轴段变化曲线');
figure(8)
plot(d*i,yingli89,'r');title('第89轴段变化曲线');
figure(9)
plot(d*i,yingli910,'r');title('第910轴段变化曲线');
figure(10)
plot(d*i,yingli1011,'r');title('第1011轴段变化曲线');
figure(11)
plot(d*i,yingli1112,'r');title('第1112轴段变化曲线');

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