代码编织梦想

续上文
下面还是采用分段的方法对上述UEL进行讲解

分段介绍1:UEL接口

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(*)

第一部分为UEL的接口,不可修改,主要是传入相关参数让我们去计算RHS和AMATRX
相关用到的参数会在后面讲解中说明,其余参数可以参考ABAQUS的帮助文档

分段介绍2:常量和变量声明定义

      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/

第二部分是开辟存储空间和声明所要用到的常量,包括 [ B ] [B] [B]矩阵,gauss(2)积分点位置

分段介绍3:单元长度和方向向量计算

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

这一部分就是计算单元的长度和方向向量,很好理解

分段介绍4:初始化 rhs 和 amatrx

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

初值均为零,外部荷载情况由abaqus考虑我们仅需要提供单元内部应变完全产生的内力。对于变量赋予初值是良好的编程习惯。

分段介绍5: [ B ] [B] [B]矩阵

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

nsvars是abaqus传入的参数,是我们在inp文件中*user element 中设置的变量variable=8
设置为8的原因后面说明。
接下来就是循环的开始,遍历每一个积分点do kintk,在上篇中文章的开头已经说明了该单元为采用两积分点的单元。采用数值积分的方法将积分转为求和。再具体的讲就是:单元纵向力 F ,转为了两积分点处 g1 和 g2 的力的和
[ F e ] = ∫ 0 l [ B ] T { F M } d l = ∑ i = 1 n W l ∗ [ B ] T { F M } [F_e]=\int_0^l[B]^T \begin{Bmatrix} F\\M \end{Bmatrix}\quad dl=\sum_{i=1}^n W_l*[B]^T \begin{Bmatrix} F\\M \end{Bmatrix}\quad [Fe]=0l[B]T{FM}dl=i=1nWl[B]T{FM}
接下来讲解 [ B ] [B] [B]矩阵,前面一直没提到的一个重要的问题就是单元的方向(局部坐标系)与整体的方向(坐标系)不一致的问题,也就是坐标转换问题,这里给出的 [ B ] [B] [B]矩阵实际为 [ B e ] ∗ [ λ ] [B_e]*[\lambda] [Be][λ]起到了将单元在整体坐标系下的位移转换为在局部坐标系下的位移的作用。具体的细节不在本文讲解的范围内,相信各位进行UEL编写时候一定对有限单元法有了一定的基础。
这里仅给出局部坐标下的 [ B e ] [B_e] [Be]和转换矩阵 [ λ ] [\lambda] [λ]以供参考,其中位移函数的插值为,纵向二次,法向三次,注意节点顺序:左节点,右节点,中间节点
[ B e ] = [ ( − 3 + 4 ξ ) 0 0 ( − 1 + 4 ξ ) 0 0 ( 4 − 8 ξ ) 0 ( − 6 + 12 ξ ) ( − 4 + 6 ξ ) 0 ( 6 − 12 ξ ) − 2 + 6 ξ 0 ] [B_e] = \begin{bmatrix} (-3+4\xi) & 0 & 0 & (-1+4\xi) & 0 & 0 & (4-8\xi) \\0 &(-6+12\xi) & (-4+6\xi) & 0 & (6-12\xi) & {-2+6\xi} & 0\end{bmatrix}\quad [Be]=[(3+4ξ)00(6+12ξ)0(4+6ξ)(1+4ξ)00(612ξ)02+6ξ(48ξ)0]
[ λ ] = [ c o s α s i n α 0 0 0 0 0 − s i n α c o s α 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 c o s α s i n α 0 0 0 0 0 − s i n α c o s α 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 ] [\lambda] = \begin{bmatrix} cos\alpha & sin\alpha & 0 & 0 & 0 & 0 & 0 \\-sin\alpha &cos\alpha & 0 & 0 & 0 & 0 & 0 \\ 0 & 0& 1& 0& 0& 0& 0 \\ 0&0&0& cos\alpha & sin\alpha & 0& 0 \\ 0& 0& 0& -sin\alpha &cos\alpha & 0 & 0 \\ 0 &0&0& 0&0&1&0 \\ 0&0&0&0&0&0&1 \end{bmatrix}\quad [λ]=cosαsinα00000sinαcosα000000010000000cosαsinα00000sinαcosα0000000100000001

分段介绍6:单元内力计算

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

eps为应变 ε \varepsilon ε 、cap为曲率 κ \kappa κ在初始化后,通过循环的方式,实际是进行了矩阵运算 { ε κ } = [ B ] [ u e ] \begin{Bmatrix} \varepsilon\\ \kappa \end{Bmatrix} = [B][u_e] {εκ}=[B][ue]
u(*)为abaqus传入的节点位移,按照单元内定义节点的顺序。

这里为了程序清晰简洁将本构关系提取为一个子程序ugenb,在调用本构之前对轴力bn、弯矩bm、当前切线轴向刚度和弯曲刚度daxial、dbend以及切线刚度耦合项都赋予初值0,之后在本构中实现了:
{ F M } = [ D ] { ε κ } \begin{Bmatrix} F\\M \end{Bmatrix} = [D]\begin{Bmatrix} \varepsilon \\ \kappa \end{Bmatrix} {FM}=[D]{εκ}
并将计算的结果返回。

前面提到inp文件中*user element 中设置的变量variable=8
是因为在本构中svars需要一定的空间存储变量,该变量不会随着uel调用结束释放,这里每个积分点需要4个变量因此设为8。

分段介绍7:rhs 和 amatrx 的集成

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

单元内力集成其实就是就是通过前述的公式实现的
[ F e ] = ∫ 0 l [ B ] T { F M } d l = ∑ i = 1 n W l ∗ [ B ] T { F M } [F_e]=\int_0^l[B]^T \begin{Bmatrix} F\\M \end{Bmatrix}\quad dl=\sum_{i=1}^n W_l*[B]^T \begin{Bmatrix} F\\M \end{Bmatrix}\quad [Fe]=0l[B]T{FM}dl=i=1nWl[B]T{FM}
这里的负号是由于以外力为正,因此内力即为负
刚度矩阵的集成则依据公式:
[ K e ] = ∫ 0 l [ B ] T [ D ] [ B ] d l [K_e]=\int_0^l[B]^T[D][B]dl [Ke]=0l[B]T[D][B]dl
只是将积分用数值积分的形式求解,其中的 [ B ] T [ D ] [B]^T[D] [B]T[D]通过bd1、bd2计算得到。

小结

以上就是整个UEL的实现了,我自己应该是理清楚了,但是不知道有没有表达清楚。
不过想要自己完完全全的完成一个UEL的编写还是非常复杂的过程,需要对有限单元法有一个深入的学习,相关书籍我使用过的是 王勖成老师的《有限单元法》 书挺厚但是讲解的非常细致推荐给大家阅读。
目前我还没有深入了解UEL的调试方法,希望后续能够通过阅读书籍和博客掌握相关方法,然后分享给大家。

摸🐟摸🐟

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/weixin_39310341/article/details/123552038

abaqus二次开发概述-爱代码爱编程

随着计算技术和计算机的快速发展,有限元软件的发展速度迅速,功能日渐强大。目前国际上被广泛采用的通用有限元软件有 ANSYS、MSC、ABAQUS 等。利用商业软件进行计算现在已是科学研究中的一项重要手段。由于工程问题的千差万别,不同的用户有不同的专业背景和发展方向,通用软件不免在具体的专业方面有所欠缺,针对这些不足,大部分的通用软件都提供了二次开发功能,以

ANSYS Help USER300-爱代码爱编程

ANSYS Help USER300 1 概述1.1 单元描述1.2 输入数据1.3 假设与限制条件2 创建新单元2.1 输入和输出缩略词2.2 通过自定义单元API创建新单元2.2.1 建议和限制2.2.2 子程序 `UserElem`2.2.2 子程序 `ElemGetMat`2.2 直接访问程序数据库创建新单元2.2.1 用户子程序概述2.2

一个简单的UEL梁单元(上)-爱代码爱编程

出于课题需要,学习研究了abaqus uel二次开发,这里通过一个简单的平面线弹性梁单元的UEL实例来复习巩固一下学到的东西,为后续钢筋和混凝土的本构修改做准备。 主要的参考资料为:王涛老师的《基于ABAQUS的有限元子程序开发及应用》和阚前华著的《非线性本构关系在ABAQUS中的实现》 推荐大家看看 图示为实例中使用的单元节点和Gauss积