当前位置：服务支持 > 软件文章 > ABAQUS子程序学习笔记精华（1）

ABAQUS子程序学习笔记精华（1）

阅读数 8

1 Fortran语言

Fortran之名来自于FORmula TRANslation，原意为数学公式编辑器，用于进行科学和工程计算。它被广泛应用于高性能计算、数值计算和科学计算中。它最初于20世纪50年代开发，于1977年发展出广为人知的Fortran 77语言标准。1990年初发展出新一代标准Fortran 90 ，打破了Fortran 77的固定格式(fixed format)，改为现代的自由格式(free format)，之后发展出Fortran 95 、Fortran 2003、Fortran 2008等标准，以符合与满足当代开发者之思维与需求。

Fortran具有许多特点，如动态存储分配、指针、递归和自定义数据类型等，而且它采用的是静态强类型语言，可以有效地处理大量数据和复杂的运算。Fortran语言的语法和结构较为简单，主要由四个部分组成：声明、执行、程序和函数，其中声明部分包括变量、常量和数组等；执行部分包括循环、条件语句和子程序等。程序和函数则是Fortran程序的基本单位，程序就是一个完整的Fortran程序，而函数则只负责完成某一特定任务。Fortran还具有丰富的库函数和工具，如数值计算库（BLAS、LAPACK）、图形库（Gnuplot）和图形用户界面库（Qt），这些库函数和工具可以大大提高Fortran的可扩展性和功能性。

2 UMAT

2.1 UMAT概述

子程序UMAT是用于更新材料点，即单元积分点的应力和其他材料参数，返回该点的材料刚度矩阵。在UMAT中，应力和应变皆表示为一维向量，在固体材料中，应力应变一维向量为：

除应力与应变为，若需要额外的材料状态参数(如塑性应变 )，需要向ABAQUS要求额外的材料储存空间。在inp文件中，在关键字[*Material]下面，以*DePVAR定义状态变量个数，或在CAE定义材料时，于General下的Depvar内填写所需状态变量的个数。

Number of solution-dependent state variables指定状态变量的数量，分配存储空间。

Variable Number Controlling Element Deletion控制“不需要”的elements数量的一种方法。在工程设计和优化过程中，模型的elements经常会随着不断的修改和优化而发生变化。如果每次修改模型都需要重新生成网格，这样会浪费很多时间和计算资源。因此，通过使用VNC element deletion技术，可以使得不必要的elements自动删除，从而加快有限元分析的速度和效率。

激活Element Deletion：

（1）控制elements删除的状态变量号；

（2）在Step中，勾选State/Field/User/Time类别中的STATUS作为字段输出变量；

（3）在Mesh模块中，在Element Type中打开元素删除。

ABAQUS求解非线性问题一般采用的是增量-迭代思想，实施方法增量步控制法，称为Newton-Raphson方法。该方法通过将非线性问题分解为一系列小的增量步来逼近其解，并确保在每个增量步可经由迭代计算求得收敛的数值解，以得到非线性问题的最终解。

具体来说，增量步控制法包括以下几个步骤：

构建初始状态：给定材料和几何结构的初始状态。
划分增量步骤：将初始状态划分为若干个增量步骤，并确定每个增量步骤的时间长度、载荷大小和载荷施加方式等参数。
进行初次加载：在第一个增量步骤中施加载荷，并使用牛顿-拉夫逊迭代方法求解。
进行下一增量步骤：在前一个增量步骤的基础上，依次施加后续的增量载荷，并使用原始牛顿-拉夫逊法或修改的牛顿-拉夫逊法进行迭代求解。
判断收敛性质：对每个增量步进行计算，观察解的收敛性质，如果未达到收敛要求，则需要调整模型参数并重复执行前面的步骤，直到达到收敛。

2.2 UMAT基本架构

      subroutine umat(
C *** 子程序设计 ********************************************************* C
     1  STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     2  RPL,DDSDDT,DRPLDE,DRPLDT,
     3  STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     4  NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     5  CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C **** 定义参数类型 ****************************************************** C
      include 'ABA_PARAM.INC'
C -------------------------------------------------------------------- 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)
C --- 材料模型主程序 ---------------------------------------------------- C
C     提供 DDSDDE, STRESS, STATEV [, RPL, DDSDDT, DRPLDE, DRPLDT, SSE, SPD, SCD, PNEWDT]
C ----------------------------------------------------------------------- C
      return
      end subroutine umat

UMAT子程序中参数可分为三类：需定义的参数、可更新的参数、不可更新的参数。需定义的参数在每一个分析步中，根据材料模型更新。可更新的参数不直接影响分析结果，根据使用者需求更新和调整。不可更新的参数仅供其他两类参数更新用，其值的改动可能会被忽略，甚至程序报错。

2.2.1 需定义的参数

DDSDDE(NTENS, NTENS)

DDSDDE为雅克比矩阵，亦称为材料刚度矩阵。NTENS*NTENS，NTENS为该单元所对应应力/应变分量个数(如固体单元NTENS=6；壳体单元NTENS=3)。DDSDDE定义在单一迭代分析步中，表示：应力增量Δσ对应的应变增量 Δε的变化率矩阵，即：

ABAQUS默认并强制DDSDDE为对称矩阵，若想使用费对称材料刚度矩阵，须在User Material下勾选Use unsymmetric material stiffness matrix。

STRESS(NTENS)

STRESS为Cauchy应力一维向量。ABAQUS会将每一分析开始前时材料点的应力全量经过刚性和坐标轴旋转量 (DROT矩阵)更新后，传入子程序，使用者仅需针对变形增量的部分更新柯西应力。

STATEV(NSTATV)

定义 STATEV 数组，来保存需要跨越时间步骤传递的材料内部变量或状态量。这些变量通常包括材料的应力、应变、塑性应变、温度、损伤参数等，它们的值会随分析步更新。

在UMAT中，ABAQUS 会自动将当前分析步的STATEV数组传递给下一个分析步，并且将其作为输入参数传递到UMAT中进行计算，从而实现了材料的状态更新和演化。因此，在定义STATEV数组时，需要确定每个状态变量的初始值、更新方程以及如何在子程序中读写这些变量的值。

RPL

代表材料的塑性流动规律。RPL参数对应于材料的应力-应变曲线在塑性区域内的形状。RPL参数通常是一个指数功率函数，其中指数控制了应力-应变曲线的斜率，在不同的应变速率条件下可以改变该参数以更好地描述材料塑性行为。这个参数在有限元模拟中被用来模拟材料的非线性、塑性刚性行为。

DDSDDT

应力增量对温度增量的变化率

2.2.2 可更新的参数

SSE, SPD, SCD

SSE表示specific elastic strain energy，(比弹性应变能 ) ，通常也称为elastic strain energy density。比弹性应变能表示在材料发生塑性变形之前，由外部载荷引起的弹性变形所储存的应变能与材料的体积之比。具体而言，它是指单位体积内由弹性形变产生的势能密度，通常用 $J/m^{3}$ J/m^{3} 或Pa表示。其计算方式为将整个模拟中所有时间步长中的弹性应变能加总，并除以材料的体积。

SPD表示specific elastic plastic dissipation，(比塑性耗散能) ，比塑性耗散能表示由材料的塑性变形所引起的内部微观摩擦、滑移等机制所产生的耗散能密度，与材料的体积有关，通常用J/m3或Pa表示。它是指单位体积内由塑性形变产生的损耗能量密度。具体而言，比弹塑性耗散能是指整个模拟中所有时间步长中由材料的塑性变形所损失的能量加总，并除以材料的体积。

SCD表示specific creep dissipation，(比蠕变耗散能) ，比蠕变耗散能表示由材料的蠕变变形所引起的内部微观摩擦、滑移等机制产生的耗散能密度，与材料的体积有关，通常用 $J/m^{3}$ J/m^{3} 或Pa表示。

PNEWDT

分析步时间增量(DTIME) 调整系数，所对应的新时间增量值为PNEWDT*DTIME。

2.2.3 不可更新的参数

STRAN(NTENS)

分析步开始时的全应变(不包括温度应变)的一维向量

DSTRAN(NTENS)

应变增量一维向量

TIME(2)

分析步开始时，单一分析的时间步

DTIME

时间增量

TEMP

增量步开始时的温度

DTEMP

温度增量

PREDEF

预定义变量，是一些可以从外部输入的常数，例如材料参数、元素几何形状、外部环境条件等。UMAT程序需要使用PREDEF变量来获取这些预定义变量的值，并在计算中使用它们。这些预定义变量通常通过输入文件的方式输入，例如在ABAQUS中，可以使用材料属性定义（*MATERIAL）和实体定义（*ELEMENT）中的参数来定义这些预定义变量。

DPRED

预定义变量的增量

CMNAME

材料名称，需要使用多个自定义材料模型时，根据CMNAME的值来呼叫对应的材料模型。

NDI

表示应力张量中，主对角线上独立的应力个数。一般实体单元有3个独立分量，平面应力单元为2。在编写UMAT函数时，NDI变量的值会影响输入和输出的矩阵维度。对于2D的应变问题，输入和输出矩阵的大小为3×3，而对于3D应变问题，输入和输出矩阵的大小为6×6。

NSHR

表示应力张量中，非对角线或剪切应力部分的独立张量个数，一般实体单元为3个独立分量，平面应力单元为1个

NTENS

NDI和NSHR的和。一般实体单元为3+3，平面应力单元为2+1

PROPS(NPROPS)

用来储存使用者定义的材料参数

NPROPS

使用者定义的材料参数的个数

COORDS

材料点坐标一维向量，若考虑几何非线性，代表形变后的坐标；若不考虑，代表初始坐标

DROT(3,3)

DROT(3,3)是一个3x3的旋转矩阵，其中每个元素表示一个坐标系之间的旋转角度。在某些情况下，UMAT函数需要在不同坐标系之间进行计算，例如在材料的本地坐标系和全局坐标系之间进行转换。在这种情况下，需要定义一个旋转矩阵DROT，并在UMAT函数中使用它将应变张量或其他相应的张量从一个坐标系旋转到另一个坐标系中。

对于某些特定的加载路径和材料，可能需要在每个时间步长更新DROT矩阵的值。在这种情况下，UMAT函数需要计算新的旋转矩阵，并将其存储在DROT(3,3)中，以便在下一个时间步长中使用。