ABAQUS UEL二次开发（用户自定义单元体力施加）

概述：研究了平面应力四边形单元的体力施加问题，分别采用四种方式给CPS4单元施加体力，为：

（1）采用ABAQUS自带的体力项施加体力

（2）采用DLOAD子程序施加体力

（3）将体力转化为节点力施加体力

（4）编写UEL，并通过JDLTYP接口施加体力

在动力隐式计算中对比上述计算结果，最终结果一致。

注：

（）本帖子参考了如下链接：ABAQUS用户定义单元UEL与VUEL从入门到放弃系列3————讲一讲的RHS和NBLOCK-技术邻 (jishulink.com)，以及该帖子评论区的讨论：

经验证，该发帖者、以及评论区这位评论者的说法是正确的，供大家参考。

（）本帖子附件中的obj文件不可编辑，仅适用于该帖子对应的inp文件。

（）模型信息

几何尺寸为2x2，弹性模量1e10，密度2000，泊松比0.25，重力加速度取9.81，网格尺寸为1x1，共计四个单元。边界条件为右端全固定，四个单元施加体力，示意图为：

（）工况说明

对于工况一：采用ABAQUS自带的体力项施加体力，可以直接在ABAQUS中操作即可。

对于工况二：采用DLOAD子程序施加体力，需要先建立单元集合，然后在DLOAD中输出单元的体力项，这里的子程序编写非常的简单，直接计算单元的质量就行，编写的子程序为：

      SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,
     1 COORDS,JLTYP,SNAME)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION TIME(2), COORDS (3)
      CHARACTER*80 SNAME
!输出单元质量
      F=-9.81*2000
      RETURN
      END

对于工况三：将体力转化为节点力施加体力，需要将单元的总质量分到节点上面，因为这里的算例非常的简单，而且单元的形状是规则的，四个单元的尺寸都是一样的，因此可以直接将单元的总质量直接均分到四个节点上，特别注意，当网格不规则或者体力复杂时候，需要借助形函数进行积分计算每个节点的质量，这里只是为了探索流程的正确性，都是采用最简单的。在这个算例中，单元的总质量为9.81*2000=19620，则每个节点的质量为4905。将质量分配到节点上，示意图为：

可以看到，一共需要计算九个节点的质量，顶层、中层、底层各三个节点。

对于底层节点，左边第一个节点：质量为4905，  中间节点为4905*2，右边第一个节点为4905

对于中层节点，左边第一个节点：质量为4905*2，中间节点为4905*4，右边第一个节点为4905*2

对于顶层节点，左边第一个节点：质量为4905，  中间节点为4905*2，右边第一个节点为4905

对应的inp文件为：

*Cload
Part-1-1.down-1, 2, -4905.
** Name: down-2   Type: Concentrated force
*Cload
Part-1-1.down-2, 2, -9810.
** Name: down-3   Type: Concentrated force
*Cload
Part-1-1.down-3, 2, -4905.
** Name: middle-1   Type: Concentrated force
*Cload
Part-1-1.middle-1, 2, -9810.
** Name: middle-2   Type: Concentrated force
*Cload
Part-1-1.middle-2, 2, -19620.
** Name: middle-3   Type: Concentrated force
*Cload
Part-1-1.middle-3, 2, -9810.
** Name: top-1   Type: Concentrated force
*Cload
Part-1-1.top-1, 2, -4905.
** Name: top-2   Type: Concentrated force
*Cload
Part-1-1.top-2, 2, -9810.
** Name: top-3   Type: Concentrated force
*Cload
Part-1-1.top-3, 2, -4905.

对于工况四：编写UEL，并通过JDLTYP接口施加体力，需要inp与UEL子程序的相互配合。

首先形成inp文件，包括用户自定义单元的定义个荷载的定义，对应的inp文件为：

*USER ELEMENT, NODES=4, TYPE=U1008, IPROPERTIES=3, PROPERTIES=8, COORDINATES=2, VARIABLES=1200
1, 2
*ELEMENT, TYPE=U1008, ELSET=SBFEM_UEL
1, 1, 2, 5, 4
2, 2, 3, 6, 5
3, 4, 5, 8, 7
4, 5, 6, 9, 8
*UEL PROPERTY, ELSET=SBFEM_UEL
1e+10, 0.25,2000,	0, 0, 0, 0, 0

然后定义荷载，这里的荷载定义至关重要，他直接体现在UEL的JDLTYP参数中，对应的inp文件为：

*Dload
Part-1-1.gravity, U1001, 1.

这里的“U1001”中的“1001”就是UEL中JDLTYP的数值，UEL中可写程序识别这个参数，然后进行体力（外力Fext）的计算，然后提现到RHS中。

（）分析步设置

直接放inp文件吧，工况四的inp分析步设置为：其他三个工况的设置见附件。

*Step, name=Step-1, nlgeom=NO, inc=1000
*DYNAMIC, direct
0.01, 1., 
*Boundary
Part-1-1.FIXED, 1, 1
Part-1-1.FIXED, 2, 2
*Dload
Part-1-1.gravity, U1001, 1.
*Restart, write, frequency=0
*Output, field, variable=PRESELECT
*Output, history, variable=PRESELECT
*End Step

可以看到采用了动力隐式分析步，增量步长为0.01，计算总时长为1，总计增量步为100。

（）计算结果

体力荷载向位移云图为：窗口1~4即为工况且1~4的计算结果。

体力荷载向支反力云图为：窗口1~4即为工况且1~4的计算结果。

不想写字了，有兴趣的直接移步附件吧。。。。。。。。。。。。。。

附件：

-------------
-------------
ABAQUS-DLOAD
  A.FOR
  JOB-1.INP
  R.BAT
ABAQUS-NODE-LOAD
  JOB-1.INP
  R.BAT
ABAQUS-PRESSURE-LOAD
  JOB-1.INP
  R.BAT
UEL-JDLTYP
  JOB-1.INP
  job.obj
  RUN.BAT
-------------
-------------