xiaoci1112wanziyoumotai.m

内燃机转子仿真

%本程序是固有频率和振型计算（jiangnan厂推进情况）
%没有分支，日期：2005年9月15日 

close all
clear all
n=16;
n1=n-1;
ma=zeros(n,n);
c=zeros(n,n);
kk=zeros(n,n);


rate=[1 1 1 1 1   1 1 1 1 1   1 1 1 1 1  1 1 1 1 1   1 1 1 1 1  ];
%惯量

E=[19.627 9.599 0.509 13.018  9.908  13.018 13.018 9.908 13.018 42.11 49.09 177.6 14007.5 11.759585 308.155203 310.113323];
F=[0.651 39.86 32.75 32.75 32.75 32.75 32.75 32.75 325 0.316 205 9.7768 11.75483 3.11507 11.10479]*1.0e6;

E1=E;

for i=1:n;
   ma(i,i)=E(i)*rate(i)*rate(i);
end
for i=1:n-1
   F(i)=F(i)*rate(i)*rate(i);%THIS IS OK
end
%************
for i=1:n-1
	kk(i,i+1)=-F(i);
	kk(i+1,i)=-F(i);
end
%******************
for i=2:n-1
	kk(i,i)=F(i-1)+F(i);
end
kk(1,1)=F(1);
kk(n,n)=F(n-1);

%*************************计算固有频率和振型*****************************
%以下两种方法都可以-----方法一
 %j1=ma^0.5;
%j=inv(j1);
%方法二
j1=chol(ma);
 %[j1 j2]=lu(ma);
j=inv(j1);
b=j*kk*j';
[v,x]=eig(b);
for i=1:n
   x1(i)=sqrt(x(i,i))/2/3.1415926*1;
end

%********计算振型********
a=j'*v;%

 for i=1:n
         a(:,i)=a(:,i)/a(1,i);
 end
      
 %c1=a(:,23);
 %c2=a(:,21);
 %c3=a(:,20);
%c4=a(:,19);

%for i=1:n
 %  c1(i)=c1(i)*rate(i);%注意振型要乘速比
 %end
%for i=1:n
 %  c2(i)=c2(i)*rate(i);
 %end
%for i=1:n
 %  c3(i)=c3(i)*rate(i);
 %end
%for i=1:n
 %  c4(i)=c4(i)*rate(i);
 %end