UEL子程序开发——梁单元

易公子 2022年9月27日浏览：1720 收藏：1

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

今天木木给大家带来的是梁单元的UEL子程序，作为杆系单元重要的组成成员，该部分是得好好讲讲了，而且要全面的讲，尽木木的知识能力所及讲的透彻，让刚入手UEL的小白们看的“起劲”！参考文档和源文件可在公众号后台回复：梁单元UEL即可。

今天分享的内容主要有：梁单元介绍、等效节点荷载处理、状态变量的自定义输出、理论解对比等等，精彩超多，请进一步往下翻阅~

模型讲解

梁单元的讲解在前几期的有限元基础编程的推文中有详细讲述，感兴趣的小伙伴可以跳转相关页面查阅。

案例简介

我们今天讲的案例如图二所示，梁长 $l=100$ ，分布荷载为 $-10$ ，将左端点的U1、U2、UR3自由度束缚，右端点的U2自由度束缚。

理论解

这里我们简单讲述一下变形协调条件，图3所示（图2分解）为一承受均布荷载的悬臂梁，在B端点产生的位移 $v_{Bq}$ 在数值上可等价于由支反力 $F_B$ 产生的位移 $v_{BB}$ ，有：
$v_{Bq}+v_{BB}=0$ $\frac{ql^4}{8EI}-\frac{F_Bl^3}{3EI}=0$ $F_B=\frac{3}{8}ql$ $F_A=\frac{5}{8}ql,M_A=\frac{1}{8}ql^2$

由以上公式可得到右端点节点力、弯矩以及左端点节点力。

子程序讲解

编制矩阵操作函数

为了程序的易读性，我们可以将经常使用的的计算部分编写成一个子程序，使用时直接call一下即可，比如在UMAT中，我们会常常用到Mises应力公式，所以我们可以将Mises的求解编写为一个子程序，在UMAT中调用即可，大大减少了代码量，更能提升易读性。

本程序中，我编写了矩阵置零函数，用于数组的初始化，大家在自己的程序中编制矩阵相乘函数和矩阵求逆函数，大大提升编写代码的效率。

call K_MATRIX_ZERO(AMATRX,ndofel,ndofel)
call K_MATRIX_ZERO(RHS,ndofel,nrhs)
call K_MATRIX_ZERO(SVARS,4,1)
C---------------------矩阵置零------------------------
subroutine K_MATRIX_ZERO (A,n,m)
INCLUDE 'ABA_PARAM.INC'
dimension A(n,m)
do i=1,n
    do j=1,m
    A(i,j) = 0.d0
    end do
end do
return
end

等效荷载处理

在本次案例中，我们注意到，梁体受到向下的均布荷载，以往的UEL案例分析中，我们没考虑过这种面力、体力的作用，可以参考一下势能泛函：
$\varPi =\frac{1}{2}\int_0^L{EI\left( \frac{d^2v}{dx^2} \right)}^2dx-\int_0^L{pv}dx-\sum_m{p_mv_m-}\sum_k{M_kv_{k}^{'}}$
具体含义可见附带文档。本案例的面力处理方式：将其转化为单元节点荷载，转化的原理可参考附带文档，也可查阅有限元相关教材，推荐曾攀老师的《有限元基础教程》。
$f^e=\left[ \begin{matrix} \frac{pl_e}{2},& \frac{pl_{e}^{2}}{12},& \frac{pl_e}{2},& -\frac{pl_{e}^{2}}{12}\\ \end{matrix} \right] ^T$

H(1)=-P*DL/TWO
H(2)=-P*DL2/TWELVE
H(3)=-P*DL/TWO
H(4)=P*DL2/TWELVE
!---注意RHS的变化------------
do I=1,NDOFEL
    RHS(I,1)=RHS(I,1)-H(I)
end do

状态变量输出

为了与理论解析解做对比，可以以状态变量的形式将右端点的节点力和弯矩以及左端点的节点力输出到指定位置，这里输出到指定文件中。涉及到Fortran语言中文件的使用，编制了计算节点力的子程序，如下所示。

CALL Node_Force(SVARS,AMATRX,U,H)
!-----------------------------
SUBROUTINE  Node_Force(SVARS,AMATRX,U,H)
INCLUDE 'ABA_PARAM.INC'
DOUBLE PRECISION,DIMENSION(4):: H,SVARS,U
DIMENSION AMATRX(4,4)
SVARS(1:4)=SVARS(1:4)+MATMUL(AMATRX,U(1:4))
SVARS(1)=SVARS(1)+H(1)
SVARS(2)=SVARS(2)+H(2)
SVARS(3)=SVARS(3)+H(3)
SVARS(4)=SVARS(4)+H(4)
RETURN
END

状态变量输出

OPEN(UNIT=16,FILE='D:\TEMP\STRESS.DAT',STATUS='UNKNOWN',
1 FORM='FORMATTED')
IF(JELEM .EQ. 1) THEN
! JELEM为单元编号
    WRITE(16,*) '-----------NEXT STEP-------------'
END IF
WRITE(16,*) 'ELEMENT NUMBER',JELEM
WRITE(16,*) 'NODE',JELEM,'FORCE',SVARS(1)
WRITE(16,*) 'NODE',JELEM,'MOMENT',SVARS(2)
WRITE(16,*) 'NODE',JELEM+1,'FORCE',SVARS(3)
WRITE(16,*) 'NODE',JELEM+1,'MOMENT',SVARS(4)
WRITE(16,*) '  '