zbyz.m

内燃机转子仿真

clear all
clc
m=[2 0;0 1];
c=[0 0;0 0];
k=[6 -2;-2 4];
F=[0;10];
n=length(F);
x10=zeros(n,1);
x20=zeros(n,1);
f0=zeros(n,1);
x30=inv(m)*(f0-c*x20-k*x10); 
b=1.4;
dt=0.00281;
n1=10000;
i=1;
a1=6/b^3/dt^2;
a2=6/b^2/dt;
a3=1-3/b;
a4=dt/2;
a5=dt^2/6;
a6=6/b^2/dt^2;
a7=6/b/dt;
a8=3/b/dt;
a9=b*dt/2;
xyz=zeros(2*n+1,n1+1);
xyz(1,i)=i;
xyz(n:n+1,i)=x10;
xyz(n+2:2*n+1,i)=x20;

for i=2:n1
fa=F+m*(a6*x10+a7*x20+2*x30)+c*(a8*x10+2.*x20+a9*x30);
 x11=inv(k)*fa;
     x3=a1*(x11-x10)-a2*x20+a3*x30;
     x2=x20+a4*(x3+x30);
     x1=x10+dt*x20+(x3+2*x30)*a5;
     x10=x1;
     x20=x2;
     xyz(1,i)=i;
     xyz(n:n+1,i)=x10;
     xyz(n+2:2*n+1,i)=x20;
end