2D梁的模态分析

以下面的问题为例，讨论FEM进行梁的模态分析：

单跨梁受约束如图所示，梁的密度ρ=3e3kg/m³，试求梁的自振频率及振型。设梁的弹性模量为E=200GPa,横截面积A=1.0e-2㎡，长度L=2，界面惯性矩为I=1E-4m⁴。

从题意，不考虑轴向拉伸，因此，刚度矩阵和质量矩阵为4*4的矩阵（对于6*6矩阵，见文章有限元法中，求解2D梁问题（http://www.cnblogs.com/heaventian/articles/2766227.html）

下面给出其计算的程序：

首先，给定单元划分，单元编号，节点的单元及整体编号的子程序

coord2DBeam.m

% THE coordinate index is from 1 to Nx(from left to right)
% 2D beam
xmin=0;xmax=2;
Nx=100;
x=linspace(xmin,xmax,Nx);
Coord=x';
save coordinates.dat Coord -ascii

% node index of each element
vertices=[(1:Nx-1)',(2:Nx)'];
save elements2DBeam.dat vertices -ascii

% boundary condition of flexure and angle of rotation
boundary=[2*Nx-1;2*Nx];
save flexure.dat boundary -ascii

然后主程序fem.m

clear
close all
%  Read the nodal coordinate data file.
load coordinates.dat;
%  Read the triangular element data file.
load elements2DBeam.dat;
E=200e9;
A=1e-2;
I=1e-4;
l=2;% length of beam
rho= 3e3;
%  Read the boundary condition data file of flexure and angle of rotation.
  load flexure.dat;
  
K = sparse ( length(coordinates)*2, length(coordinates)*2);
M = sparse ( length(coordinates)*2, length(coordinates)*2);
% Assemble global stiffness
for j = 1 : size(elements2DBeam,1)
    eleindex=reshape([2*elements2DBeam(j,:)-1;2*elements2DBeam(j,:)],4,1);
    lel=abs(coordinates(elements2DBeam(j,2))-coordinates(elements2DBeam(j,1)));
        K(eleindex,eleindex) = K(eleindex,eleindex)...
      + Beam2D2Node_Stima(E,I,lel);
        M(eleindex,eleindex) =M(eleindex,eleindex)...
      + Beam2D2Node_Mass(rho,A,lel);
end

%  Determine which nodes are associated with flexure and angle of rotation conditions.
  BoundNodes = unique ( flexure );
%  Compute the solutionfor the remaining unknown values of U.
%
  FreeNodes = setdiff ( 1:size(coordinates,1)*2, BoundNodes );

%eig(full(K(FreeNodes,FreeNodes)*inv(M(FreeNodes,FreeNodes))))
[V,D]=eig(full(K(FreeNodes,FreeNodes)),full(M(FreeNodes,FreeNodes)));

主程序调用下面两个代码：

生成单元刚度矩阵的程序Beam2D2Node_Stima.m：

function k = Beam2D2Node_Stima(E,I,L)
%determines the local stiffness matrix
%Young's module E,moment of inertia I, length L of beam element
% Notice that displacement along axis direction is not considered here
%Output: k(4×4)
k=E*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 M = Beam2D2Node_Stima(rho,A,l)
%determines the local stiffness matrix
%density rho, Area A,length l of beam element
% Notice that displacement along axis direction is not considered here
%Output: M(4×4)
M=rho*l*A/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
];

对于上面的程序，给出了前五阶固有频率：

a=sqrt(diag(V));

a(1:5)

1.0e+04 *

    0.0718

    0.4498

    1.2594

    2.4679

    4.0796

并且给出了其对应的振型：

挠度的振型：

subplot(5,1,1), plot(coordinates,[V(1:2:end,1);0])
subplot(5,1,2), plot(coordinates,[V(1:2:end,2);0])
subplot(5,1,3), plot(coordinates,[V(1:2:end,3);0])
subplot(5,1,4), plot(coordinates,[V(1:2:end,4);0])
subplot(5,1,5), plot(coordinates,[V(1:2:end,5);0])

弯角的振型：

 subplot(5,1,1), plot(coordinates,[V(2:2:end,1);0])
 subplot(5,1,2), plot(coordinates,[V(2:2:end,2);0])
subplot(5,1,3), plot(coordinates,[V(2:2:end,3);0])
subplot(5,1,4), plot(coordinates,[V(2:2:end,4);0])
subplot(5,1,5), plot(coordinates,[V(2:2:end,5);0])