niumake.m

内燃机转子仿真

clear all;
clc;
close all;
format long;
n=2;
%m=[1 1 0;1 2 0;0 0 3];
%c=[1 0 0;0 0 1;0 1 0];
%k=[2 1 0;1 3 0;0 0 4];
%f=[2;4;6];
m=[2 0;0 1];
c=[0 0;0 0];
k=[6 -2;-2 4];
f=[0;10];
w=1;
Bs=inv(m)*k;
[px pn]=eigs(Bs);
wf=pn^0.5;
%wf=2*pi*wf;
%for i=1:3
    %px(:,i)= px(:,i)/ px(1,i);
%end
p=diag(wf);
fmax=max(p);
T=1/fmax;
%h=T/1000;
h=0.028;
tspan=28;
st=rem(tspan,h);
tt=tspan-st;
nn=tt/h;
de=0.5;%de=3.5085955
da=0.25*(0.5+de)^2;
i=1;
hf=f*sin(w*(i-1)*h);
x1=[0;0];
x2=[0;0];
x3=inv(m)*(hf-k*x1-c*x2);
a0=1/da/h^2;
a1=de/da/h;
a2=1/da/h;
a3=(1/2/da)-1;
a4=(de/da)-1;
a5=0.5*h*(de/da-2);
a6=h*(1-de);
a7=de*h;
numak=zeros(n+1,nn);
numal(:,1)=[0;x1];

for i=2:nn
    t=(i-1)*h;
    k1=k+a0*m+a1*c;
    hf=f;
    r=hf+m*(a0*x1+a2*x2+a3*x3)+c*(a1*x1+a4*x2+a5*x3);
    y1=inv(k1)*r;
    y3=a0*(y1-x1)-a2*x2-a3*x3;
    y2=x2+a6*x3+a7*y3;
    x1=y1;
    x2=y2;
    x3=y3;
    numak(:,i)=[t;y1];
end
save('sjjl1','numak');
plot(numak(1,:),numak(2,:),numak(1,:),numak(3,:));