有限元基础编程——梁单元

引言

本次推文分享的是梁单元的有限元基本格式以及如何使用Matlab对其编程。另外说一下，后处理中的网格绘制以及云图显示不属于有限元核心部分，将可能会在后期安排内容推送，或者读者自行在网上下载相关代码进行改动(Copy)。推文主要有以下内容：1. 平面纯弯梁单元描述(节点无轴向位移)及其Matlab编程；2. 一般平面梁单元描述(节点无轴向位移)及其Matlab编程；3. 常见梁单元等效节点荷载。特别说明：理论内容及主要代码来自曾攀——《有限元分析基础教程》。

平面纯弯梁单元描述

有限元基本格式描述

节点位移    、节点力    

位移场：

其中，    为梁单元内任一点位移，    为单元形函数，注意到这里已经用到了“归一化”思想，为以后讲解等参单元做基础。 

应变场：

其中，    为单元的几何矩阵，    为所测点以中性层为起点的y方向的坐标。 

应力场：

其中，    为弹性模量，    为单元的应力矩阵。 

单元刚度矩阵：

单元刚度方程：

模型描述

Matlab编程：

function k =Beam1D2Node_Stiffness(E,I,L)
% 直接组装梁单元刚度
k = E*I/(L*L*L)*[12 6*L -12 6*L
                6*L 4*L*L -6*L 2*L*L
                   -12 -6*L 12 -6*L
                   6*L 2*L*L -6*L 4*L*L];
function z = Beam1D2Node_Assembly(KK,k,i,j)
% 该函数进行整体刚度矩阵的组装
% 输入单元刚度矩阵k，单元的节点编号i、j
% 输出整体刚度矩阵KK 
DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
for n1=1:4
   for n2=1:4
      KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
   end
end
z=KK;
function  v = Beam1D2Node_Deflection(x,L,u)
% 该函数计算单元内某点的挠度
% 输入所测点距梁单元左节点的水平距离x 
% 输入梁单元的长度L，节点位移列阵u  
e=x/L;
N1=1-3*e*e+2*e*e*e;
N2=L(e-2*e*e+e*e*e);
N3=3*e*e-2*e*e*e;
N4=L(e*e*e-e*e);
N=[N1,N2,N3,N4];
v=N*u

% -----------主程序----------------
format compact
E = 200*10^9;I = 118.6*10^-6;
L1 = 5;L2 = 2.5;
% cal element stiffness
k1 = Beam1D2Node_Stiffness(E,I,L1);
k2 = Beam1D2Node_Stiffness(E,I,L2);
% stiffness assembly
KK = zeros(6,6);
KK = Beam1D2Node_Assembly(KK,k1,1,2);
KK = Beam1D2Node_Assembly(KK,k2,2,3)
% cal displacement
k = KK(4:6,4:6);
p = [39062;-31250;13021];
u = k\p
% cal node force
U = [0;0;0;u];
P = KK*U
% cal RF
F = [-62500;-52083;-93750;39062;-31250;13021];
RF = P-F

最终计算得到的支反力、节点位移值与书中解析解一致，可自行验证。

一般平面梁单元描述

一般平面梁单元在节点自由度上比纯弯梁多了一个轴向位移，故一个节点三个自由度，单元刚度矩阵是一个    的矩阵(两节点)。

总的节点荷载列阵：

单元刚度：

模型描述

首先应该先将分布荷载等效为节点荷载，然后参与计算，常见的梁单元等效节点荷载形式将会在下一小节展示。

总的节点荷载列阵：

   解题步骤与上个例题类似：单元刚度矩阵组装、整体刚度矩阵组装、位移计算、支反力计算。 

常见梁单元等效节点荷载