wilson12.m

内燃机转子仿真

%function xy=wilson1((M,C,K,dK,F,x1,x2,b,n);
clear all
clc
close all
M=[1.000 1.000 0.000;1.000 2.000 0.000;0.000 0.000 3.000];%3
C=[1.000 0.000 0.000;0.000 0.000 1.000;0.000 1.000 0.000];%3
K=[2.000 1.000 0.000;1.000 3.000 0.000;0.000 0.000 4.000];%3
F=[2.000;4.000;6.000];%3
t=0;
b=1.4;
S=inv(M)*K;
[px pn]=eig(S);
pn=pn^0.5;
p=diag(pn);
w=1;
fmin=min(p);
if w==0
Tmax=1/2/pi/fmin;
else
Tmax=2*pi/w;
end    
h=Tmax/320.0000;
dn=rem(Tmax,h);
dnn=Tmax-dn;
n1=Tmax/h;
f0=F*sin(t);
n=length(F);
x1=[-1;-3;-2];%3
x2=[1;2;3];%3
i=1;
x3=inv(M)*(f0-C*x2-K*x1);
xy1=zeros(n+1,5000);
xy1(1,i)=t;
xy1(2:n+1,i)=x1;
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;
 K1=K+a0*M+a1*C;
for i=2:5000
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;
xy1(1,i)=t;
xy1(2:n+1,i)=x1;
end
save('dingk','xy1');
figure(1)
plot(xy1(1,:),xy1(2,:),xy1(1,:),xy1(3,:),xy1(1,:),xy1(4,:));