ABAQUS 静/动力隐式UEL 二次开发

概述：结合HHT时程积分法，推导了ABAQUS 静/动力隐式算法中的UEL关键矩阵表达式，并将公式应用到自编CPS4/CPE4、C3D8 BBAR和C3D20用户自定义单元中，计算结果均与ABAQUS自带单元保持一致。其中，静力计算中，关键矩阵AMARTX和RHS等可直接按照刚度矩阵和方程右端不平衡力输出。动力隐式计算中，这两者的输出较为复杂，需要结合HHT时程积分法进行推导，将刚度矩阵、质量矩阵和阻尼矩阵依据LFLAGS数组的数值进行组合，RHS同样需要进行推导计算，并以合适的方法将解相关的状态量储存在SVARS中，供后面的增量步调用。

关于UEL的程序设计，只支持静力通用计算分析步的资料有很多，本帖子内容适不但适用于静力通用，而且适用于动力隐式、频率分析等分析步计算，采用模块化程序设计，所有的矩阵求解均被封装，调用方便，读者可以针对自己的需求对相应函数进行改编即可，尤其是涉及到动力隐式计算部分，适用于任何运动方程的动力隐式求解，可直接移植使用。

----------------------------

（）UEL接口基本参数介绍

----------------------------

用户自定义单元(USER DEFINED ELEEMNT, UEL)适合进阶的工程师/学者使用，u、UEL的目的是实现一个单元的力学行为，即力-位移关系，这部分的功能需要的编程工作量和理论功底都比较高， 涉及到的知识大致包括：

有限元理论：形函数插值、应力-应变关系，刚度矩阵组装，质量矩阵求解，阻尼矩阵求解，

数值算法：高斯积分（全积分、缩减积分）、非线性方程组迭代求解流程（增量迭代，常刚度迭代等）、雅各比矩阵求解

弹/塑性力学：应力-应变关系，位移-应变关系

熟练的FORTRAN编程操作是基础技能。

UEL接口参数众多，包括自己本身的参数和与其他子程序联合的参数，下面介绍我比较重要的几个参数。

<>需要编程定义的参数

RHS（right hand side）：这个命名是从方程组的角度来的，顾名思义，他就是方程的右端量，其本质是：外力-内力，外力部分程序的编写涉及到与其他子程序的联合使用，包括DLOAD和UTRACLOAD等等，这部分内容是给用户自定义单元施加复杂的广义力，UEL接口为其提供了相应的参数，如JDLTYP、NDLOAD等等，目前没有做过尝试。至于内力的求解，在静力线性计算中，数值上等于-KU（后面有这个公式的推导），即刚度与位移乘积负数。在动力隐式分析中，这个矩阵需要结合HHT时程积分法推导具体的表达式。

AMATRX：这个参数往往被认为是刚度，但并不这样，他只有在静力计算的时候才是刚度。他的具体取值依据分析类型的不同而不同，具体的表达形式不唯一。在动力隐式分析中，这个矩阵需要结合HHT时程积分法推导具体的表达式。

SVARS：这个参数官方文档说是取决于结果的状态量，他的具体意义和数据由我们自己确定，而我把它理解为一个小仓库，可以存放我们的数据。需要注意的是里面的数据可以在不同的增量步之间传输，就是说，里面的数据会传到下一个增量步，只要不更新他，他可以一直被使用。这个数据我们自己决定更新与否。

ENERGY：与用户自定义单元相关的一些动量等能量，也是依据不同的分析类型而不同。

 () CPS4/CPE4公式推导

等参单元中，母单元向笛卡尔坐标的转换为：（注意：这里的箭头标反了）

要再现位移场，可假设单元内部位移为坐标的函数，即：

将四个节点的位移和坐标代入上式，有：

如此，求解系数，并用节点位移表示单元内部位移：

上式即为单元节点位移与坐标的插值关系。考虑到等参单元，即位移的插值函数与母单元-笛卡尔单元变换采用相同插值函数。则等参变化的插值函数为：

上式的矩阵表达为：

其中：

对形函数求导，得形函数对母单元坐标的导数：

由弹性力学几何方程，应变矩阵为：

合并为：

JOCABIN矩阵为：

其中：

则JOCABIN的具体表达式为：

刚度矩阵的求解：

质量矩阵分为两种，集中质量矩阵和协调质量矩阵，下面给出的是一致质量矩阵（又叫协调质量矩阵）矩阵，该种质量矩阵因其采用了与位移插值相同的插值函数而被称为一致质量矩阵。但是协调质量矩阵为满阵，不利于数值计算，所以有限元中一般采用集中质量矩阵，即只有对角元素不为零。集中质量矩阵通常由协调质量矩阵处理而来，比如将协调质量矩阵每一行的质量集中到对角元素，然后将非对角元素数值置为零。

阻尼矩阵，采用比例阻尼：

UEL主程序采用常刚度法设计，刚度矩阵和质量矩阵只在第一个增量步计算一次，然后储存在变量SVARS中，后面一直到分析结束之间的所有增量步不再计算，均直接读取SVARS数组中的数据并输出会ABAQUS主程序。该主程序将各个关键量的求解分为不同的子程序，这些子程序基本是所有用户自定义单元都会涉及的，因此该主程序的可移植性很强，在线弹性问题中，即可用常刚度迭代法求解的运动方程中，读者完全可以按照这个直接编写自己的程序。

    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),PROPS(*),
     1 SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL),
     2 DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*),
     3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
     4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)
      double precision DMATX(3,3)
      double precision KK(2*nnode,2*nnode),MM(2*nnode,2*nnode)  
      KK=0.D0   
      MM=0.D0       
!      Print work state on screen and ANA_STAT.txt
      CALL PRINTSTATE(JPROPS,LFLAGS,DTIME,JELEM) 
      IF(KSTEP.EQ.1.AND.KINC.EQ.1) THEN           
!        material martix   
        CALL DMATERIAL(DMatx,PROPS)    
!        stiffness martix     
        call KKmartix(KK,coords,DMatx,mcrd,nnode,jelem)  
!        mass martix  
        call MMmartix(MM,coords,props,mcrd,nnode,jelem,dmatx)
!        Store pre-calculated stiff/mass matrices
        CALL STOREMATRICES (SVARS,NSVARS,KK,MM,NDOFEL)
      else
        CALL READMATRICES(SVARS,NSVARS,KK,MM,NDOFEL)
      endif 
!      Output required variables
      CALL OUTPUTVARIABLE(RHS,AMATRX,SVARS,PROPS,ENERGY,U,V,A,
     1 LFLAGS,DTIME,NDOFEL,NRHS,NSVARS,MLVARX,JELEM,PARAMS,KK,MM)
      RETURN
      END

不想写文字了，麻烦，直接把付费内容放在附件里面了：

-------------------------------------------------

-------------------------------------------------

（1）CPE4/CPS4单元UEL，适用于静/动力隐式计算

（2）C3D8 BBAR单元UEL，，适用于静/动力隐式计算，包含UMAT给应力应变可视化部分

（3）C3D20 单元UEL，适用于静/动力隐式计算

（4）WORD文档，推导、讲解动力隐式计算中的AMARTX、RHS

赠送内容：

（5）收集的一些UEL-UMAT程序

（6）一对一辅导，包括理论、公式和编程等......

-------------------------------------------------

-------------------------------------------------