wilson5.m

内燃机转子仿真

%function xy2=wilson2(M,C,K,dK,F,x1,x2,b,n,w);
clear all
clc
close all
format long;
n=1;
%M=[1 1 0;1 2 0;0 0 3];%3
%C=[1 0 0;0 0 1;0 1 0];%3
%K=[2 1 0;1 3 0;0 0 4];%3
%F=[2;4;6];%3
M=[2 0;0 1];
C=[0 0;0 0];
K=[6 -2;-2 4];
F=[0;10];
%M=2;
%K=20;
%C=0.5;
%F=1;
w=1;
S=inv(M)*K;
[px pn]=eig(S);
pn=pn^0.5;
p=diag(pn);
fmin=max(p);
Tmax=1/2/pi/fmin;
h=0.028;
tspan=28;
dn=rem(tspan,h);
dnn=tspan-dn;
n1=dnn/h;
t=0;
f0=F*sin(w*t);
n=length(F);
i=1;
x1=[0;0];%3
x2=[0;0];%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;
%xy(n+2:2*n+1,i)=x2;
%xy(2*n+2:3*n+1,i)=x3;
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;
 K1=K+a0*M+a1*C;
for i=2:n1
t=t+h;
%fb=F*sin(w*(t-h))+b*(F*sin(w*(t))-F*sin(w*(t-h)));
%ff(:,i)=F*sin(w*t);
fa=F+M*(a0*x1+a2*x2+2*x3)+C*(a1*x1+2*x2+a3*x3);
xb=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;
xy(1,i)=t;
xy(2:n+1,i)=xx1;
xy(n+2:2*n+1,i)=xx2;
xy(2*n+2:3*n+1,i)=xx3;
end
figure(1);
plot(xy(1,:),xy(2,:));
save('duibi','xy');