首页/文章/ 详情

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

易木木响叮当

1年前浏览630

本次给大家带来的主要内容是：如何使用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)
enddo

上式是对应变矩阵的编译；

对 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)
enddo
enddo
enddo

计算单元内力矢量：，目前阶段的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