pinyujisuan.m

内燃机转子仿真

%在400——1100rpm转速下,从0.5－10谐次的 振动幅值，扭矩，应力 及合成值
close all
clear all
n=10;
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];        %inertia(kg*m**2)
F=[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]; %stiffness(n*m/rad)
ABS_DAMP=[0 0 230 230 230 230 230 230 0 0];    %绝对阻尼 N*M/RAD/S
REIA_DAMP=[2700 0 0 0 0 0 0 0 0 0];            %相对阻尼  N*M/RAD/S
DIAM=[0.1 210 210 210 210 210 210 300 220];   %mm   ！！第一个直径赋值0.1是为了在求应力时分母不为零
lamd=0.2544;%曲柄连杆比
stroke=29;   %cm
r=stroke/2;  %曲柄半径
d=28;     %气缸直径cm
speed=[400:50:1100]*2*pi/60;% rad/s
jiaodu=[360 600 120 480 240 0]*pi/180;  %转换为弧度，气缸的顺序为简化模型编号顺序的倒序
%*********************转动惯量矩阵
for i=1:n
    jj(i,i)=E(i);
end
%*********************刚度系数矩阵
for i=2:n
    kk(i,i)=F(i-1)+F(i);                  
end
for i=1:n-1
    kk(i,i+1)=-F(i);
    kk(i+1,i)=-F(i);
end
kk(1,1)=F(1);
%***********************阻尼矩阵
%绝对阻尼
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;  %总阻尼
%********************气体压力产生的切向干扰力矩Mt(g)
for i=1:15
    MEP(i)=20.153*(speed(i)*60/2/pi/1000)^0;      %在不同转速下的平均有效压力，注意单位。  发电机组特性，乘方为0；船舶为螺旋桨特性，乘方为2   
end 
load('qPA6.dat');
%load matlabpa6
for m=1:15                        %共划分为15个转速值
    for i=1:20            
        re(i)=(qPA6(2*i-1,1)*MEP(m)^3+qPA6(2*i-1,2)*MEP(m)^2+qPA6(2*i-1,3)*MEP(m)+qPA6(2*i-1,4))*1e-4;   %气缸压力切向力cos部分系数，即实数系数
        im(i)=(qPA6(2*i,1)*MEP(m)^3+qPA6(2*i,2)*MEP(m)^2+qPA6(2*i,3)*MEP(m)+qPA6(2*i,4))*1e-4;       %气缸压力切向力sin部分系数，即虚数系数
    end
    Ptg=re-j*im;   %气缸压力切向力kg/cm**2  注意符号，化为复数形式以后，虚部为负  见论文NO.40
    Pti(2)=73.5*4/pi/d^2*r*(speed(m))^2*(0.25*lamd);   %i*2－系数*2 曲柄连杆机构往复惯性力所产生的干扰力矩 只有第1，2，3，4谐次下赋值，其余为零
    Pti(4)=-73.5*4/pi/d^2*r*(speed(m))^2*(0.5);
    Pti(6)=-73.5*4/pi/d^2*r*(speed(m))^2*(0.75*lamd);
    Pti(8)=-73.5*4/pi/d^2*r*(speed(m))^2*(0.25*lamd^2);
    Pti(1)=0;
    Pti(3)=0;
    Pti(5)=0;
    Pti(7)=0;
    Pti(9:20)=0;
    for i=0.5:0.5:10        %i为谐次，k为质量号  给定简谐次数
        for k=3:8           %气体压力激振扭矩只在气缸处的元素为非零
            Mgg(k,i*2)=pi/4*d^2*r*Ptg(i*2)*1e-2*9.8*(cos(i*jiaodu(k-2))-j*sin(i*jiaodu(k-2)));%气缸激励力矩  ！！！注意 *9.8
            Mgi(k,i*2)=pi/4*d^2*r*(-j*Pti(i*2)*1e-4)*(cos(i*jiaodu(k-2))-j*sin(i*jiaodu(k-2)));%往复惯性力矩  %'-'表示系数相减，见论文NO.40  
            Mg(k,i*2)=Mgg(k,i*2)+Mgi(k,i*2); %总激振力矩  
        end
        Mc(:,i*2)=real(Mg(:,i*2));  %总激振力矩化为实部和虚部 见论文NO.42
        Ms(:,i*2)=imag(Mg(:,i*2));
        FF=kk-(speed(m)*i)^2*jj;
        GG=speed(m)*i*c;
        %**********求解
        A=[FF -GG;GG FF];
        A1=inv(A);
        M1(:,i*2)=[Mc(:,i*2);Ms(:,i*2)];
        y(:,2*i)=A1*M1(:,i*2);
        y1(:,2*i)=y(1:10,2*i);      %α
        y2(:,2*i)=y(11:20,2*i);     %β
        AA(:,i*2)=(y1(:,i*2).^2+y2(:,i*2).^2).^0.5;  %在不同谐次下的振幅A
        for t=1:9
            MM(t,i*2)=F(t)*((y1(t,2*i)-y1(t+1,2*i)).^2+(y2(t,2*i)-y2(t+1,2*i)).^2).^0.5;  %不同谐次下 轴段振动力矩M  单位 N*m
        end
    end
%************************求解时域NO.2 质量的振动位移幅值合成值，列为转速400-1100，行为转角10度-720度，转换为弧度
    for t=1:72         %时域角度从10度－720度  该段程序中涉及的公式见论文NO.43
    for k=1:20
        if y1(2,k)>=0
            v(k)=atan(y2(2,k)/y1(2,k));
        else
            v(k)=pi+atan(y2(2,k)/y1(2,k));
        end
        if k==1
            AAA(t,k)=AA(2,k).*cos(0.5*t*10*pi/180+v(k));
        else
            AAA(t,k)=AAA(t,k-1)+AA(2,k).*cos(0.5*k*t*10*pi/180+v(k));          %论文公式中的W0*t为曲轴转角，即时域中的10度－720度
        end
    end     %for k=1:20 循环语句的结束语
    end     %for t=1:72 循环语句的结束语
    AAA2(:,m)=AAA(:,20);                      %！！！NO.2的角位移振幅合成值
%************************求解时域各轴段振动力矩幅值各谐次的叠加，行为转角10度-720度，转换为弧度
    for q=1:9          %q为轴段编号
        for t=1:72         %时域角度从10度－720度  该段程序中涉及的公式见论文NO.46
            for k=1:20
                if (y1(q,k)-y1(q+1,k))>=0
                    u(k)=atan((y2(q,k)-y2(q+1,k))/(y1(q,k)-y1(q+1,k)));
                else 
                    u(k)=pi+atan((y2(q,k)-y2(q+1,k))/(y1(q,k)-y1(q+1,k)));
                end
                if k==1
                    MMM(t,k)=MM(q,k).*cos(0.5*t*10*pi/180+u(k));
                else
                    MMM(t,k)=MMM(t,k-1)+MM(q,k).*cos(0.5*k*t*10*pi/180+u(k));%0.5:0.5:10谐次合力矩
                end
            end     %for k=1:20 循环语句的结束语
        end     %for t=1:72 循环语句的结束语
        QMMM(:,q)=MMM(:,20);                     %各谐次按时域叠加后的值，列为轴段编号，行为曲轴转角度数10度－720度  
    end
%***************轴段振动力矩幅值的合成值 每一转速下的赋值 ：最大最小值的绝对值的算术平均值！！！
    for k=1:9
        ss(:,k)=QMMM(:,k);
        Mfuzhi(k,m)=(max(ss(:,k))-min(ss(:,k)))/2;      %！！！各轴段振动力矩的合成值，行为轴段编号，列为转速400－1100
    end
%**********************曲轴最大应力，按时域计算
    for t=2:7
        W(t)=pi*DIAM(t)^3/16;                  %曲轴各轴段的抗扭界面模量 mm**3
        TTT(t,m)=Mfuzhi(t,m)./W(t)*1e3;      %曲轴各轴段在各个转速下的应力  单位 N/mm**2  行为轴段编号，列为各个转速
    end
%**********************频域值 各谐次下的振幅 力矩及应力
    AA1(m,:)=AA(1,:);         %列为谐次20列，行为转速：No.M2角位移振幅，从400转开始，每隔50转，至1100转 
    AA2(m,:)=AA(2,:);
    MM7(m,:)=MM(7,:);         %列为谐次20列，行为转速：No.7轴段振动力矩，从400转开始，每隔50转，至1100转
    MM9(m,:)=MM(9,:);         %列为谐次20列，行为转速：No.9轴段振动力矩，从400转开始，每隔50转，至1100转
end              %！！！！不同转速值for m=1:15 循环语句的结束语
%***************NO.2 质量的振动位移幅值合成值 各转速下的赋值 ：最大最小值的绝对值的算术平均值！！！
for k=1:15
    s(:,k)=AAA2(:,k);
    fuzhi(k)=(max(s(:,k))-min(s(:,k)))/2;
end
%***************曲轴最大应力
for p=1:15                                %承受最大应力的轴段在各个转速下的应力值  单位 N/mm**2  
    Tfuzhi(p)=max(TTT(:,p));
end
%***************
x=[400:50:1100];
figure(1);
  plot(x,AA1(:,3),x,AA1(:,6),x,AA1(:,7),x,AA1(:,9),x,AA1(:,12),x,fuzhi(:),'.-')    %No.2角位移振幅 1.5谐次 3谐次 3.5谐次 4.5谐次 6谐次 及合成值
  axis([400 1100 0 8*1e-3])
figure(2);
  plot(x,MM9(:,3),x,MM9(:,6),x,MM9(:,9),x,MM9(:,12),x,Mfuzhi(7,:),'.-')             %No.7轴段振动力矩 1.5谐次 3谐次 4.5谐次 6谐次 及合成值
  axis([400 1100 0 75*1e3])
figure(3);
  plot(x,MM9(:,6),x,MM9(:,9),x,MM9(:,12),x,Mfuzhi(9,:),'.-')                        %No.9轴段振动力矩 3谐次 4.5谐次 6谐次
  axis([400 1100 0 55*1e3])
figure(4);
  plot(x,Tfuzhi(:))                                                            %曲轴最大应力
  axis([400 1100 0 70])
  figure(5);
  plot(x,AA2(:,3),x,AA2(:,6),x,AA2(:,7),x,AA2(:,9),x,AA2(:,12),x,fuzhi(:),'.-')    %No.2角位移振幅 1.5谐次 3谐次 3.5谐次 4.5谐次 6谐次 及合成值
  axis([400 1100 0 8*1e-3])