UEL单元开发（2）——一维杆单元（内含大福利！）

易公子 2022年5月20日浏览：1975 收藏：1

木木做单元开发也有一段日子了，在编写XFEM的单元时，遇到了瓶颈，打算先告一段落，将自己的一些经验分享出来给大家，希望帮助大家快速入门，在以后的推文中也会以单元开发为主，就像之前讲解Umat一样，用自己通俗易懂的语言讲述下来。

本篇推文由两部分组成：「其一」，以「一维杆单元」UEL子程序为主，从理论到代码Step-by-Step向大家讲解，并于Abaqus中标准单元做对比，验证代码的正确性；「其二」，木木受邀于北鲲云超算平台做主讲人，将在5月24日晚上七点于B站进行大约一个小时的讲演汇报，带着大家熟悉Abaqus的基本操作，基于Abaqus的扩展有限元仿真以及如何正确使用超算平台，直播过程中还有大量惊喜礼品、丰富算力金进行抽奖~欢迎大家积极参与。

一维杆单元UEL子程序

一维杆单元子程序的代码应该是UEL中最简单的了吧，就好像线弹性Umat那样，囊括了一些UEL的基本思想，可对比之前的二维弹簧单元子程序来学习。

理论讲解

弹性力学中利用「最小势能原理」（或虚功原理）来得到单元刚度方程，具体操作可查看：

❝

《有限元基础教程》——曾攀

《北航弹性力学博后讲透弹性力学》——@兵心依旧022（B站Up）

❞

单元刚度矩阵：

其中：为杨氏模量，为杆的横截面积，为杆的长度。UEL需要的「核心」就是求得单元刚度矩阵，「残余力RHS」可以认为是0（Cohesive单元除外），以后会给大家普及国际上流行的PPR内聚力单元。

代码讲解

篇幅原因，INP文件可在后台回复“杆单元”，即可获得，这里仅给出「关键行」：

*User element, nodes=2, type=U1001, properties=2, coordinates=2, variables=4
1,2,6
*Element, type=U1001
 1,1,2
*Elset, elset=Set-1
 1
*Uel property, elset=Set-1
 2.0E11,  0.02

程序代码如下：

      SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
     1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,
     2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,
     3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),
     1 SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL),
     2 DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*),
     3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
     4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)

      INTEGER I,J
C     获取材料属性
      Y_M = PROPS(1)
      AREA = PROPS(2)
C     获取节点坐标计算单元长度
      LEN_EL = ABS(COORDS(1,2)-COORDS(1,1))
      write(6,*) LEN_EL
C     构造刚度矩阵
      K_EL = Y_M*AREA/LEN_EL
      AMATRX(1,1) = K_EL
      AMATRX(1,2) = -K_EL
      AMATRX(2,1) = -K_EL
      AMATRX(2,2) = K_EL
C     求RHS

      DO I = 1,NDOFEL
          RHS(I,1) = 0
          DO J = 1,NDOFEL
              RHS(I,1) = RHS(I,1) - AMATRX(I,J)*U(J)
          ENDDO
      ENDDO
      
      RETURN
      END