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

zhaoqingfifa

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

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

1. %%
   
2. clc;
   
3. clear;
   
4. Nshaft=8;                                                                  %轴段数量；
   
5. %%
   
6. RotorE = 2.095e11;                                                         %转子弹性模量；
   
7. RotorM = 7.85e3;                                                           %转子材料密度
   
8. ShaftL = [0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25];                        %各轴段长度
   
9. ShaftDI = ones(1,Nshaft)*0.1;                                              %各轴段外径
   
10. ShaftDO = ones(1,Nshaft)*0.0;                                              %各轴段内径；
    
11. %%
    
12. LocationF=[1,9];                                                           %支承所在节点编号；
    
13. SPstiffness=[1e9,1e9];                                                     %支承刚度;
    
14. SPtype=1;SPtype= 1;                                                        %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
    
15. %%
    
16. addtionN = [];                                                             %附加轮盘编号
    
17. addtionM = [];                                                             %附加轮盘质量
    
18. addtionJ = [];                                                             %附加轮盘转动惯量
    
19. %%
    
20. U=pi*(ShaftDI.^2-ShaftDO.^2)*RotorM;                                       %单位长度质量
    
21. a=[U;ShaftL;ShaftDI;ShaftDO];
    
22. A=a.';
    
23. %%
    
24.                                                                            %整体质量矩阵计算
    
25. for i=1:Nshaft                                                                                 
    
26.     u=A(i,1);
    
27.     l=A(i,2);
    
28.     r=A(i,3);
    
29.     M(:,:,i)=Ms(u,l,r);
    
30. end;
    
31. 
    
32. B=zeros(Nshaft*2+2,Nshaft*2+2);
    
33. %整体质量矩阵第1、2行
    
34. B(1:2,1:2)=B(1:2,1:2)+M(1:2,1:2,1);
    
35. B(1:2,3:4)=B(1:2,3:4)+M(1:2,3:4,1);
    
36. %整体质量矩阵第3至Nshaft*2行
    
37. for i=3:2:Nshaft*2
    
38.     j=i-2;
    
39.     B(i:(i+1),j:(j+1))=B(i:i+1,j:j+1)+M(3:4,1:2,(i-1)/2);
    
40.     j=i;
    
41.     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);
    
42.     j=i+2;
    
43.     B(i:i+1,j:j+1)=B(i:i+1,j:j+1)+M(1:2,3:4,(i+1)/2);
    
44. end
    
45. %整体质量矩阵第Nshaft*2+1、Nshaft*2+2行
    
46. 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);
    
47. 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);
    
48. %%                                                               
    
49.                                                                            %整体回转矩阵计算
    
50. for i=1:Nshaft
    
51. u=A(i,1);
    
52. l=A(i,2);
    
53. r=A(i,3);
    
54. J(:,:,i)=Js(u,l,r);
    
55. end;
    
56. C=zeros(Nshaft*2+2,Nshaft*2+2);
    
57. %整体回转矩阵第1、2行
    
58. C(1:2,1:2)=C(1:2,1:2)+J(1:2,1:2,1);
    
59. C(1:2,3:4)=C(1:2,3:4)+J(1:2,3:4,1);
    
60. %整体回转矩阵第3-Nshaft*2行
    
61. for i=3:2:Nshaft*2                                                         
    
62.     j=i-2;
    
63.     C(i:(i+1),j:(j+1))=C(i:i+1,j:j+1)+J(3:4,1:2,(i-1)/2);
    
64.     j=i;
    
65.     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);
    
66.     j=i+2;
    
67.     C(i:i+1,j:j+1)=C(i:i+1,j:j+1)+J(1:2,3:4,(i+1)/2);
    
68. end
    
69. %整体回转矩阵第Nshaft*2+1、Nshaft*2+2行
    
70. 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);
    
71. 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);
    
72. %%                                                               
    
73.                                                                            %整体刚度矩阵的计算     
    
74. for i=1:Nshaft
    
75.     l=A(i,2);
    
76.     r=A(i,3);
    
77.     K(:,:,i)=Ks(l,r,RotorE);
    
78. end;
    
79. D=zeros(Nshaft*2+2,Nshaft*2+2);
    
80. %整体刚度矩阵第1、2行
    
81. D(1:2,1:2)=D(1:2,1:2)+K(1:2,1:2,1);
    
82. D(1:2,3:4)=D(1:2,3:4)+K(1:2,3:4,1);
    
83. %整体刚度矩阵第3至Nshaft*2行
    
84. for i=3:2:Nshaft*2
    
85.     j=i-2;
    
86.     D(i:(i+1),j:(j+1))=D(i:i+1,j:j+1)+K(3:4,1:2,(i-1)/2);
    
87.     j=i;
    
88.     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);
    
89.     j=i+2;
    
90.     D(i:i+1,j:j+1)=D(i:i+1,j:j+1)+K(1:2,3:4,(i+1)/2);
    
91. end
    
92. %整体刚度矩阵第Nshaft*2+1、Nshaft*2+2行
    
93. 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);
    
94. 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);
    
95.                                                                            %叠加支承刚度的影响
    
96. Size=size(LocationF);
    
97. for i=1:1:Size(2)
    
98.     k=2*LocationF(i)-1;
    
99.     D(k,k)=D(k,k)+SPstiffness(i);
    
100. end
     
101. %%                                                               
     
102.                                                                            %论坛程序计算结果
     
103. CSN = 5;
     
104. CriticalSpeeds_forum=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)
     
105. %%
     
106. DB=inv(B)*D;                                                               %特征矩阵
     
107. [P,w2]=eig(DB);                                                            %求解特征值和特征向量
     
108. v=sqrt(w2)*60/(2*pi);                                                      %临界转速
     
109. for i=1:1:5
     
110.     CriticalSpeed(i)=v(2*Nshaft+2-i+1,2*Nshaft+2-i+1);
     
111. end
     
112. CriticalSpeed'

复制代码

                               
登录/注册后可看大图

zhaoqingfifa

但是求解的结果与论坛给定的相差比较的多，希望各位能帮忙看看问题出在哪里？
论坛结果：
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

问题已找到，
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!

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!

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