zhoujunbo 发表于 2010-1-13 21:46

求助——钟一谔转子动力学有限元法


当轴承为各向同性的时候,临界转速可以从频率方程求得。但是在Matlab中编程当节点数超过6个时就求不出符号行列式,得不到方程更无从求解临界转速。
遇到这种情况改怎么处理。

zhoujunbo 发表于 2010-1-13 21:53

主程序:
%%
clc
clear all;
Nshaft=8;                                                                                 %轴段数量;
%%
RotorE = 2.095e11;                                                                        %转子弹性模量;
RotorM = 7.85e3;                                                                        %转子材料密度
ShaftL = ;                                       %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1;                                                             %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0;                                                             %各轴段内径;
%%
LocationF=;                                                                        %支承所在节点编号;
SPstiffness=;                                                                     %支承刚度;

%%
U=pi*(ShaftDI.^2-ShaftDO.^2)*RotorM;                                                      %单位长度质量
a=;
A=a.';
%%                                                                           质量矩阵的计算
for i=1:Nshaft                                                                                 
    u=A(i,1);
    l=A(i,2);
    r=A(i,3);
    M(:,:,i)=Ms(u,l,r);
end;
                     
B=zeros(Nshaft*2+2,Nshaft*2+2);
%整体质量矩阵第1、2行
B(1:2,1:2)=B(1:2,1:2)+M(1:2,1:2,1);
B(1:2,3:4)=B(1:2,3:4)+M(1:2,3:4,1);
%整体质量矩阵第3至Nshaft*2行
for i=3:2:Nshaft*2
    j=i-2;
    B(i:(i+1),j:(j+1))=B(i:i+1,j:j+1)+M(3:4,1:2,(i-1)/2);
    j=i;
    B(i:i+1,j:j+1)=B(i:i+1,j:j+1)+M(3:4,3:4,(i-1)/2)+M(1:2,1:2,(i+1)/2);
    j=i+2;
    B(i:i+1,j:j+1)=B(i:i+1,j:j+1)+M(1:2,3:4,(i+1)/2);
end
%整体质量矩阵第Nshaft*2+1、Nshaft*2+2行
B(Nshaft*2+1:Nshaft*2+2,Nshaft*2-1:Nshaft*2)=B(Nshaft*2+1:Nshaft*2+2,Nshaft*2-1:Nshaft*2)+M(3:4,1:2,Nshaft);
B(Nshaft*2+1:Nshaft*2+2,Nshaft*2+1:Nshaft*2+2)=B(Nshaft*2+1:Nshaft*2+2,Nshaft*2+1:Nshaft*2+2)+M(3:4,3:4,Nshaft);
%%                                                                  回转矩阵的计算
for i=1:Nshaft
u=A(i,1);
l=A(i,2);
r=A(i,3);
J(:,:,i)=Js(u,l,r);
end;
C=zeros(Nshaft*2+2,Nshaft*2+2);
%整体回转矩阵第1、2行
C(1:2,1:2)=C(1:2,1:2)+J(1:2,1:2,1);
C(1:2,3:4)=C(1:2,3:4)+J(1:2,3:4,1);
%整体回转矩阵第3-Nshaft*2行
for i=3:2:Nshaft*2                                                                                          %节点号为(i+1)/2
    j=i-2;
    C(i:(i+1),j:(j+1))=C(i:i+1,j:j+1)+J(3:4,1:2,(i-1)/2);
    j=i;
    C(i:i+1,j:j+1)=C(i:i+1,j:j+1)+J(3:4,3:4,(i-1)/2)+J(1:2,1:2,(i+1)/2);
    j=i+2;
    C(i:i+1,j:j+1)=C(i:i+1,j:j+1)+J(1:2,3:4,(i+1)/2);
end
%整体回转矩阵第Nshaft*2+1、Nshaft*2+2行
C(Nshaft*2+1:Nshaft*2+2,Nshaft*2-1:Nshaft*2)=C(Nshaft*2+1:Nshaft*2+2,Nshaft*2-1:Nshaft*2)+J(3:4,1:2,Nshaft);
C(Nshaft*2+1:Nshaft*2+2,Nshaft*2+1:Nshaft*2+2)=C(Nshaft*2+1:Nshaft*2+2,Nshaft*2+1:Nshaft*2+2)+J(3:4,3:4,Nshaft);
%%                                                                      整体刚度矩阵的计算   
for i=1:Nshaft
    l=A(i,2);
    r=A(i,3);
    K(:,:,i)=Ks(l,r,RotorE);
end;
D=zeros(Nshaft*2+2,Nshaft*2+2);
%整体刚度矩阵第1、2行
D(1:2,1:2)=D(1:2,1:2)+K(1:2,1:2,1);
D(1:2,3:4)=D(1:2,3:4)+K(1:2,3:4,1);
%整体刚度矩阵第3至Nshaft*2行
m=1;
Size=size(LocationF);
for i=3:2:Nshaft*2
    j=i-2;
    D(i:(i+1),j:(j+1))=D(i:i+1,j:j+1)+K(3:4,1:2,(i-1)/2);
    j=i;
    if Size(1,2)==0;                                                                         %判断支承位置,添加支承刚度
      for i=3:2:Nshaft*2
            j=i-2;
            D(i:(i+1),j:(j+1))=D(i:i+1,j:j+1)+K(3:4,1:2,(i-1)/2);
            j=i;
            D(i:i+1,j:j+1)=D(i:i+1,j:j+1)+K(3:4,3:4,(i-1)/2)+K(1:2,1:2,(i+1)/2);
            j=i+2;
            D(i:i+1,j:j+1)=D(i:i+1,j:j+1)+K(1:2,3:4,(i+1)/2);
      end
    else
      if i==2*LocationF(m)-1;
            D(i:i+1,j:j+1)=D(i:i+1,j:j+1)+K(3:4,3:4,(i-1)/2)+K(1:2,1:2,(i+1)/2);
            D(i,j)=D(i,j)+SPstiffness(m);
            while m<Size(1,2)
                m=m+1;
            end
      else
            D(i:i+1,j:j+1)=D(i:i+1,j:j+1)+K(3:4,3:4,(i-1)/2)+K(1:2,1:2,(i+1)/2);
      end
    end
    j=i+2;
    D(i:i+1,j:j+1)=D(i:i+1,j:j+1)+K(1:2,3:4,(i+1)/2);
end
%整体刚度矩阵第Nshaft*2+1、Nshaft*2+2行
D(Nshaft*2+1:Nshaft*2+2,Nshaft*2-1:Nshaft*2)=D(Nshaft*2+1:Nshaft*2+2,Nshaft*2-1:Nshaft*2)+K(3:4,1:2,Nshaft);
D(Nshaft*2+1:Nshaft*2+2,Nshaft*2+1:Nshaft*2+2)=D(Nshaft*2+1:Nshaft*2+2,Nshaft*2+1:Nshaft*2+2)+K(3:4,3:4,Nshaft);
%%
syms w;
pin=-B*w^2+C*w^2+D;
Speed=det(pin);
digits(5);
ssss=vpa(solve(vpa(Speed))*60/(2*pi));         临界转速

zhoujunbo 发表于 2010-1-13 22:00

调用函数
1.
function y=Ms(u,l,r)
y=Msr(u,l,r)+Mst(u,l);

%u,l,r,r1分别为轴段单位长度质量、轴段长度、轴段外径、内径
2.
function Msr=Msr(u,l,r)
Msr=u*r^2/(120*l)*[36 3*l -36 3*l
    3*l 4*l^2 -3*l -l^2
    -36 -3*l 36 -3*l
    3*l -l^2 -3*l 4*l^2];
3.
function Mst=Mst(u,l)
Mst=u*l/420*[156 22*l 54 -13*l
    22*l 4*l^2 13*l -13*l^2
    54 13*l 156 -22*l
    -13*l -3*l^2 -22*l 4*l^2];
4.
function Js=Js(u,l,r)
Js=u*r^2/(60*l)*[36 3*l -36 3*l
    3*l 4*l^2 -3*l -l^2
    -36 -3*l 36 -3*l
    3*l -l^2 -3*l 4*l^2];
5.
function Ks=Ks(l,r,RotorE)
I=pi*(r^4)/4;
Ks=RotorE*I/l^3*[12 6*l -12 6*l
    6*l 4*l^2 -6*l 2*l^2
    -12 -6*l 12 -6*l
    6*l 2*l^2 -6*l 4*l^2];

zhoujunbo 发表于 2010-1-13 22:22

程序是参考钟一谔 转子动力学 有限单元法---第八章
如图所示,程序求解的是8个轴段(无圆盘),弹性支承分别位于第2、8个节点的轴系临界转速
在Matlab中当节点数超过6个的时候,带符号的频率方程就求不出来,这时候该怎么处理?(主程序最后几行)

另外,当不考虑回转矩阵的时候,程序最后几行改为
DB=B\D;                                                                                 %动力矩阵
=eig(DB);                                                                      %求解特征值 和特征向量
v=sqrt(W2)*60/(2*pi)                                                            %临界转速
求出来的跟论坛出品的求临界转速的P文件得到的结果不一致
请各位帮忙看看程序本身有什么问题,谢谢!

zhoujunbo 发表于 2010-1-13 22:28

Q 5752830
希望多跟大家一起交流!

zhoujunbo 发表于 2010-1-15 15:51

世态炎凉啊……

也可能是问题弱智,高手懒得搭理:@(

ChaChing 发表于 2010-1-15 18:54

LZ的问题应属专业方面, 要刚好遇到同专业的高手吧!
不然那么长的程序看起来非常花时间吧

[ 本帖最后由 ChaChing 于 2010-1-15 18:58 编辑 ]

woppy 发表于 2010-1-18 21:54

回复 楼主 zhoujunbo 的帖子

楼主用传递矩阵法,试试看吧

zhoujunbo 发表于 2010-1-19 17:20

原帖由 woppy 于 2010-1-18 21:54 发表
楼主用传递矩阵法,试试看吧
传递矩阵法不能考虑陀螺效应吧?
求临界转速可以,无法求出对数衰减率,该怎么表征系统稳定性呢?
望指点!

woppy 发表于 2010-1-19 19:35

回复 9楼 zhoujunbo 的帖子

楼主的问题和我的一样啊,我现在想用张文的《转子动力学理论基础》里面提到的广义Jacobi法编个程序试试了,MATLAB的立面的EIG函数,对反对称矩阵求出的特征值精度不高。

zhoujunbo 发表于 2010-1-19 20:30

回复 9楼 zhoujunbo 的帖子

张文的书没看过,广义Jacobi法可以解决这问题吗?

toes 发表于 2010-1-20 12:10

回复 9楼 zhoujunbo 的帖子

很久没做这方面东西了,我可以肯定的是:传递矩阵法考虑了陀螺效应(你指的是陀螺力矩引起的效应吗?)。
求对数衰减率,如何表征系统稳定性我没做过。

davesnw 发表于 2010-1-30 22:56

用 ansys直接解决 不用那么费劲

vehicle 发表于 2010-3-5 11:47

原帖由 davesnw 于 2010-1-30 22:56 发表
用 ansys直接解决 不用那么费劲
ansys怎么求考虑陀螺效应的模态呢?跟不考虑陀螺效应的求法区别是什么?

davesnw 发表于 2010-3-5 18:26

ansys 在计算的时候已经考虑陀螺矩阵了
页: [1] 2 3
查看完整版本: 求助——钟一谔转子动力学有限元法