求助：传递矩阵法求解固有频率和振型

ME!

%求解转子系统前三个临界转速和主振型的传递矩阵法

a=1;%截面剪切系数
u=0.3;%泊松比
rou=7800;
E=2.0e11;
L1=0.187;%轴的总长
G=E/(2*(1+u));
K1=2.0e7;%弹性支撑刚度
%
K=[K1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K1  ]; %弹性支撑的刚度
L=[16.5 13.5   4  10 5  6.5   3   6.5    15 12 20  5  20.5  5  10   5  6.5  13    10  ]*1e-3;%轴的分段
d1=[15   16    17 17 17 20.5  29  21.7   21 21 21  20  20.5 17 17   17  16  12.5  12.5 ]*1e-3; %%%光轴尺寸
I=pi*d1.^4/64;
A=pi*d1.^2/4;
v=6*E*I./(a*G*A.*L.^2);
M=zeros(1,length(L));
Jp=zeros(1,length(L));
Jd=zeros(1,length(L));
for jk=1:length(L)
m(jk)=rou*pi*d1(jk)^2*L(jk)/8 ;%轴的质量分布
jp(jk)=rou*A(jk)*L(jk)^2*d1(jk)^2/16;     %轴的极转动惯量Jp=mi*r^2/2;
end
M(:,2:length(L)-1)=m(:,2:length(L)-1)+m(:,3:length(L));
M(1)=m(1);
M(length(L))=m(length(L));
Jp(:,2:length(L)-1)=m(:,2:length(L)-1)+m(:,3:length(L));
Jp(1)=jp(1);
Jp(length(L))=jp(length(L));
Jd=Jp./2;
J=Jp-Jd;%直径转动惯量

k=0;
for w=0:10:20000;
for i=1:length(L)
T(:,:,i)=[1+(L(i)^3)*(1-v(i))*(M(i)*w^2-K(i))/(6*E*I(i))    L(i)+L(i)^2*J(i)*w^2/(2*E*I(i))    L(i)^2/(2*E*I(i))    L(i)^3*(1-v(i))/(6*E*I(i));
   (L(i)^2)*(M(i)*w^2-K(i))/(2*E*I(i))     1+L(i)*J(i)*w^2/(E*I(i))      L(i)/(E*I(i))     L(i)^2/(2*E*I(i));
   L(i)*(M(i)*w^2-K(i))    J(i)*w^2    1    L(i);
   M(i)*w^2-K(i) 0 0 1];   %传递矩阵
end
H=T(:,:,1);
for i2=2:length(L);
       H=T(:,:,i2)*H;
end
    F=H(3,1)*H(4,2)-H(3,2)*H(4,1) ; %剩余量
    if F*(-1)^k < 0 %求解临界转速
       k=k+1;
       wi(k)=w ;    %固有圆频率
       w=wi(k);
        ni(k)=wi(k)*30/pi ; %临界转速

end
end
ni=ni'  
wi=wi'
f=ni/60

%绘制前三阶模态振型
for i1=1:3;
    w=wi(i1);
      for j=1:length(L)-1;
T(:,:,j)=[1+(L(j)^3)*(1-v(j))*(M(j)*w^2-K(j))/(6*E*I(j)) L(j)+L(j)^2*J(j)*w^2/(2*E*I(j)) L(j)^2/(2*E*I(j)) L(j)^3*(1-v(j))/(6*E*I(j));
   (L(j)^2)*(M(j)*w^2-K(j))/(2*E*I(j)) 1+L(j)*J(j)*w^2/(E*I(j)) L(j)/(E*I(j)) L(j)^2/(2*E*I(j));
   L(j)*(M(j)*w^2-K(j)) J(j)*w^2 1 L(j);
   M(j)*w^2-K(j) 0 0 1];
end
H=T(:,:,1);
for j=2:length(L);
       H=T(:,:,j)*H;
   end
    b=-H(4,1)/H(4,2);
  X(:,1)=([1 b 0 0]');
for n=2:length(L)+1;
    X(:,n)=T(:,:,n-1)*X(:,n-1) ; %相邻两质点右边的传递关系
end
    x=zeros(1,length(L)+1);
for j1=1:length(L);
     y(j1)=X(1,j1);     %挠度   
     z(j1)=X(3,j1);   %弯矩
   x(j1+1)=L(j1)+x(j1);
end


y(length(L)+1)=X(1,length(L)+1);
z(length(L)+1)=X(3,length(L)+1);
y=y/max(abs(y)) ;  %归一化
z=z/max(abs(z))  ;  %归一化
subplot(3,1,i1)
plot(x,y,'b-',x,z,'r:');
Tit=['第一阶频率的振型和弯矩图';'第二阶频率的振型和弯矩图';'第三阶频率的振型和弯矩图'];
title(Tit(i1,:))   
xlabel('轴长'),ylabel('不平衡值')
%axis([0,L1,-1,1])
grid on
z;
end
legend('振型','弯矩')

振型图怎么是这个德行啊！
程序是按照论坛里面的一个例子做的，但是我没有加圆盘，就是一根光轴，计算的固有频率和ansys计算的、有限元法计算的结果相差都很大



补充内容 (2013-7-10 21:42):
转帖！

dw04116

振型图怎么能是“不平衡度”的纵坐标呢。
我没坐过转子方面的，但是如果是弯曲，纵坐标肯定是节点挠度，或者位移加速度什么的。
如果是扭转，那么肯定是扭转角度，速度之类……
仅供参考

ME!

dw04116 发表于 2013-3-25 11:14

                               
登录/注册后可看大图

振型图怎么能是“不平衡度”的纵坐标呢。
我没坐过转子方面的，但是如果是弯曲，纵坐标肯定是节点挠度，或 ...

状态向量的第一行就是挠度啊！

dw04116

我定义的计算，第一行就是第一个节点的广义位，第二行是第二个点的……
第一列是一阶频率……

jackyyjy

   b=-H(4,1)/H(4,2);
  X(:,1)=([1 b 0 0]');
for n=2:length(L)+1;
    X(:,n)=T(:,:,n-1)*X(:,n-1) ; %相邻两质点右边的传递关系
end

为什么   b=-H(4,1)/H(4,2);

Charles_Zh

你这M和JP的计算公式错了。M的个数比L的多一个。

补充内容 (2013-9-7 14:43):
Jp(:,2:length(L)-1)=m(:,2:length(L)-1)+m(:,3:length(L));计算极转动惯量的公式也错了。