一个简单的UEL梁单元（上）

一个简单的ABAQUS的UEL实例，介绍了inp文件的编写和UEL的编写主要是一个平面线弹性梁单元的UEL实现

sUaN_SuAn____

2685人浏览 · 2022-03-17 15:44:42

sUaN_SuAn____ · 2022-03-17 15:44:42 发布

出于课题需要，学习研究了abaqus uel二次开发，这里通过一个简单的平面线弹性梁单元的UEL实例来复习巩固一下学到的东西，为后续钢筋和混凝土的本构修改做准备。

主要的参考资料为：王涛老师的《基于ABAQUS的有限元子程序开发及应用》和阚前华著的《非线性本构关系在ABAQUS中的实现》推荐大家看看

图示为实例中使用的单元节点和Gauss积分点位置

添加内部节点C的目的是为了避免切向位移插值采用线性插值（ $u_a,u_b$ ）产生恒定的轴向应变，这与法向的三次插值( $v_a,\phi_a,v_b,\phi_b$ )产生的线性变化曲率不一致。内部节点C仅有一个切向自由度（ $u_c$ ）

inp文件的构成

先从inp文件讲起，代码如下：

inp代码

*Heading
** Job name:Job-beam-uel Model name: Model-beam
*Preprint, contact=NO, echo=NO, history=NO, model=NO
**
** PART
**
*Part, name=Part-1
*Node
	1,		0.,		0.
	21,		0.,		1.
*NGEN, NSET=all_node
	1,		21,		1
*User Element, type=u1, node=3, coordinates=2, variables=8
	1,		2,		6
	3,		1
*UEL Property, elset=all_element
	0.01,	0.02,	2.1e11
*Element, type=u1
	1,		1,	3,	2
*ElGEN, elset=all_element
	1,		10,		2,		1
*End Part
**
** ASSEMBLY
**
*Assembly, name=Assembly
*Instance, name=Part-1-1, part=Part-1
*End Instance
*Nset, nset=Set-Fix, instance=Part-1-1
	1,
*Nset, nset=Set-Load, instance=Part-1-1
	21,
*End Assembly
**
** STEP:step-1
**
*step, name=step-1, nlgeom=NO, inc=1000
*static
	0.01,	1.,	1e-5,	0.01
**
** BOUNDARY CONDITION
**
*Boundary
Set-Fix, 1, 2
Set-Fix, 6, 6
*Boundary
Set-Load, 1, 1, 0.1
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0 
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
U,
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, frequency=1
*node output, nset=Set-Load
rf1, rf2
*node output, nset=Set-Load
u1, u2
*End Step

有关inp文件的解读也可以参考方自虎老师的教学视频，免费而且非常细致

接下来分部分对inp文件做一个介绍：

分段介绍1：Heading

*Heading
** Job name:Job-beam-uel Model name: Model-beam
*Preprint, contact=NO, echo=NO, history=NO, model=NO
**

abaqus不区分大小写，通过 $\ lines$ 的形式构成基本语句
通过 $* *$ 表示注释行
*Heading: 表示定义分析的标题，后续的第一行的前80个字符会作为标题打印出来，也可省略为空
*Preprint: 用于选择将从分析输入文件处理器获得的打印输出，对于分析计算无影响

分段介绍2：Part

** PART
**
*Part, name=Part-1
*Node
	1,		0.,		0.
	21,		0.,		1.
*NGEN, NSET=all_node
	1,		21,		1
*User Element, type=u1, node=3, coordinates=2, variables=8
	1,		2,		6
	3,		1
*UEL Property, elset=all_element
	0.01,	0.02,	2.1e11
*Element, type=u1
	1,		1,	3,	2
*ElGEN, elset=all_element
	1,		10,		2,		1
*End Part
**

与通过abaqus创建一个完整的job一样需要定义part，这里的部件名为 Part-1 ，首先通过 *Node 方式定义梁的两端节点，再通过 *NGEN 方式，以节点 1 和节点 21 作为首尾节点，以1为节点号数字增量生成节点。
在定义节点后通过 *User Element 定义用户单元，包括：类型u1；单元包含节点数3；单元节点自由度2（平面梁单元）；UEL中需要作为数据的属性值的数量，也称状态变量8，根据积分点需要确定，下面两行data lines 第一行代表所有节点的自由度即平动1、平动2、转动6，第二行表示例外的节点，第三个节点，自由度为平动1；
通过 *UEL Property 定义用户单元的属性，梁截面的宽w、高h和材料的弹性模量E。
在定义完UEL的相关参数后，通过 *Element 生成单元1，使用节点为1、3、2，注意这里的编号顺序，该顺序与传入subroutine UEL中的节点顺序一致。
通过 *ElGEN 的方式生成单元，即以单元1为参考，以1单元内节点编号加2的方式生成10个单元，单元编号以1增长。
即生成：单元 2 节点包含3，4，5，单元 3 节点包含5，6，7…
完成部件定义，*End Part

分段介绍3：Assembly

** ASSEMBLY
**
*Assembly, name=Assembly
*Instance, name=Part-1-1, part=Part-1
*End Instance
*Nset, nset=Set-Fix, instance=Part-1-1
	1,
*Nset, nset=Set-Load, instance=Part-1-1
	21,
*End Assembly

由于部件只有一个，因此装配环节较为简单，在这里定义了两个节点集Nset，分别是支座固定位置1和施加荷载位置21。通过定义节点集有利于后续的荷载、边界条件的施加以及最后计算结果的导出。

分段介绍4：Step-Boundary-Output

** STEP:step-1
**
*step, name=step-1, nlgeom=NO, inc=1000
*static
	0.01,	1.,	1e-5,	0.01
**
** BOUNDARY CONDITION
**
*Boundary
Set-Fix, 1, 2
Set-Fix, 6, 6
*Boundary
Set-Load, 1, 1, 0.1
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0 
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
U,
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, frequency=1
*node output, nset=Set-Load
rf1, rf2
*node output, nset=Set-Load
u1, u2
*End Step

最后一部分的代码较为简单，与abaqus可视化操作一致，定义分析步，边界条件，荷载情况和场变量的输出。结合下图应该能够理解每句话的意思，便不再赘述。

建模效果图

一端固定一端施加水平位移
在这里插入图片描述

UEL 的实现

UEL原理和功能简述

在讲UEL代码前先简单阐述一下UEL的原理和功能
简单的说使用UEL就是依据上一迭代步的结果和当前的单元节点位移，通过使用UEL来计算右手残差向量RHS和单元刚度矩阵AMATRX

单元 $[B]$ 矩阵将轴向应变和曲率与单元的的位移 $u_e]$ 关联
材料本构关系矩阵 $[D]$ 将轴向力和力矩与轴向应变和曲率相关联
单元刚度矩阵可由 $[K_e]=\int_0^l[B]^T[D][B]dl$ 求得
实际计算采用数值积分的方式实现

UEL部分代码

通过上述步骤就可以完成 RHS 和 AMATRX 的计算
当然在实际编写UEL过程复杂程度远超上述的几句话，尤其是对于复杂材料的本构实现，下面介绍本实例中采用的平面线弹性梁单元，代码如下：（出于对参考书的知识尊重，没包含ugenb）

      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), svars(*), props(*),
     1 energy(8), coords(mcrd,nnode), u(ndofel), du(mlvarx, *), v(ndofel),
     2 a(ndofel), time(2), params(*), jdltyp(mdload, *), adlmag(mdload, *),
     3 ddlmag(mdload, *),predef(2, npredf, nnode), lflags(*), jprops(*)

      dimension b(2,7), gauss(2)
      parameter(zero=0.d0, one=1.d0, two=2.d0, three=3.d0, four=4.d0,
     1          six=6.d0, eight=8.d0, twelve=12.d0)
      data gauss/.211324865d0,    .788675135d0/
c
c     calculate length and direction cosines
c
      dx = coords(1, 2) - coords(1, 1)
      dy = coords(2, 2) - coords(2, 1)
      dl2 = dx**2 + dy**2
      dl = sqrt(dl2)
      hdl = dl/two
      acos = dx/dl
      asin = dy/dl
c
c     initialize rhs and amatrx
c
      do k1 = 1, 7
          rhs(k1, 1) = zero
          do k2 = 1, 7
              amatrx(k1, k2) = zero
          end do
      end do
c
c nsvars: User-defined number of solution-dependent state variables associated with the element, with a value of 8
c nsvint: For each integral point has four state variables
c
      nsvint=nsvars/2
c
c loop over integral point
c
      do kintk = 1, 2
          g=gauss(kintk)
c
c make b-matrix, in global coordinates
c
          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
c
c calculate increment/total strains and curvatures
c eps: axial total strain, deps: incremental axial strain 
c cap: curvature, dcap: incremental curvature
c
          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
c
c call constitutive routine ugenb
c
          isvint = 1+(kintk-1)*nsvint
          bn = zero
          bm = zero
          daxial = zero
          dbend = zero
          dcouple = zero
          call ugenb(bn, bm, daxial, dbend, dcoupl, eps, deps, cap, dcap,
     1               svars(isvint), nsvint, props, nprops)
c
c assemble rhs ans amatrx
c
          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
c
      return
      end
            subroutine ugenb(bn,bm,daxial,dbend,dcoupl,eps,deps,cap,dcap,
     1                 svint,nsvint,props,nprops)
c
      include 'aba_param.inc'
c
      parameter(zero=0.d0,twelve=12.d0)
c
      dimension svint(*),props(*)
c

到这篇幅似乎有些太长了，我放到下一篇文章介绍UEL的内部实现，以上。

结语

课题进展缓慢，模拟也跑的不好，写博客摸摸🐟
最近疫情好严重，希望大家都健健康康的…

2048 AI社区

有“AI”的1024 = 2048，欢迎大家加入2048 AI社区

更多推荐

LLM开发必备：2025年主流推理框架深度解析与智能体构建全攻略(建议收藏)

2048 AI社区

39_统一多模态理解与生成：2025年的技术融合与创新_深度解析

2048 AI社区

别手动搭n8n了，n8n邪修出手，1个mcp工具让搭建部署全自动

这两天在研究n8n,我发现，n8n的数据，包括工作流的配置，全部都是json格式的，这意味什么？过程非常顺利哈，就是第一遍它没有自动帮我把工作流创建到我的n8n.然后用提示词强调了一下，接着它就调用工具把工作流部署到我的n8n了，看到这里生成了工作流ID说明成功了。因为我抓的网址是动态加载数据的，这个爬取还是有一定难度的，要使用模拟浏览器什么的。不过，我只测试了一两遍，不够说明哈，各位道友，有兴趣