UEL单元开发教程 | 三节点梁单元
知乎、B站:[易木木响叮当]
关注可了解更多的有限元数值仿真技巧。问题或建议,请公众号留言;
如果你觉得木木同学对你有帮助,欢迎赞赏。
本次给大家带来的主要内容是:如何使用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案例,本构关系以增量形式编写,便于对非线性截面行为一般化。
想要深入了解的同学也可将程序推广到三维分析中,至于非线性分析可能有些许困难,诸君努力,木木先撤了~~
本次分享仅限于此了,欢迎大家点赞收藏转发!
谢谢你看完木木同学的分享,今日份阅读花费的流量+1M哈哈哈哈哈哈。
-End-
易木木响叮当
想陪你一起度过短暂且漫长的科研生活