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

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

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

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

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

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

分段介绍1：Heading

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

分段介绍2：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

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

分段介绍4：Step-Boundary-Output

最后一部分的代码较为简单，与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/cc     calculate length and direction cosinesc 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/dlcc     initialize rhs and amatrxc do k1 = 1, 7  rhs(k1, 1) = zero  do k2 = 1, 7      amatrx(k1, k2) = zero  end do end docc nsvars: User-defined number of solution-dependent state variables associated with the element, with a value of 8c nsvint: For each integral point has four state variablesc nsvint=nsvars/2cc loop over integral pointc do kintk = 1, 2  g=gauss(kintk)cc make b-matrix, in global coordinatesc  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) = zerocc calculate increment/total strains and curvaturesc eps: axial total strain, deps: incremental axial strain c cap: curvature, dcap: incremental curvaturec  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 docc call constitutive routine ugenbc  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,1svars(isvint), nsvint, props, nprops)cc assemble rhs ans amatrxc  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))

                                                        声明：本站部分文章及图片源自用户投稿，如本站任何资料有侵权请您尽早请联系jinwei@zod.com.cn进行处理,非常感谢！