UEL单元开发教程 | 三节点梁单元

知乎、B站[易木木响叮当]
关注可了解更多的有限元数值仿真技巧。问题或建议,请公众号留言;
如果你觉得木木同学对你有帮助,欢迎赞赏。

本次给大家带来的主要内容是:如何使用UEL子程序中开发三节点平面梁单元?

关于梁单元的介绍,我在之前的推文中介绍过如何用Matlab编制梁单元有限元程序,及两节点的梁单元UEL子程序编制,感兴趣的朋友可以点击跳转阅读浏览。

要点

为防止轴力过大出现“过度的刚性行为”,常常在梁单元内部增加一个节点,发挥“缓和”作用,增加的这个节点只有切线自由度,即切线位移,单元刚度矩阵的求解需要用到高斯数值积分

UEL单元开发教程 | 三节点梁单元的图1
三节点梁单元示意图

如上图所示,三节点梁单元共有 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
  1. 计算单元长度 ,单元转角

  2. 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即可

UEL单元开发教程 | 三节点梁单元的图2
DU官方培训手册解读

    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)
  1. isvint:在当前积分点获取SDV;

  2. 使用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
  1. 计算单元内力矢量: ,目前阶段的UEL分析中不涉及等效节点荷载处理的情况, 均为0;

  2. 离散化:

    为高斯积分系数。
 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 个自由度,左端固支,右端承受向上的集中载荷

UEL单元开发教程 | 三节点梁单元的图3
梁模型

Abaqus 计算结果如下,由云图可知RF2的最大值为 1,也就是施加的载荷 大小。

UEL单元开发教程 | 三节点梁单元的图4
Abaqus RF2值

资源获取:如需UEL源程序及INP文件,可联系公众号【易木木响叮当】后台。

以上就是本次分享的三节点梁单元UEL案例,本构关系以增量形式编写,便于对非线性截面行为一般化。

想要深入了解的同学也可将程序推广到三维分析中,至于非线性分析可能有些许困难,诸君努力,木木先撤了~~


本次分享仅限于此了,欢迎大家点赞收藏转发!


谢谢你看完木木同学的分享,今日份阅读花费的流量+1M哈哈哈哈哈哈UEL单元开发教程 | 三节点梁单元的图5



-End-


易木木响叮当

想陪你一起度过短暂且漫长的科研生活

(2条)
默认 最新
学习了
评论 点赞 1
感谢分享
评论 点赞
点赞 2 评论 2 收藏
关注