xierbaote.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=0.028;
%h=pi/400;
tspan=28;
st=rem(tspan,h);
tt=tspan-st;
nn=tt/h;
i=1;
hf=f*sin(w*(i-1)*h);
x1=[0;0];
x2=[0;0];
x3=inv(m)*(hf-k*x1-c*x2);
xierbote=zeros(n+1,nn);
xierbote(:,i)=[0;x1];
de=0.5;
oe=1.2;
maxbe=oe/2/(oe+1);
minbe=(2*oe^2-1)/4/(2*oe^3-1);
fprintf('%f,%f',minbe,maxbe);
be=input('the data be:');
%be=0.25;
a0=de*oe^2*h;
a1=be*oe^3*h^2;
a2=1-oe;
a3=(1-de*oe)*oe*h;
a4=oe*h;
a5=oe^2*h^2*(0.5-be*oe);
a6=(1-de)*h;
a7=de*h;
a8=(0.5-be)*h^2;
a9=be*h^2;
m1=oe*m+a0*c+a1*k;
for i=2:nn
    t=(i-1)*h;
    hf=a2*f+oe*f;
    r=hf-a2*m*x3-c*(x2+a3*x3)-k*(x1+a4*x2+a5*x3);
    y3=inv(m1)*r;
    y2=x2+a6*x3+a7*y3;
    y1=x1+h*x2+a8*x3+a9*y3;
    x1=y1;
    x2=y2;
    x3=y3;
    xierbote(:,i)=[t;x1];
end
save('xierbotedata','xierbote');
plot(xierbote(1,:),xierbote(2,:),xierbote(1,:),xierbote(3,:));