UEL单元开发教程 | 三节点梁单元
本次给大家带来的主要内容是：如何使用UEL子程序中开发三节点平面梁单元？
关于梁单元的介绍，我在之前的推文中介绍过如何用Matlab编制梁单元有限元程序，及两节点的梁单元UEL子程序编制，感兴趣的朋友可以点击跳转阅读浏览。
要点
为防止轴力过大出现“过度的刚性行为”，常常在梁单元内部增加一个节点，发挥“缓和”作用，增加的这个节点只有切线自由度，即切线位移，单元刚度矩阵的求解需要用到高斯数值积分。
如上图所示，三节点梁单元共有 7 个自由度，需要两个高斯积分点保证数值结果精度，节点顺序按照上图表示的来，中间节点为第三个节点。
有限元格式
为应变关系矩阵， 为轴向应变， 为曲率，轴向力 ，力矩 ，单元长度 。应变矩阵：
局部坐标系下的B矩阵为：
坐标转换矩阵：
通过坐标转换矩阵将局部坐标系下的 转换为整体坐标系下的矩阵 。本构关系：
单元刚度矩阵：
单元内力矢量：
代码解读
dx=coords(1,2)-coords(1,1)
dy=coords(2,2)-coords(2,1)
dl2=dx**2+dy**2
dl=sqrt(dl2)fortran
hdl=dl/two
acos=dx/dl
asin=dy/dl计算单元长度 ，单元转角 和 ；
coords(1,2)表示第 2 个节点的 坐标，相应的coords(1,1)表示第 1 个节点的 坐标。
b(1,1)=(-three+four*g)*acos/dl
b(1,2)=(-three+four*g)*asin/dl
b(1,3)=zero
b(1,4)=(-one+four*g)*acos/dl
b(1,5)=(-one+four*g)*asin/dl
b(1,6)=zero
b(1,7)=(four-eight*g)/dl
b(2,1)=(-six+twelve*g)*(-asin)/dl2
b(2,2)=(-six+twelve*g)*acos/dl2
b(2,3)=(-four+six*g)/dl
b(2,4)=(six-twelve*g)*(-asin)/dl2
b(2,5)=(six-twelve*g)*acos/dl2
b(2,6)=(-two+six*g)/dl
b(2,7)=zero以上代码用于计算应变关系 矩阵，由于Fortran语言不方便矩阵的乘积运算，可直接运算好写入程序。
如上面代码片所示，还可以使用我之前做的Fortran矩阵运算相关推文里面的矩阵相乘子程序，效果是一样的。
eps=zero ! 当前应变
deps=zero ! 应变增量
cap=zero ! 当前曲率
dcap=zero ! 曲率增量
do k=1,7
eps=eps+b(1,k)*u(k)
deps=deps+b(1,k)*du(k,1)
cap=cap+b(2,k)*u(k)
dcap=dcap+b(2,k)*du(k,1)
end do上式是对应变矩阵的编译；
对 7 个自由度进行遍历，u(k)表示单元第 个自由度的值;
du(k,1)表示第 个自由度的增量值，有关du的解释，官方培训手册给出如下解读，一般情况下KRHS取1即可
isvint=1+(kintk-1)*nsvint
bn=zero ! 轴力
bm=zero ! 弯矩
daxial=zero ! 切线轴向刚度
dcoupl=zero ! 切线刚度耦合项
call ugenb(bn,bm,daxial,dbend,dcoupl,eps,deps,cap,dcap,
1 svars(isvint),nsvint,props,nprops)isvint：在当前积分点获取SDV；使用
ugenb本构关系子程序获得轴向应变、曲率、轴向力和弯矩，存储在高斯积分点上。
do k1=1,7
rhs(k1,1)=rhs(k1,1)-hdl*(bn*b(1,k1)+bm*b(2,k1))
bd1=hdl*(daxial*b(1,k1)+dcoupl*b(2,k1))
bd2=hdl*(dcoupl*b(1,k1)+dbend*b(2,k1))
do k2=1,7
amatrx(k1,k2)=amatrx(k1,k2)+bd1*b(1,k2)+bd2*b(2,k2)
end do
end do
end do计算单元内力矢量： ，目前阶段的UEL分析中不涉及等效节点荷载处理的情况， 均为0；
离散化：
， 为高斯积分系数。
subroutine ugenb(bn,bm,daxial,dbend,dcoupl,eps,deps,cap,dcap,
1 svint,nsvint,props,nprops)
include 'aba_param.inc'
parameter(zero=0.d0,twelve=12.d0)
dimension svint(*),props(*)
h=props(1)! 截面高度
w=props(2)! 截面宽度
E=props(3)! 弹性模量
daxial=E*h*w
dbend=E*w*h**3/twelve
dcoupl=zero
bn=svint(1)+daxial*deps
bm=svint(2)+dbend*dcap
svint(1)=bn
svint(2)=bm
svint(3)=eps
svint(4)=cap
return
end以上代码片为ugenb本构子程序，具体的一些公式我还没找到，就不做细节讨论了，这里只给出子程序供大家学习参考。
问题描述
结合悬臂梁案例，运行一下本文所讲的三节点梁单元。将整根梁划分为 6 个单元，13 个节点，共计 20 个自由度，左端固支，右端承受向上的集中载荷 。
Abaqus 计算结果如下，由云图可知RF2的最大值为 1，也就是施加的载荷 大小。
资源获取：如需UEL源程序及INP文件，可联系后台。
以上就是本次分享的三节点梁单元UEL案例，本构关系以增量形式编写，便于对非线性截面行为一般化。
想要深入了解的同学也可将程序推广到三维分析中，至于非线性分析可能有些许困难，诸君努力，木木先撤了~~
本次分享仅限于此了，欢迎大家点赞收藏转发！
