求助有限元法求解转子临界转速

zhaoqingfifa · 发表于 2013-1-23 23:49

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

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

x

利用有限元法对光轴进行求解，采用的程序是在论坛之前水友的程序基础上修改的，源程序可见http://forum.vibunion.com/thread-89771-1-1.html

在其基础之上进行的修改。
程序求解的是8个轴段（无圆盘），弹性支承分别位于第2、8个节点的轴系临界转速
程序如下：

%%
clc;
clear;
Nshaft=8; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = [0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25]; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
%%
LocationF=[1,9]; %支承所在节点编号；
SPstiffness=[1e9,1e9]; %支承刚度;
SPtype=1;SPtype= 1; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
%%
addtionN = []; %附加轮盘编号
addtionM = []; %附加轮盘质量
addtionJ = []; %附加轮盘转动惯量
%%
U=pi*(ShaftDI.^2-ShaftDO.^2)*RotorM; %单位长度质量
a=[U;ShaftL;ShaftDI;ShaftDO];
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
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行
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
%整体刚度矩阵第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);
%叠加支承刚度的影响
Size=size(LocationF);
for i=1:1:Size(2)
k=2*LocationF(i)-1;
D(k,k)=D(k,k)+SPstiffness(i);
end
%%
%论坛程序计算结果
CSN = 5;
CriticalSpeeds_forum=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)
%%
DB=inv(B)*D; %特征矩阵
[P,w2]=eig(DB); %求解特征值和特征向量
v=sqrt(w2)*60/(2*pi); %临界转速
for i=1:1:5
CriticalSpeed(i)=v(2*Nshaft+2-i+1,2*Nshaft+2-i+1);
end
CriticalSpeed'

复制代码

登录/注册后可看大图

zhaoqingfifa · 发表于 2013-1-23 23:51

但是求解的结果与论坛给定的相差比较的多，希望各位能帮忙看看问题出在哪里？
论坛结果：
CriticalSpeeds_forum =
1.0e+004 *
0.3041
1.2138
2.7234
4.8289
7.5341
改写的程序结果：
ans =
1.0e+004 *
0.5834
2.0561
3.7697
5.5592
8.1540

zhaoqingfifa · 发表于 2013-1-24 22:39

本帖最后由 zhaoqingfifa 于 2013-1-24 22:47 编辑

问题已找到，
CriticalSpeeds_forum =

1.0e+004 *

0.3041
1.2138
2.7234
4.8289
7.5341

改写的程序结果
CriticalSpeed =

1.0e+004 *

0.3033
1.2016
2.6620
4.6314
7.0321
ANSYS计算结果为：
3030
11980
26450
45855
69444
结果基本算是对应上了，但是貌似有限元法计算的结果偏小，希望高手给予解答！

ME! · 发表于 2013-6-5 10:19

本帖最后由 ME! 于 2013-6-5 10:33 编辑

zhaoqingfifa 发表于 2013-1-24 22:39

登录/注册后可看大图

问题已找到，
CriticalSpeeds_forum =

请问错误处在哪里能说明一下吗，我运行你的程序还是第一轮的结果，那应该是没有改过的吧好像是因为论坛程序的内径是直径，而你的是u为半径，

但是我用ansys计算好像得不到你的结果，请问你是用什么单元吗？
我没有考虑转动惯性矩阵的影响matlab程序和ansys beam3单元的计算结果能对应，直径为0.1，半径为0.05
matlab计算结果：固有频率f：
f =  1.0e+002 *
1.149803836884734
3.171213383661560
6.227643357869332

ansys计算结果：
SET TIME/FREQ LOAD STEP SUBSTEP  CUMULATIVE
   1  114.87          1       1       1
   2  315.98          1       2       2
   3  618.00          1       3       3
   4  1019.8          1       4       4
   5  1299.8          1       5       5
   6  1523.0          1       6       6
   7  2128.7          1       7       7
   8  2649.8          1       8       8
   9  2802.3          1       9       9
10  3908.5          1       10       10

ME! · 发表于 2013-6-5 10:46

请看http://forum.vibunion.com/forum. ... mp;extra=#pid719765不知道是不是同一个人，好像是叠加轴承刚度的那里有问题

账号		自动登录	找回密码
密码			我要加入

[转子动力学] 求助有限元法求解转子临界转速

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