设为首页收藏本站联系我们
切换到宽版

声振论坛

 找回密码
 我要加入

QQ登录

只需一步，快速开始

声振论坛»论坛 › 基础理论 › 一般力学和振动理论 › 求大虾帮看一下求临界转速的程序，结果怎么都是NAN啊？
返回列表 发新帖
查看: 1270|回复: 3

[转子动力学] 求大虾帮看一下求临界转速的程序，结果怎么都是NAN啊？

[复制链接]
bingye0724
发表于 2008-9-25 09:36 | 显示全部楼层 |阅读模式

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要 登录 才可以下载或查看，没有账号？我要加入

x
本人用子结构传递矩阵法，用频率搜索法编制的程序。不知道怎么回事，结果全是NAN！费了1个多月时间找问题，还是没结果啊！急！不甚感激！！！:handshake
clc
clear
%计算wp7发动机的前二阶临界转速
M=[1.1074,40.6234,1.3623,0.97489,29.9160,0.74109,25.3895,0.93389,1.9347,3.024,...
3.5131,0.92958,8.4052,0.74066,1.3265,0.91047,61.4905,0,2.4238,21.2588,...
22.7067,1.2725,2.8241,2.8521,2.2117,9.6168,1.2497,1.468,4.0509,1.523,...
43.9880,0,3.8008,37.8808,0];

L=[0.06032,0.06032,0.045865,0.046035,0.055015,0.055485,0.09083,0.05107,0.0501,0.12771,...
0.27479,0.059,0.5111,0.0399,0.06,0.045,0.045,0,0.07309,0.07309,...
0.052212,0.019088,0.0815,0.0984,0.086173,0.43483,0.085,0.04,0.0705,0.0305,...
0.0305,0,0.0713,0.0713,0];

EI=[3.23E+03,3.23E+03,3.08E+05,4.06E+05,3.42E+05,3.42E+05,3.42E+05,4.44E+05,4.40E+05,3.24E+05,...
1.54E+05,2.67E+05,3.80E+05,4.13E+05,9.83E+05,1.96E+06,1.96E+06,1.96E+06,2.44E+07,2.44E+07,...
4.07E+07,6.59E+06,7.24E+05,7.04E+05,5.63E+05,7.53E+05,4.63E+05,1.25E+06,3.69E+05,3.47E+05,...
3.47E+06,3.47E+06,2.20E+07,2.20E+07,2.20E+07];

Jp=[0,0.3718,0,0,0.6944,0,0.6238,0,0,0,...
0,0,0,0,0,0,1.0078,0,0,0.4998,...
0.6476,0,0,0,0,0,0,0,0,0,...
0.7394,0,0,1.2561,0];

K1=9.8270E7;
K2=6.9128E8;
K3=1.8475E8;

K=[0,0,K1,0,0,0,0,0,0,0,...
    0,0,0,0,0,0,0,0,0,0,...
    0,0,0,K2,0,0,0,0,K3,0,...
    0,0,0,0,0];

PP=[0,0,0,0,0,0,0,0,1E15,0,...
    0,0,0,0,0,0,0,0,0,0,...
    0,0,0,0,0,0,0,0,0,0,...
    0,0,0,0,0];

U=[0,0,0,0,0,0,0,0,0,0,...
   0,0,0,0,0,0,0,0,0,0,...
   0,0,0,0,0,0,0,0,0,0,...
   0,0,0,0,0];
%求解高压转子激起的前二阶临界转速
XX=[1;0;0;0;0];SS=[0;1;0;0;0];MM=[0;0;1;0;0];QQ=[0;0;0;1;0]; %单位状态参数列阵

i=0;
for ww=1:250:2000;
     w=1;
     i=i+1;
    for j=1:2;
    s=0;
    q=0;
    while (s*q>=0)
T=TL7(9,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;M9X1=T(3);%M9X1
T=TL7(9,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M9S1=T(3);%M9S1

T=TL7(10,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X10X1=T(1);%X10X1
T=TL7(10,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X10S1=T(1);%X10S1
T=TL7(10,9,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X10P9=T(1);%X10P9

T=TH7(22,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X22X19=T(1);%X22X19
T=TH7(22,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X22S19=T(1);%X22S19

T=TH7(22,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;X22M33=T(1);%X22M33
T=TH7(22,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;X22Q33=T(1);%X22Q33

T=TL7(14,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X14X1=T(1);%X14X1
T=TL7(14,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X14S1=T(1);%X14S1
T=TL7(14,9,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X14P9=T(1);%X14P9

T=TH7(30,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X30X19=T(1);%X30X19
T=TH7(30,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X30S19=T(1);%X30S19

T=TH7(30,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;X30M33=T(1);%X30M33
T=TH7(30,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;X30Q33=T(1);%X30Q33
T=TL7(14,10,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;X14R10=T(1);%X14R10

T=TH7(21,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X21X19=T(1);%X21X19
T=TH7(21,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X21S19=T(1);%X21S19

T=TH7(21,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;S21X19=T(2);%S21X19
T=TH7(21,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;S21S19=T(2);%S21S19

T=TL7(18,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;M18X1=T(3);%M18X1
T=TL7(18,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M18S1=T(3);%M18S1
T=TL7(18,9,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M18P9=T(3);%M18P9

T=TL7(18,10,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M18R10=T(3);%M18R10
T=TL7(18,14,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M18R14=T(3);%M18R14

T=TL7(18,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;Q18X1=T(4);%Q18X1
T=TL7(18,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;Q18S1=T(4);%Q18S1
T=TL7(18,9,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;Q18P9=T(4);%Q18P9

T=TL7(18,10,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q18R10=T(4);%Q18R10
T=TL7(18,14,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q18R14=T(4);%Q18R14


T=TH7(32,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;M32X19=T(3);%M32X19
T=TH7(32,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M32S19=T(3);%M32S19

T=TH7(32,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;M32M33=T(3);%M32M33
T=TH7(32,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M32Q33=T(3);%M32Q33

T=TH7(32,22,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M32R10=T(3);%M32R10
T=TH7(32,30,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M32R14=T(3);%M32R14

T=TH7(32,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;Q32X19=T(4);%Q32X19
T=TH7(32,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;Q32S19=T(4);%Q32S19

T=TH7(32,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;Q32M33=T(4);%Q32M33
T=TH7(32,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q32Q33=T(4);%Q32Q33

T=TH7(32,22,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q32R10=T(4);%Q32R10
T=TH7(32,30,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q32R14=T(4);%Q32R14

T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;M35X33=T(3);%M35X33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M35S33=T(3);%M35S33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;M35M33=T(3);%M35M33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M35Q33=T(3);%M35Q33

T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;Q35X33=T(4);%Q35X33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;Q35S33=T(4);%Q35S33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;Q35M33=T(4);%Q35M33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q35Q33=T(4);%Q35Q33

HH=[M9X1,M9S1,0,0,0,0,0,0,0,0,0;
     X10X1,X10S1,X10P9,-X22X19,-X22S19,0,0,X22M33,X22Q33,0,0;
     X14X1,X14S1,X14P9,-X30X19,-X30S19,0,0,X30M33,X30Q33,X14R10,0;
     0,0,0,X21X19,X21S19,-1,0,0,0,0,0;
     0,0,0,S21X19,S21S19,0,-1,0,0,0,0;
     M18X1,M18S1,M18P9,0,0,0,0,0,0,M18R10,M18R14;
     Q18X1,Q18S1,Q18P9,0,0,0,0,0,0,Q18R10,Q18R14;
     0,0,0,M32X19,M32S19,0,0,-M32M33,-M32Q33,-M32R10,-M32R14;
     0,0,0,Q32X19,Q32S19,0,0,-Q32M33,-Q32Q33,-Q32R10,-Q32R14;
     0,0,0,0,0,M35X33,M35S33,M35M33,M35Q33,0,0;
     0,0,0,0,0,Q35X33,Q35S33,Q35M33,Q35Q33,0,0];


    s=det(HH);
    w=w+1;
    T=TL7(9,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;M9X1=T(3);%M9X1
T=TL7(9,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M9S1=T(3);%M9S1

T=TL7(10,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X10X1=T(1);%X10X1
T=TL7(10,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X10S1=T(1);%X10S1
T=TL7(10,9,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X10P9=T(1);%X10P9

T=TH7(22,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X22X19=T(1);%X22X19
T=TH7(22,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X22S19=T(1);%X22S19

T=TH7(22,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;X22M33=T(1);%X22M33
T=TH7(22,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;X22Q33=T(1);%X22Q33

T=TL7(14,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X14X1=T(1);%X14X1
T=TL7(14,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X14S1=T(1);%X14S1
T=TL7(14,9,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X14P9=T(1);%X14P9

T=TH7(30,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X30X19=T(1);%X30X19
T=TH7(30,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X30S19=T(1);%X30S19

T=TH7(30,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;X30M33=T(1);%X30M33
T=TH7(30,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;X30Q33=T(1);%X30Q33
T=TL7(14,10,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;X14R10=T(1);%X14R10

T=TH7(21,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;X21X19=T(1);%X21X19
T=TH7(21,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;X21S19=T(1);%X21S19

T=TH7(21,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;S21X19=T(2);%S21X19
T=TH7(21,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;S21S19=T(2);%S21S19

T=TL7(18,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;M18X1=T(3);%M18X1
T=TL7(18,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M18S1=T(3);%M18S1
T=TL7(18,9,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M18P9=T(3);%M18P9

T=TL7(18,10,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M18R10=T(3);%M18R10
T=TL7(18,14,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M18R14=T(3);%M18R14

T=TL7(18,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;Q18X1=T(4);%Q18X1
T=TL7(18,1,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;Q18S1=T(4);%Q18S1
T=TL7(18,9,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;Q18P9=T(4);%Q18P9

T=TL7(18,10,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q18R10=T(4);%Q18R10
T=TL7(18,14,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q18R14=T(4);%Q18R14


T=TH7(32,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;M32X19=T(3);%M32X19
T=TH7(32,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M32S19=T(3);%M32S19

T=TH7(32,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;M32M33=T(3);%M32M33
T=TH7(32,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M32Q33=T(3);%M32Q33

T=TH7(32,22,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M32R10=T(3);%M32R10
T=TH7(32,30,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M32R14=T(3);%M32R14

T=TH7(32,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;Q32X19=T(4);%Q32X19
T=TH7(32,19,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;Q32S19=T(4);%Q32S19

T=TH7(32,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;Q32M33=T(4);%Q32M33
T=TH7(32,21,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q32Q33=T(4);%Q32Q33

T=TH7(32,22,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q32R10=T(4);%Q32R10
T=TH7(32,30,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q32R14=T(4);%Q32R14

T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;M35X33=T(3);%M35X33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;M35S33=T(3);%M35S33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;M35M33=T(3);%M35M33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;M35Q33=T(3);%M35Q33

T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*XX;Q35X33=T(4);%Q35X33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*SS;Q35S33=T(4);%Q35S33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*MM;Q35M33=T(4);%Q35M33
T=TH7(35,33,M,L,K,EI,Jp,PP,w,ww,U);T=T*QQ;Q35Q33=T(4);%Q35Q33

HH=[M9X1,M9S1,0,0,0,0,0,0,0,0,0;
     X10X1,X10S1,X10P9,-X22X19,-X22S19,0,0,X22M33,X22Q33,0,0;
     X14X1,X14S1,X14P9,-X30X19,-X30S19,0,0,X30M33,X30Q33,X14R10,0;
     0,0,0,X21X19,X21S19,-1,0,0,0,0,0;
     0,0,0,S21X19,S21S19,0,-1,0,0,0,0;
     M18X1,M18S1,M18P9,0,0,0,0,0,0,M18R10,M18R14;
     Q18X1,Q18S1,Q18P9,0,0,0,0,0,0,Q18R10,Q18R14;
     0,0,0,M32X19,M32S19,0,0,-M32M33,-M32Q33,-M32R10,-M32R14;
     0,0,0,Q32X19,Q32S19,0,0,-Q32M33,-Q32Q33,-Q32R10,-Q32R14;
     0,0,0,0,0,M35X33,M35S33,M35M33,M35Q33,0,0;
     0,0,0,0,0,Q35X33,Q35S33,Q35M33,Q35Q33,0,0];

   q=det(HH);
   
    end
    wn(j)=(w*abs(s)+(w-1)*abs(q))/(abs(s)+abs(q));
    end
  a(i,:)=wn;
end
xx=1:250:2000;
xxx=0:2000;

plot(xx,a(:,1),xx,a(:,2),xxx,0.8319*xxx+229.5)

%传递矩阵
function THH=TH7(mm,nn,M,L,K,EI,Jp,PP,w,ww,U)
T=1;
for i=mm-1:-1:nn;
T=T*[1+L(i)^3*(M(i)*w^2-K(i))/(6*EI(i)),L(i)+L(i)^2*(Jp(i)*w^2/2)/(2*EI(i)),(L(i)^2/(2*EI(i)))+(L(i)/PP(i)),L(i)^3/(6*EI(i)),L(i)^3*U(i)*w^2/(6*EI(i));
    L(i)^2*(M(i)*w^2-K(i))/(2*EI(i)),1+L(i)*(Jp(i)*w^2/2)/(EI(i)),(L(i)/EI(i))+(L(i)/PP(i)),L(i)^2/(2*EI(i)),L(i)^2*U(i)*w^2/(2*EI(i));
L(i)*(M(i)*w^2-K(i)),Jp(i)*w^2/2,1,L(i),L(i)*U(i)*w^2;
M(i)*w^2-K(i),0,0,1,U(i)*w^2;
0,0,0,0,1];
end

THH=T;

function TLL=TL7(mm,nn,M,L,K,EI,Jp,PP,w,ww,U)
T=1;
for i=mm-1:-1:nn;
T=T*[1+L(i)^3*(M(i)*w^2-K(i))/(6*EI(i)),L(i)+L(i)^2*(((2*ww/w)-1)*Jp(i)*w^2/2)/(2*EI(i)),(L(i)^2/(2*EI(i)))+(L(i)/PP(i)),L(i)^3/(6*EI(i)),L(i)^3*U(i)*w^2/(6*EI(i));
    L(i)^2*(M(i)*w^2-K(i))/(2*EI(i)),1+L(i)*(((2*ww/w)-1)*Jp(i)*w^2/2)/(EI(i)),(L(i)/EI(i))+(L(i)/PP(i)),L(i)^2/(2*EI(i)),L(i)^2*U(i)*w^2/(2*EI(i));
L(i)*(M(i)*w^2-K(i)),((2*ww/w)-1)*Jp(i)*w^2/2,1,L(i),L(i)*U(i)*w^2;
M(i)*w^2-K(i),0,0,1,U(i)*w^2;
0,0,0,0,1];
end

TLL=T;
回复
分享到：

使用道具 举报

zhanghongxing
发表于 2008-9-25 11:35 | 显示全部楼层
bu dong le !
回复 支持 反对

使用道具 举报

ChaChing
发表于 2008-9-26 09:00 | 显示全部楼层

回复 楼主 bingye0724 的帖子

debug时找找那些地方有除以0, 就不难找出那儿出问题
应该自己练习找找, 日後功力自然提升

还有你的程式写法, 我直觉应该可以有空间精进
...
...
TH7 & TL7 function中, 矩阵中PP(i)很多为零,当然会出现inf及nan
回复 支持 反对

使用道具 举报

bingye0724
 楼主| 发表于 2008-10-1 14:57 | 显示全部楼层

回复 板凳 ChaChing 的帖子

谢谢啦！的却如此:handshake
回复 支持 反对

使用道具 举报

返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 13:51 , Processed in 0.062788 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表