【JY】有限元基础编程——梁单元

最新推荐文章于 2024-04-12 10:00:31 发布

建源之光

最新推荐文章于 2024-04-12 10:00:31 发布

阅读量396

点赞数

原文链接：https://mp.weixin.qq.com/s?__biz=MzIyNjI5NTAzNA==&mid=2247493455&idx=1&sn=2ae1a226facbdb14433051616a7955c8&chksm=e99a1396c893fc674c8ae40c9307715af2fb70cd9c5b0b175e3a12a46c4f634677d4e5fe3b44&scene=126&sessionid=0

版权

Dr.JY有话说：

本文转载自易木木的公众号易木木响叮当，欢迎关注！

易木木在有限元领域展现出深厚的专业知识和卓越的理解能力。他们对有限元分析原理和方法有深入的掌握，能够熟练运用各种有限元软件，如ABAQUS，进行复杂结构的模拟和分析。

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

引言

本次推文分享的是梁单元的有限元基本格式以及如何使用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

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