wilson2.m

内燃机转子仿真

%function xy2=wilson2(M,C,K,dK,F,x1,x2,b,n,w);
clear all
clc
close all
M=[1 1 0;1 2 0;0 0 3];%3
C=[1 0 0;0 0 1;0 1 0];%3
K2=[2 1 0;1 3 0;0 0 4];%3
F=[2;4;6];%3
ci=2;
n=length(F);
t=0;
dK=[0 0 0;0 0.0 0;0 0 0];
K=K2+dK;
S=inv(M)*K;
[px pn]=eig(S);
pn=pn^0.5;
p=diag(pn);
fmin=min(p);
w=500;
if w==0
Tmax=1/2/pi/fmin;
else
Tmax=2*pi/w;
end  
h=Tmax/50;
tspan=10;
dn=rem(tspan,h);
dnn=tspan-dn;
n1=tspan/h;
f0=F*sin(w*t);
n=length(F);
i=1;
x1=[-1;-3;-2];%3
x2=[1;2;3];%3
x3=inv(M)*(f0-C*x2-K*x1);
xy=zeros(3*n+1,n1);
xy(1,i)=t;
xy(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=2:640
t=t+h;
%f(i)=jili(F,dtF);
dK=biangd(mid,D,de,ni,E,n)
K=K2+dK*sin(w*t);
K1=K+a0*M+a1*C;
xh=x1+h*(x2+x3*h);
xh1=x1;
while abs(x1-xh)>0.0001
fb=(F*sin(w*(t-h))-dK*sin(w*(t-h))*xh1)+b*(F*sin(w*t)-dK*sin(w*t)*xh-F*sin(w*(t-h))+dK*sin(w*(t-h))*xh1);
fa=fb+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;
xh=x1;
x2=xx2;
x3=xx3;
end
xy(1,i)=t;
xy(2:n+1,i)=x1;
xy(n+2:2*n+1,i)=x2;
xy(2*n+2:3*n+1,i)=x3;
end
figure(1);
plot(xy(1,:),xy(2,:));
save('biank','xy');