Abaqus子程序之UMAT

0. 前言

这次主要是简单介绍下umat子程序的帮助文档介绍和简单方法，最后的提供的案例也是我学会的第一个umat子程序，较为清晰的展示了umat编写的基本流程：增量形式本构公式推导，雅可比矩阵的推导，将前两步的成果写入umat。实际上umat就是实现更新应力，核心就是增量本构的推导。

1.简介

 可用于定义材料的力学本构行为；
建议在具有规定牵引（简单拉伸）载荷的单个元素模型上进行初始试验。
 将在材料定义包括用户定义的材料行为的元素的所有材料计算点调用；
可用于包括力学行为的任何程序；
 可以使用依赖于求解的状态变量；
 对于力学本构模型，必须提供材料雅可比矩阵，C=(1/J)∂Δ(Jσ)/∂Δε
 可以与用户子程序USDFLD一起使用，以便于在任何场变量传入（uamt）之前重新定义它们；

2. 子程序界面设置
​

3. 用户子程序接口

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
     4 JSTEP(4)
 
      user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
      and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
 
      RETURN
      END

4. 要定义的变量

（1）所有情况都要定义的变量
DDSDDE(NTENS,NTENS)：
本构模型的雅可比矩阵，时基尔霍夫(Kirchhoff )应力增量；表征参考构型体积变化的变形梯度的行列式；是柯西应力；是应变增量；如果体积变化很小，雅可比矩阵可以近似为∂Δσ/∂Δε，Δσ是柯西应力增量。DDSDDE（I，J）定义了由应变增量阵列的第J个分量的无穷小扰动引起的时间增量结束时第I个应力分量的变化。除非为用户定义的材料调用非对称方程求解功能，否则Abaqus/Standard将仅使用DDSDDE的对称部分。矩阵的对称部分是通过取矩阵及其转置之和的一半来计算的。
对于频域中的粘弹性行为，雅可比矩阵的尺寸必须为DDSDDE（NTENS，NTENS，2）。刚度贡献（储能模量）必须在DDSDDE（NTENS，NTENS，1）中提供，而阻尼贡献（损耗模量）必须提供在DDSDD（NTENS、NTENS，2）中。
STRESS(NTENS)：
该数组在增量开始时作为应力张量传入，并且必须在此程序中更新为增量结束时的应力张量。如果指定了初始应力（初始条件），则该数组将包含分析开始时的初始应力。此数组的大小取决于下面定义的NTENS值。在有限应变问题中，在调用UMAT之前，应力张量已经旋转，以说明刚体在增量中的运动，因此在UMAT中只应对应力的共旋转部分的进行积分。所用的应力度量是“真实”（柯西）应力。如果UMAT使用了一种总体的混合公式（与默认的增量行为相反），则应力数组将扩展到NTENS之外。数组的第一个NTENS条目包含所有应力（数组的前NTENS个元素都是应力），如上所述（8中不可压缩材料）。额外数量如下：
STRESS(NTENS+1): ，该值只可读取（在程序中只能从外界读取，不能修改和定义），
STRESS(NTENS+2): 该值只可写入（在程序中只能定义和修改，不能从外界获取），
STRESS(NTENS+3): ，该值只可写入（在程序中只能定义和修改，不能从外界获取），这里U是应变能密度势的体积部分（见8）。
STATEV(NSTATV)：
包含求解相关的状态变量（sdv）的数组。这些值将作为增量开始时的值传入，除非在用户子程序USDFLD或UEXPAN中更新，在这种情况下，将传入更新的值。在所有情况下，STATEV都必须返回增量结束的值（以便于下一增量初始时使用或者其他子程序使用）。数组的大小定义如“ Allocating Space for Solution-Dependent State Variables”中所述。
在有限应变问题中，除了更新与本构行为相关的所有值之外，任何向量值或张量值的状态变量都必须旋转以说明材料的刚体运动。为此，提供了旋转增量矩阵DROT。
SSE, SPD, SCD：
分别为弹性应变能、塑性耗散（能）和“蠕变”耗散（能）。这些值作为增量开始时的值传入，并应在增量结束时更新为相应的特定能量值。除了用于能量输出外，它们对求解没有影响。
（2）仅在完全热—应力耦合或者完全热—电—结构耦合分析中定义的变量
RPL：
由材料力学加工导致的，在增量结束时单位时间内产生的体积热量。
DDSDDT(NTENS)：
与温度相关的的应力增量的变化。
DRPLDE(NTENS)：
RPL相对于应变增量的变化。
DRPLDT：
RPL相对于温度增量的变化。
（3）仅在土压应力程序或孔隙压力粘性单元(cohesive)的孔隙流体扩散/应力耦合分析中
RPL：
RPL用于指示粘性单元(cohesive)是否对孔隙流体的切向流动开放。如果没有切向流，将RPL设置为0；否则，如果单元开放，则为RPL分配一个非零值。一旦开放，粘性单元(cohesive)将保持对流体流动开放。

5. 可更新的变量

PNEWDT：
建议的新时间增量与使用的时间增量的比率（DTIME，即使用的时间增量，见本节后面的讨论）。此变量允许您为Abaqus/Standard中的自动时间增量算法提供输入（如果选择了自动时间增量）；对于准静态程序，Abaqus/Standard使用的自动时间步是基于积分标准蠕变定律的技术（参见 Quasi-Static Analysis），不能在UMAT子程序内进行控制。
在每次调用UMAT之前，将PNEWDT设置为大值。
如果PNEWDT被重新定义为小于1.0，则Abaqus/Standard必须放弃时间增量（使用的时间增量），并以较小的时间增量再次尝试。为自动时间积分算法提供的建议新时间增量为PNEWDT×DTIME，其中使用的PNEWDT是允许该迭代重新定义PNEWDT的、所有调用的用户子程序的最小值。
如果对于本次迭代的所有用户子程序调用，PNEWDT的值大于1.0，且增量在本次迭代中收敛，则Abaqus/Standard可增加时间增量。为自动时间积分算法提供的建议新时间增量为PNEWDT×DTIME，其中所使用的PNEWDT是该迭代中所有调用用户子程序的最小值。
如果在分析过程中未选择自动时间增量，则大于1.0的PNEWDT值将被忽略，小于1.0的PNEWDT值将导致作业终止。

6. 传递信息的变量

STRAN(NTENS)：
包含增量开始时的总应变的数组。如果热膨胀包含在同一材料定义中，则传递到UMAT中的应变仅为力学应变（即，基于热膨胀系数计算的热应变已从总应变中减去）。这些应变可作为“弹性”应变输出。
在有限应变问题中，在调用UMAT之前，应变分量已被旋转以考虑增量中的刚体运动，并且是对数应变的近似值。
DSTRAN(NTENS)：
应变增量数组。如果热膨胀包含在同一材料定义中，则这些是力学应变增量（总应变增量减去热应变增量）。
TIME(1)：当前增量或频率开始时的分析步时间值。
TIME(2)：当前增量开始时的总时间值。
DTIME：时间增量。
TEMP：当前增量开始时的温度。
DTEMP：温度增量。
PREDEF：
基于在节点处读入的值，在增量开始时的这一点上预定义场变量的插值数组。
CMNAME：
用户定义的材料名称，左对齐。一些内部材料模型的名称以“ABQ_”字符串开头。为避免冲突，不应使用“ABQ_”作为CMNAME的前缀字符串。
NDI：一点处直接应力分量的数量。
NSHR：一点处工程剪切应力分量的数量。
NTENS：应力或应变分量数组的尺寸（NDI+NSHR）。
NSTATV：
与此材料类型关联的求解相关状态变量的数量（如“Allocating Space for Solution-Dependent State Variables”中所述定义）。
PROPS(NPROPS)：与此用户材料关联的用户指定的材料常数数组。
NPROPS：与此用户材料关联的用户定义的材料常数个数。
COORDS（3）：
包含此点坐标的数组。如果在步骤中考虑几何非线性，则这些是当前坐标（请参见定义分析）；否则，数组包含点的原始坐标。如果是二维模型，则数组为COORDS（2）。
DROT(3,3):
旋转增量矩阵。该矩阵表示存储应力（应力）和应变（STRAN）分量的基础系统（基本坐标系）的刚体旋转增量。提供这种方法是为了在该子程序中适当地旋转矢量或张量值状态变量：在调用UMAT之前，应力和应变分量已经被这个量旋转了。如果材料点的基本系统（基本坐标系）随材料旋转（如在壳单元中或使用局部定向时），则该矩阵在小位移分析和大位移分析中作为单位矩阵传递。
CELENT:
单元特征长度，这是一阶单元中穿过单元的线的典型长度；它是二阶单元相同典型长度的一半。对于梁和桁架，它是沿单元轴线的特征长度。对于膜和壳，它是参考表面中的特征长度。对于轴对称单元，它仅是平面(r,z)的特征长度。对于内聚力单元（cohesive），它等于基本单元的厚度。
DFGRD0(3,3):
数组，包含增量开始时的变形梯度。如果在材料点定义了局部方向，在增量开始时，变形梯度分量将在由（局部）方向定义的局部坐标系中表示。有关各种单元类型的变形梯度可用性的讨论，请参见变形梯度。
DFGRD1(3,3):
数组，包含增量结束时的变形梯度。如果在材料点定义了局部方向，则变形梯度分量将在由（局部）方向定义的局部坐标系中表示。如果与此增量相关的步长定义中不包括非线性几何效应，则将此阵列设置为单位矩阵。有关各种图元类型的变形梯度可用性的讨论，请参见变形梯度。
NOEL : 单元编号。
NPT : 积分点编号。
LAYER : 层号（用于复合壳和分层实体）。
KSPT : 当前层内的截面点编号。
JSTEP(1) ：分析步编号。

7. 示例：使用多个用户定义的力学材料模型

要使用多个用户定义的材料模型，变量CMNAME可以在用户子程序UMAT中测试不同材料名称，如下所示：

IF (CMNAME(1:4) .EQ. 'MAT1') THEN
CALL UMAT_MAT1(argument_list)
ELSE IF(CMNAME(1:4) .EQ. 'MAT2') THEN
CALL UMAT_MAT2(argument_list)
END IF

8. 示例：简单线性粘弹性材料

这是用户子程序UMAT编码的一个简单示例，请考虑图1中所示的线性粘弹性模型。它可以说明如何对umat子程序进行编码的完整过程。

 

▲ 图1 简单的线性粘弹性模型

这是用户子程序UMAT编码的一个简单示例，请考虑图1中所示的线性粘弹性模型。它可以说明如何对umat子程序进行编码的完整过程。

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),
2 DDSDDT(NTENS),DRPLDE(NTENS),
3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
DIMENSION DSTRES(6),D(3,3)
C
C EVALUATE NEW STRESS TENSOR
C
EV = 0.
DEV = 0.
DO K1=1,NDI
EV = EV + STRAN(K1)     !! 直接应变和
DEV = DEV + DSTRAN(K1)  !!直接增量应变和
END DO
C
TERM1 = .5*DTIME + PROPS(5) !!
TERM1I = 1./TERM1            !!
TERM2 = (.5*DTIME*PROPS(1)+PROPS(3))*TERM1I*DEV   !!
TERM3 = (DTIME*PROPS(2)+2.*PROPS(4))*TERM1I  !!
C 计算正应力增量和更新正应力
C 正应力增量用下面的本构方程
 
C     ++
DO K1=1,NDI
DSTRES(K1) = TERM2+TERM3*DSTRAN(K1)
1     +DTIME*TERM1I*(PROPS(1)*EV
2+2.*PROPS(2)*STRAN(K1)-STRESS(K1))
C        
STRESS(K1) = STRESS(K1) + DSTRES(K1)
END DO
C  更新剪应力
C   TERM2=
C 本构方程
C  
C *+**
TERM2 = (.5*DTIME*PROPS(2) + PROPS(4))*TERM1I
I1 = NDI
DO K1=1, NSHR
I1 = I1+1
DSTRES(I1) = TERM2*DSTRAN(I1)+
1     DTIME*TERM1I*(PROPS(2)*STRAN(I1)-STRESS(I1))
STRESS(I1) = STRESS(I1)+DSTRES(I1)
END DO
C 计算雅可比距阵元素并创建雅可比矩阵
C  TERM2=*
TERM2 = (DTIME*(.5*PROPS(1)+PROPS(2))+PROPS(3)+
1  2.*PROPS(4))*TERM1I
​
TERM3 =*
TERM3 = (.5*DTIME*PROPS(1)+PROPS(3))*TERM1I
DO K1=1,NTENS
DO K2=1,NTENS
DDSDDE(K2,K1) = 0.
END DO
END DO
C
DO K1=1,NDI
DDSDDE(K1,K1) = TERM2
END DO
C
填充距阵如下位置
​
DO K1=2,NDI
N2 = K1–1
DO K2=1,N2
DDSDDE(K2,K1) = TERM3
DDSDDE(K1,K2) = TERM3
END DO
END DO
TERM2 =*
TERM2 = (.5*DTIME*PROPS(2)+PROPS(4))*TERM1I
I1 = NDI !!l1=3
C  按新TERM2，填充距阵如下位置（实际上没必要这么麻烦。因为雅可比矩阵的每一个元素都有确定的表达式，参见非线性本构关系笔记P38）
​
DO K1=1,NSHR       !!!按剪切力运行3次
I1 = I1+1
DDSDDE(I1,I1) = TERM2
END DO
C     以上步骤，雅可比矩阵填充完毕
C 该增量步的总比能的改变量
C TOTAL CHANGE IN SPECIFIC ENERGY
C    当前增量步总比能量改变计算公式：
C 上面的表达式是张量表达式，实际上是六个应力、应变分量对应的比能量之和：
TDE = 0.
DO K1=1,NTENS
TDE = TDE + (STRESS(K1)+.5*DSTRES(K1))*DSTRAN(K1)
END DO
C 
C     该增量步的比弹性能的改变量
C  D矩阵是四阶弹性刚度矩阵，下面并没有直接创建D矩阵，而是根据3个正应力与三个正应变的关系，得到简化的正应变与正应力之间满足的刚度矩阵，本质上就是公式：，将正应力—正应变和剪应力—剪应变的关系分开写，可以这样写是因为正应力与剪应变不耦合，剪应力与正应变不耦合。
C CHANGE IN SPECIFIC ELASTIC STRAIN ENERGY
 
TERM1 = PROPS(1) + 2.*PROPS(2) !!TERM1=
C   填充矩阵如下位置
​
DO  K1=1,NDI
D(K1,K1) = TERM1
END DO
C  填充距阵如下位置
​
DO K1=2,NDI
N2 = K1-1
DO K2=1,N2
D(K1,K2) = PROPS(1)
D(K2,K1) = PROPS(1)
END DO
END DO
​
C 分析以下循环，实质为能量变化，力*位移
DEE = 0.
DO K1=1,NDI
TERM1 = 0.
TERM2 = 0.
DO K2=1,NDI
TERM1 = TERM1 + D(K1,K2)*STRAN(K2)
TERM2 = TERM2 + D(K1,K2)*DSTRAN(K2)
END DO
DEE = DEE + (TERM1+.5*TERM2)*DSTRAN(K1)
END DO
​
C 剪切应力部分
I1 = NDI
DO K1=1,NSHR
I1 = I1+1
DEE = DEE + PROPS(2)*(STRAN(I1)+.5*DSTRAN(I1))*DSTRAN(I1)
END DO
C 最终的DDE为直接应力与剪切应力发生能量改变之和
最终能量消散所改变的值SSE,SCD
SSE = SSE + DEE         !! 总弹性比能。右边SSE是更新前（或上一增量步）的总弹性比能，它是自动传入子程序的值；
!!这里是更新前（或上一增量步）的总比弹性能加上当前比弹性能增量。
SCD = SCD + TDE – DEE   !!更新前（或上一增量步）的蠕变消散加上“总比能增量与弹性比能增量的差值，即蠕变耗散增量，这里未考虑塑性”。
RETURN
END

结语

       案例主要是展示了一个清晰完整的umat编写过程，是入门的一个非常好的案例，主要是可以用来参考学习，我分享了一些其他案例供大家学习，公众号冬生亦东生回复umat子程序即可获取，我本人也学习子程序时间不久，可能会有不当之处，欢迎大家一起学习交流。