wilson11.m

内燃机转子仿真

clear all
clc
close all
n=3;
m=[1 1 0;1 2 0;0 0 3];
c=[1 0 0;0 1 0;0 0 1];
k=[2 1 0;1 3 0;0 0 4];
F=[2;4;6];
%[M,C,K]=arryf(m,c,k,n);
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/100;
t=0;
%h=pi/40;
dn=rem(T,h);
dnn=T-dn;
n1=T/h;
%f0=F*sin(t);
f0=F*sin(t);
i=1;
x1=[-1 -3 -2];%3
x1=x1';
x2=diag(eye(n))*0.02;%3
x3=inv(m)*(f0-c*x2-k*x1);
speed=60;
w=speed/60;
v=zeros(n+1,55000);
v(1,i)=t;
v(2:n+1,i)=x1;
b=1.4;
a0=6/b^2/h^2;
a1=3/b/h;
a2=2*a1;
a3=b*h/2;
a4=a0/b;
a5=-a2/b;
a6=1-3/b;
a7=h/2;
a8=h^2/6;
for i=1:55000
    K1=k+a0*m+a1*c;
    t=t+h;
fa=F*sin(t-h)+b*(F*sin(t)-F*sin(t-h))+m*(a0*x1+a2*x2+2*x3)+c*(a1*x1+2*x2+a3*x3);
xb=inv(K1)*fa;
xx3=a4*(xb-x1)+a5*x2+a6*x3;
xx2=x2+a7*(xx3+x3);
xx1=x1+h*x2+a8*(xx3+2*x3);
x1=xx1;
x2=xx2;
x3=xx3;
v(1,i)=t;
v(2:n+1,i)=x1;
end
save('sjjldata','v');
figure(1);
plot(v(1,:),v(2,:),v(1,:),v(3,:),'g',v(1,:),v(4,:),'r');