CPE8/CPS8用户自定义单元二次开发（1）_Abaqus_二次开发

CPE8/CPS8用户自定义单元二次开发（1）

零

概述

有限元法编写二维四边形二阶用户自定义单元，采用abaqus提供的uel二次开发平台，包括单元的刚度矩阵K、质量矩阵M和阻尼矩阵C，uel需要向abaqus主程序输出AMARTIX、RHS等矩阵。

编写的uel子程序适用于静动力计算、频域计算等等，帖子分为两个部分

->理论部分CPE8/CPS8 UEL用户自定义单元开发（1）结果部分CPE8/CPS8 UEL用户自定义单元开发（2）

壹

形函数/位移插值函数

在经典的三角形单元中，位移的插值是直接在子单元中进行，但是随着单元形状的复杂，直接在子单元进行计算回带来计算的困难的公式推导的困难，因此等参单元法应运而生。在等参元中，单元形状映射中的参数和位移场插值的参数数目相等，因此得名等参元。

采用等参单元法计算二维八节点单元，其中，母单元为

子单元，即我们绘制的网格真实形状为

假设的坐标映射函数为

将母单元中八个节点的坐标导入上述式子，求解出16个系数a0~a15，然后求得八节点等参元的形函数

这是形状映射中x坐标的映射，对于y方向得映射，将上述所有x替换为y即可。

写成求和形式

写成矩阵形式为

其中，对于边角点有

对于边中点有

上述公式，对于位移插值同样适用。

贰

雅各比矩阵

等参元中的数值积分中涉及到自然坐标和局部坐标的转换，这种转换通过雅各比矩阵实现。

根据链式求导法则有

写成矩阵形式

雅各比矩阵为

将式

代入雅各比矩阵有

形函数对母单元坐标的导数，可以直接求导，因此，等参元的雅各比矩阵可以求得。

则有

因此，雅各比矩阵的逆矩阵也可求得。

叁

B矩阵

有限元中，单元的刚度矩阵为

D矩阵是弹性矩阵，在线弹性问题中，D矩阵为常数矩阵，根据物理方程可以推导出。B矩阵是应变-位移矩阵，用作将位移转化为应变，下面重点求解B矩阵。

应变与位移有如下关系

其中算子L为

单元的位移u为

则应变与位移的关系可表示为

则B矩阵为B=LN，具体表达式为

将B矩阵写成一个矩阵为

其中B_i为

在求B矩阵的时候，我们的目标是求得形函数对xy的倒数，这时候我们已经求得了形函数对局部坐标的倒数，因此，通过雅各比矩阵我们就可有求得形函数对xy的倒数，首先，引入以下矩阵

通过雅各比矩阵，上述矩阵可以表示为

上面等式右面的雅各比逆矩阵，我们已经求得，形函数对st坐标的导数也已经求得，因此，上面引入的矩阵就可以求出来，然后再按照元素对应的规律将上述矩阵以此存放到B矩阵中。

在计算质量矩阵的时候，首先计算单元的额协调质量矩阵，然后将行质量集中到对角元，让后将非对角元素置为零，具体可参考帖子：

质量矩阵的求解
‍
有限元先生，公众号：有限元先生C3D8 BBAR UEL用户自定义单元二次开发（1）

单元的阻尼矩阵可按照比例阻尼求解，即刚度矩阵和质量矩阵的线性组合。

求解刚度矩阵的主程序为

      subroutine KKmartix(KK,coords,DMatx,mcrd,nnode,jelem)          INCLUDE 'ABA_PARAM.INC'      double precision coords(mcrd,nnode),dmatx(3,3)      double precision INTEGERPOINT(3),wgt(3)      double precision KK(2*nnode,2*nnode)      double precision B(3,2*NNODE),J(2,2)      double precision DETJ,G,H,COFF      double precision NPXY(2,NNODE)                 INTEGERPOINT(1)= 0.D0      INTEGERPOINT(2)= 0.7745966692D0      INTEGERPOINT(3)=-0.7745966692D0                  wgt(1)=0.8888888889D0      wgt(2)=0.5555555556D0      wgt(3)=0.5555555556D0                  B=0.d0      J=0.d0      DETJ=0.d0      G=0.d0      H=0.d0      COFF=0.d0      NPXY=0.d0                               do I=1,3        do II=1,3            G=INTEGERPOINT(I)            H  =INTEGERPOINT(II)            call CalBJ(B,J,DETJ,G,H,MCRD,NNODE,JELEM,COORDS)            COFF=DETJ*WGT(I)*WGT(II)                                          KK=KK+MATMUL(MATMUL(TRANSPOSE(B),DMatx),B)*COFF        enddo          enddo               return       end

求解质量矩阵的主程序为

      subroutine MMmartix(MM,COORDS,PROPS,MCRD,NNODE,JELEM)      INCLUDE 'ABA_PARAM.INC'            INTEGER MCRD,NNODE,JELEM      DOUBLE PRECISION MM(2*NNODE,2*NNODE),COORDS(MCRD,NNODE)      DOUBLE PRECISION NN(2,16)      DOUBLE PRECISION INTEGERPOINT(3),WGT(3)      DOUBLE PRECISION G,H,DETJ,DENSITY,COFF      DOUBLE PRECISION J(2,2),PROPS(*)            MM=0.D0      wgt=0.D0      JJ=0.D0        NN=0.D0            G=0.D0      H=0.D0        j=0.D0                 DETJ=0.D0            Coff=0.D0        DENSITY=0.D0                                  INTEGERPOINT(1)= 0.D0      INTEGERPOINT(2)= 0.7745966692D0      INTEGERPOINT(3)=-0.7745966692D0                  wgt(1)=0.8888888889D0      wgt(2)=0.5555555556D0      wgt(3)=0.5555555556D0               DENSITY=PROPS(3)                              DO I=1,3        DO II=1,3          G  =INTEGERPOINT(I)          H  =INTEGERPOINT(II)              CALL CalDETJ(J,DETJ,G,H,MCRD,NNODE,JELEM,COORDS)          CALL CalNN(NN,G,H)                                COFF=DETJ*WGT(I)*WGT(II)*DENSITY                      MM  =MM+(MATMUL(TRANSPOSE(NN),NN))*COFF        ENDDO       ENDDO            return      end

至此，单元的刚度矩阵、质量矩阵和阻尼矩阵已经求得，uel程序的核心部分已经完成，下一个帖子会更新自编CPE8/CPS8 uel程序的计算结果。

来源：有限元先生

abaqus 界面美化科研绘图

零概述有限元结果展示是后处理中的重要步骤，一图胜万言，在与人讨论分析计算结果的时候，一个清晰美观的abaqus可视化界面使工作事半功倍；在正式的学术场合公开自己的研究成果的时候，一个符合学术要求的图片更加重要。因此，无论是展示结果还是公开研究成果，对abaqus的可视化界面进行适当的调整是有必要的，本着这个目的，帖子分享一些具体的调整技巧。帖子以一个涉及到接触非线性的问题为例，讲解一系列的操作。算例为两个含有槽口的模型受动荷载作用下接触的odb文件，单个槽口的尺寸为整体模型的几何图形为荷载和边界条件为其中，S1界面为固定边界，S2为施加面荷载的边界，S3为接触界面。对两个块体采用非匹配网格离散，计算模型为下面以该算例的odb为例子讲解若干个abaqus可视化界面的调整技巧。背景调整 abaqus界面默认的是灰色渐变颜色，没有调整之前为如果需要调整背景颜色，在上部菜单栏依次点击1-Views，选择下拉菜单中2-Graphics Options，弹出3-graphics options 可以看到，abaqus中默认的颜色是4-gradient，即灰色渐变。我们将颜色调整为纯白色，需要点击上面图片里的5-solid，弹出界面6-select color 选7中的白色，点击ok即将背景调整为白色。如图界面文字调整显示该算例的应力云图（显示什么不重要，只为了展示）为首先说1部分的文字内容修改，我们将1部分修改为罗马字体，并且修改字体的大小。在上部菜单栏点击5-viewport，选择下拉框的6-viewport...，弹出7-viewport...，再选择8-legend，调整下面的内容即可，包括字体种类、字体大小等等。这里展示了调整后的效果。将字体修改为罗马字体，修改了字的大小，隐藏了原有的框框，加粗了字体。效果还是可以的。对于云图下方的一些文字，有一些是不需要的，如下图红线表示了对应的关系，上图是显示了数据，如果将不勾选相应的框框，下面的内容也会隐藏。如图迭代信息是需要查看的，因此上图保留了迭代计算信息。至于左下角和右上角的坐标轴，也是在viewport...界面里面显示和隐藏的。坐标轴的xyz字体也是可以调节的，如图把坐标轴的xyz调整的很大，修改了颜色且设置了粗斜体。变形因子下图显示的是未经调整的云图（文字经过了调整）这部分主要解决变形因子的调整。变形因子是abaqus缩放真实变形的一个参数，有时候模型的变形非常的大或小，就得对模型进行缩放以进行更好的观察模型的变形情况。在odb窗口的左侧选择1-common plot options，弹出2-common plot options窗口，可以看到3-auto-compute是默认的选择，abaqus给我们的模型设置的默认变形因子为249410（这个数字随模型不同而不同）。我们选择4-unifom，在下面自己设置变形因子，这里将变形因子设置为1，即真实变形，设置后的变形为云图连续显示 abaqus默认的云图不是连续显示，相当于在等值线的基础上给不同的区域设置了颜色，不同数值的界面颜色不是渐变的，如在odb界面左侧点击1-contour options，弹出contour plot options窗口，可以看到abaqus的默认设置为2-discrete，云图中的数值3显示的是不连续的云图，默认设置不是很美观，abaqus提供了选项设置渐变云图，如在contour plot options中选择4-continue，再查看数字5指示的云图，已经变为连续，相较美观一些。隐藏网格 abaqus默认的是显示网格，如在odb界面的左侧选择1-common options，弹出窗口2-common plot options，abaqus默认的选择是3-all edges，下面也给出了很多的选项，我们点击下面的feature edges选项，即只显示模型的结构线，隐藏所有内部的网格线，修改后的效果为多窗口显示上面显示的都是一个odb模型，当需要多个odb模型对比数据的时候，abaqus提供了多个odb模型同时显示的功能，如下图点击上方的1-create viewport，创建一个窗口，然后点击2-tile viewport vertically将所有的窗口按照平行排列，这样就能比较不同odb的计算结果。也可以显示多个odb，如总结帖子分享了一些调整abaqus可视化界面的一些小技巧，包括调整背景颜色、文字方面的调整、网格的显示与隐藏、云图的连续设置以及变形因子调整等等，类似的操作还有很多，其中每一个二级菜单里面都有很多的功能，有待探索，只要花时间调整，abaqus基本能满足大部分用户的需求。点击卡片关注我们来源：有限元先生