[FEM][MATLAB][有限元] FEM Modal Analysis Programming with MATLAB (Frame Elements) (框架单元模态分析编程)

坚持实干、坚持一线、坚持积累、坚持思考,坚持创新。

接着博文《[动力学] 模态分析基础 (Basics of modal analysis)》继续介绍框架单元用于模态分析的方法。该部分内容也是 书本 《有限单元法:编程与软件应用》模态分析章节的部分内容节选。

10.3.1 梁单元质量矩阵

集中质量矩阵(局部坐标)

设杆材料密度为 ,单元长度为 ,截面积为 ,每个节点分担单元1/2的平动质量,无转动惯量,则单元质量矩阵

           (10.3‑1)

10.3.2 算例:2D框架模态分析

算例采用与2D框架结构静力分析中相同的结构,采用欧拉梁单元,材料密度为2.5493e-9t/mm3,采用集中质量矩阵,将单元质量集中于两端节点,且不考虑集中质量后的节点转动质量。

由于算例模态分析的Matlab代码与前面章节中静力分析的Matlab代码大体相同,因此这里仅给进行模态分析所需的新增代码。

% Truss 2D modal analysis

% Author JiDong Cui(崔济东),Xuelong Shen(沈雪龙)

% Website : www.jdcui.com

% 20151114

……

此前代码与静力分析代码相同,不再赘述。

% Properties definition

……

p=2.5493     %Density

……

定义单元属性时,增加单元密度ρ=2.5493,即重度25kN/m3,用于计算单元质量。

% Global stiffness matrix and mass matrix

stiffness = zeros(numDOF);

mass=zeros(numDOF);

for i=1:numEle;

    % Index of the element nodes

    noindex = EleNode(i,:);

    % Element length

    deltax = xx( noindex(2)) - xx( noindex(1));

    deltay = yy( noindex(2)) - yy( noindex(1));

    L = sqrt( deltax*deltax + deltay*deltay );

    % Element stiffness matrix

    C = deltax / L;

    S = deltay / L;

    T = [ C S 0 0 0 0;

         -S C 0 0 0 0;

          0 0 1 0 0 0;

          0 0 0 C S 0;

          0 0 0 -S C 0;

          0 0 0 0 0 1

    ];

    % Section properties

    E = EList(EleSect(i));

    A = AList(EleSect(i));

    I = IList(EleSect(i));

    EAL = E*A/L;

    EIL1 = E*I/L;

    EIL2 = 6.0*E*I/(L*L);

    EIL3 = 12.0*E*I/(L*L*L);

    % Element stiffness matrix in local coordinate system

    eleKo = [ EAL   0    0    -EAL   0   0;

             0   EIL3  EIL2    0  -EIL3  EIL2;

             0   EIL2  4*EIL1  0  -EIL2  2*EIL1;

             -EAL   0    0     EAL  0     0;

             0   -EIL3 -EIL2    0   EIL3  -EIL2;

             0   EIL2  2*EIL1  0  -EIL2 4*EIL1

             ];   

    % Element mass matrix in local coordinate system

    eleMo=p*A*L/2*diag([1,1,0,1,1,0]);

    % Corresponding freedom of the element nodes

    eleDof = [noindex(1)*3-2 noindex(1)*3-1 noindex(1)*3 noindex(2)*3-2 noindex(2)*3-1 noindex(2)*3];

    % Integration of element stiffness matrix and mass matrix

    stiffness( eleDof , eleDof) = stiffness( eleDof,eleDof) + T'*eleKo*T;

    mass( eleDof , eleDof) = mass( eleDof , eleDof) + T'*eleMo*T;

end

本段代码的作用是集成结构整体刚度矩阵和整体质量矩阵。本结构共有numDOF个自由度,因此首先采用zeros()方法定义了一个numDOF×numDOF的零矩阵mass。单元的质量矩阵eleM0采用集中质量形式求得,并采用与刚度矩阵相同的“对号入座法”将单元质量矩阵组装到结构整体质量矩阵mass中。

% Get the eigenvalues and eigenvectors

activeDof = setdiff( [1:numDOF]',restrainedDof);

if det(mass(activeDof,activeDof))>0

    A=mass(activeDof,activeDof)\stiffness(activeDof,activeDof);

    [Vo,Do] = eig(A);

    fo=diag(sqrt(Do)/(2*pi));

    Vf=sortrows([Vo',fo],size(Vo,1)+1);

    V=Vf(:,1:size(Vo,1))';   % Sort eigenvectors

elseif det(stiffness(activeDof,activeDof))>0

    A=stiffness(activeDof,activeDof)\mass(activeDof,activeDof);

    [Vo,Do] = eig(A);

    fo=power(diag(sqrt(Do)),-1)/(2*pi);

    Vf=sortrows([Vo',fo],size(Vo,1)+1);

    V=Vf(:,1:size(Vo,1))';   % Sort eigenvectors

end

f=sort(fo)

本段代码通过求解标准特征值问题,求得特征值和特征向量。Matlab内置函数eig([A])可直接求解特征值和特征向量。由于本例不考虑节点的转动质量,因此 det(mass(activeDof,activeDof)) 等于零,所以本例做了一个判断“if…else if…”,最终采用的是 [A]=[K]-1[M] 。求得的特征值放在向量{Do}中,相应的特征向量放在矩阵[Vo]中。其后的代码用于处理求解得到的数据,其中代码“fo=”用于根据特征值求得频率,代码“V=”用于将特征向量排序,代码“f=”用于将结构自振频率进行排序。图 10‑7所示是该桁架结构的前三阶振型。

(a)第1振型

(b)第2振型

(c)第3振型

  • 相关内容 ( Related Topics )

[00]. [Book][书] 有限单元法:编程与软件应用

[01]. 2D Truss by Matlab [Matlab 2D桁架有限元] (code by myself)

[02]. 3D Truss by MATLAB[MATLAB 3D 桁架静力分析] (code by myself)

[03]. [FEM][有限元][编程][Matlab][Code by myself] FEM Analysis: 2D Truss Element [有限元分析: 2D桁架单元]

[04]. [FEM][有限元][编程][Matlab][Code by myself] FEM Analysis: 2D Euler Beam Element [有限元分析: 2D欧拉梁单元]

[05]. [FEM][有限元][编程][Matlab][Code by myself] 2D剪切梁单元

[06]. [FEM][有限元][编程][Matlab][Code by myself] 2D Timoshenko梁单元

[07]. [FEM][有限元][编程][Matlab][Code by myself] 三角形常应变单元(CST)

[08]. [FEM][有限元][编程][Matlab][Code by myself] 平面4变形单元(Q4)

[09]. [FEM][有限元][编程][Matlab][Code by myself] 4节点四面体单元(TET4)

[10]. [FEM][有限元][编程][Matlab][Code by myself] 8节点六面体单元(C3D8)

[11]. [FEM][有限元][编程][Matlab][Code by myself] 平面8节点二次“完全积分”单元(CPS8)

[12]. [FEM][有限元][编程][Matlab][Code by myself] 平面6节点二次“完全积分”单元(CPS6)

[13]. [FEM][MATLAB][有限元] FEM Modal Analysis Programming with MATLAB (Truss Element) (桁架单元模态分析编程)

[14]. [FEM][MATLAB][有限元] FEM Modal Analysis Programming with MATLAB (Frame Elements) (框架单元模态分析编程)

[15]. [力学][有限元]Basics of Linear Buckling Analysis [线性曲屈分析基础]

[16]. [FEM][MATLAB][有限元][编程] FEM Buckling Analysis Programming with MATLAB (Truss Elements) (桁架单元曲屈分析编程)

[17]. [FEM][MATLAB][有限元][编程] FEM Buckling Analysis Programming with MATLAB (Frame Elements) (框架单元曲屈分析编程)


5 0

  • 注释 ( Comments )

  ( 如果您发现有错误,欢迎批评指正。邮箱:jidong_cui@163.com . 如果您喜欢这篇博文,请在上面给我 点个赞 吧! :-)   :-)   

  ( If you found any mistakes in the post, please let me know. Email : jidong_cui@163.com. If you like this posts, please give me a thumbs up rating on the above button! )

  • 微信公众号 ( Wechat Subscription)

WeChat_QRCode

欢迎关注 “结构之旅” 微信公众号

1 thought on “[FEM][MATLAB][有限元] FEM Modal Analysis Programming with MATLAB (Frame Elements) (框架单元模态分析编程)

  1. Pingback: [FEM][MATLAB][有限元] FEM Modal Analysis Programming with MATLAB (Truss Element) (桁架单元模态分析编程) | 崔济东的博客 - JiDong Cui's Blog

Leave a Reply

Your email address will not be published. Required fields are marked *

:wink: :-| :-x :twisted: :) 8-O :( :roll: :-P :oops: :-o :mrgreen: :lol: :idea: :-D :evil: :cry: 8) :arrow: :-? :?: :!:

This site uses Akismet to reduce spam. Learn how your comment data is processed.