ME! 发表于 2013-6-5 16:33

本帖最后由 ME! 于 2013-6-5 16:36 编辑

这是改好的程序,两端加上轴承刚度2e7,和论坛程序一样,我看了下楼主的,是支撑刚度没有加上去,终于找到楼主的错误了,花了我好长时间Nshaft=8;                                                                  %轴段数量;
%%
RotorE = 2.095e11;                                                         %转子弹性模量;
RotorM = 7.85e3;                                                         %转子材料密度
ShaftL = ;                        %各轴段长度
ShaftDI = ones(1,Nshaft)*0.05;                                              %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0;                                              %各轴段内径;
%%
LocationF=;                                                         %支承所在节点编号;
SPstiffness=;                                                   %支承刚度;
SPtype=1;SPtype= 1;                                                      %支承类型选择, 1 -- 刚性支承, 2 -- 弹性支承;
%%
addtionN = [];                                                             %附加轮盘编号
addtionM = [];                                                             %附加轮盘质量
addtionJ = [];                                                             %附加轮盘转动惯量
%%
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                                                         
    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



DB=inv(B)*D;                                                               %特征矩阵
=eig(DB);                                                            %求解特征值和特征向量
n=sqrt(w2)*60/(2*pi);                                                      %临界转速
for i=1:1:5
    CriticalSpeed(i)=n(2*Nshaft+2-i+1,2*Nshaft+2-i+1);
end
CriticalSpeed'

子程序如下:



function Mst=Mst(u,l)
Mst=u*l/420*[156 22*l 54 -13*l
    22*l 4*l^2 13*l -3*l^2
    54 13*l 156 -22*l
    -13*l -3*l^2 -22*l 4*l^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];

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];
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];


ME! 发表于 2013-6-5 16:50

zhoujunbo 发表于 2010-1-13 22:22 static/image/common/back.gif
程序是参考钟一谔 转子动力学 有限单元法---第八章
如图所示,程序求解的是8个轴段(无圆盘),弹性支承分 ...

%%考虑陀螺力矩时,将程序改为:
E=C-B;                         %动力矩阵
=eig(D,E);               %求解特征值 和特征向量
n1=sort(diag(sqrt(W2)*60/(2*pi)))    %临界转速
得到前三阶临界转速如下:
n1 =1.0e+006 *
                  0 + 0.002705558134827i
                  0 + 0.007862569972723i
                  0 + 0.013175908450620i
稍微大于不考虑陀螺力矩情况下的临界转速:
n2 =1.0e+005 *
   0.027029841694741
   0.078489627807782
   0.130985058432539

1713573225 发表于 2015-12-2 21:04

ME! 发表于 2013-6-5 16:50
%%考虑陀螺力矩时,将程序改为:
E=C-B;                         %动力矩阵
=eig(D,E);                ...

求涡动频率是符号解算不出来 18*18 怎么用数值解呀

空空道人 发表于 2016-1-11 20:02

ME! 发表于 2013-6-5 16:50
%%考虑陀螺力矩时,将程序改为:
E=C-B;                         %动力矩阵
=eig(D,E);                ...

能否请教下算瞬态的程序啊,谢谢

阿H 发表于 2016-5-11 10:39

请问 阻尼如何考虑?谢谢!
页: 1 2 [3]
查看完整版本: 求助——钟一谔转子动力学有限元法