wilson3.m

内燃机转子仿真

clear all
clc
close all
m=[5.98 1.08 1.04 2.913 2.913 2.913 2.913 2.913 51.463 0.6 1.115 3.944];
k=[0.0001 8.2 392.2 150 112.78 112.78 112.78 112.78 169.66 0.5 0.5 50.29 0.0001]*1e2;
c=[0 0 0 0 0 0 0 0 0 0 0 0 0];
F=zeros(12,1);
n=length(F);
[M,C,K]=arryf(m,c,k,n);
Bs=inv(M)*K;
[px pn]=eig(Bs);
wf=pn^0.5;
wf=2*pi*wf;
for i=1:12
    px(:,i)= px(:,i)/ px(1,i);
end
p=diag(wf);
fmax=min(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=[0.1 -0.2 0.3 0.2 -0.14 0.21 -0.13 0.31 0.05 -0.12 0.22 0.15];%3
x1=x1';
x2=diag(eye(n))*0.02;%3
x3=inv(M)*(f0-C*x2-K*x1);
speed=500;
w=speed/60;
v=zeros(n+1,100000);
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:100000
    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,:),'y');