zhongxinchafen1.m

内燃机转子仿真

clear all;
clc;
close all;
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];
%[M K C]=arryf(m,k,c,n);
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);
x10=[0;0];
x20=[0;0];
x30=inv(m)*(hf-k*x10-c*x20);
a0=1/h^2;
a1=1/2/h;
a2=2*a0;
a3=1/a2;
x0=x10-h*x20+a3*x30;
m1=a0*m+a1*c;
x=zeros(n+1,nn);
x(:,1)=[0;x10];
for i=2:nn
    t=i*h;
    hf=f;
r=hf-(k-a2*m)*x10-(a0*m-a1*c)*x0;
x11=inv(m1)*r;
x(:,i)=[t;x11];
x0=x10;
x10=x11;
end
save('zhongxinchafen','x');
plot(x(1,:),x(2,:),x(1,:),x(3,:));