用户子程序和应用程序（）

4.3.2　用户子程序UMAT和应用程序（Utilities）

1）用户子程序（User Subroutines）和应用程序（Utilities）

虽然ABAQUS为用户提供了大量的单元库和材料模型（如金属、橡胶、塑料、混凝土、岩土等），使用户能够利用这些模型处理绝大多数的问题，但上述道路工程（包含部分岩土工程）中常见材料模型，并没有包含在ABAQUS软件中。为了弥补这一不足，ABAQUS提供了大量的用户子程序（User Subroutines）和应用程序（Utilities），用户可以自行定义符合特定问题的模型。

（1）用户材料子程序UMAT

用户材料子程序UMAT（User－Defined Material Mechanical Behavior）是ABAQUS提供给用户进行材料本构模型二次开发的一个用户子程序接口，可以定义用户所需要的各类材料本构模型，这大大增强了ABAQUS的应用面和灵活性。

用户子程序UMAT具有如下特点：

·用来定义材料的本构关系；

·当材料的定义包含用户自定义材料模型时，每一个计算单元的材料积分点都可以调用UMAT；

·可以用于力学行为分析的任何分析过程；

·可以使用状态变量；

·对于力学本构关系，必须在UMAT中提供材料本构模型的雅可比矩阵

·可以和用户子程序USDFLD联合使用，通过USDFLD重新定义任何常变量值并传递到UMAT。

UMAT子程序的核心内容就是给出定义材料本构模型的雅可比矩阵（Jacobian矩阵，即应力增量对应变增量的变化率aΔσ

aΔε，Δσ是应力增量，Δε是应变增量），并更新应力提供给ABAQUS主程序。例如，已知第n步的结果σ_n及ε_n等，然后ABAQUS主程序给出一个应变增量dε_n＋1，UMAT根据提供的雅可比矩阵DDSDDE计算出新的应力σ_n＋1。

UMAT子程序采用Fortran语言编写，从主程序获取数据，计算单元的材料积分点的雅可比矩阵，并更新应力张量和状态变量，UMAT的接口格式如下：

SUBROUTINE UMAT（STRESS，STATEV，DDSDDE，SSE，SPD，SCD，RPL，DDSDDT，

1DRPLDE，DRPLDT，STRAN，DSTRAN，TIME，DTIME，TEMP，DTEMP，PREDEF，DPRED，

2CMNAME，NDI，NSHR，NTENS，NSTATV，PROPS，NPROPS，COORDS，DROT，

3PNEWDT，CELENT，DFGRD0，DFGRD1，NOEL，NPT，LAYER，KSPT，KSTEP，KINC）

C

INCLUDE′ABA＿PARAM.INC′

C

CHARACTER＊80CMNAME

C

DIMENSIONSTRESS（NTENS），STATEV（NSTATV），DDSDDE（NTENS，NTENS），

1DDSDDT（NTENS），DRPLDE（NTENS），STRAN（NTENS），DSTRAN（NTENS），

2TIME（2），PREDEF（1），DPRED（1），PROPS（NPROPS），COORDS（3），DROT（3，3），

3DFGRD0（3，3），DFGRD1（3，3）

User coding to define DDSDDE，STRESS，STATEV，SSE，SSD，SCD

and，if necessary，RPL，DDSDDT，DRPLDE，DRPLDT，PNEWDT

RETURN

END

UMAT子程序主要变量说明：

·DDSDDE（NTENS，NTENS）

是一个NTENS×NTENS维的方阵，称为本构关系的雅可比矩阵，DDSDDE（I，J）表示增量步结束时第J个应变分量的改变引起的第I个应力分量的变化。通常雅可比矩阵是一个对称矩阵。

·STRESS（NTENS）

应力张量数组，对应NDI个直接分量和NSHR个剪切分量。增量步开始时，该数组从ABAQUS主程序获取数据σ_n并作为已知量传入UMAT；增量步结束时，在子程序中更新应力张量为σ_n＋1。如果定义了初始应力，那么分析开始时该应力张量的数值即为初始应力。对于包含刚体转动的有限应变问题，一个增量步调用UMAT之前就已经对应力张量进行了刚体转动，因此在UMAT中只需处理应力张量的共旋部分。UMAT中应力张量的度量为柯西（真实）应力。

·STATEV（NTENS）

状态变量数组，用于存储状态变量的数组，数值在增量步开始时从主程序传递到UMAT，也可以在USDFLD或UEXPAN中先更新数据，然后增量步开始时将更新的数据传递到UMAT子程序中。在增量步结束时，必须更新STATEV中的数据返回给主程序。状态变量数组的维数通过ABAQUS输入文件中的关键字“＊DEPVAR”定义，关键字下面数据行的数值即为状态变量数组的维数。状态变量数组是用来保存用户自己定义的一些变量，如累计塑性应变、黏弹性应变等。

·STRAN（NTENS）

总应变数组，存储增量步开始时的总应变ε_n，由ABAQUS主程序自动更新。

·DSTRAN（NTENS）

总应变增量dε_n＋1。

·PROPS（NPROPS）

材料常数矩阵，用来保存用户自定义材料参数，即用户自定义材料本构关系的模型参数。材料常数的个数NPROPS等于关键字“＊USER MATERIAL”中的“CONSTANTS”参数设定的值，矩阵中元素的数值对应关键字“＊USER MATERIAL”下面的数据行。

·SSE，SPD，SCD

分别定义每一增量步的弹性应变能、塑性耗散和蠕变耗散。它们对计算结果没有影响，仅仅作为能量输出。

·DTIME

时间增量dt。

·NDI

法向应力、应变个数，三维问题、轴对称问题是3（11，22，33），平面问题是2（11，22）。

·NSHR

剪切应力、应变个数，三维问题是3（12，13，23），轴对称问题是1（12）。

·NTENS

应力和应变分量数组的大小（NTENS＝NDI＋NSHR）。

·DROT

对Finite strain问题，应变应该排除旋转部分，该矩阵提供了旋转矩阵。

·PNEWDT

可用来控制时间步的变化。如果设置为小于1的数，则程序放弃当前计算，并用新的时间增量DTIME×PNEWDT作为新的时间增量计算；对时间相关的材料如聚合物等有用；如果设为大于1的数，则下一个增量步加大DTIME为DTIME×PNEWDT。

·KSTEP，KINC

ABAQUS传到UMAT的当前分析步和增量步数值。

·TIME（1），TIME（2）

当前分析步时间（STEP TIME）和增量步时间（INCREMENT TIME）值。

（2）应用程序（Utilities）

ABAQUS提供了一些实用的应用程序（Utility Routines），供用户在开发用户子程序UMAT时调用，这些程序包括应力不变量的计算，主应力的计算等。

·SINV：用于计算应力不变量

CALLSINV（STRESS，SINV1，SINV2，NDI，NSHR）

STRESS：应力张量；

NDI：法向应力分量的数量；

NSHR：剪切应力分量的数量；

SINV1：第一应力不变量，SINV1＝traceσ，其中σ为应力张量；

SINV2：第二应力不变量，SINV2＝，其中S为偏应力张量，S＝σ－traceσI；

·SPRINC：用于计算主应力值

CALLSPRINC（S，PS，LSTR，NDI，NSHR）

S：应力或应变张量；

PS（I），I＝1，2，3：三个主应力值；

LSTR：标识，LSTR＝1表示S为应力张量，LSTR＝2表示S为应变张量。

·SPRIND：用于计算主应力值和主应力方向

CALLSPRIND（S，PS，AN，LSTR，NDI，NSHR）

S：应力或应变张量；

AN（K1，I），I＝1，2，3：PS（K1）法向应力的方向余弦。

2）ABAQUS主程序与UMAT子程序协同工作过程

ABAQUS主程序与UMAT子程序之间是一个动态交互传递数据、协同工作的过程。UMAT子程序作为ABAQUS主程序的一个接口，在单元的积分点上调用，增量步开始时，主程序通过UMAT的接口进入UMAT子程序，单元当前积分点必要变量的初始值将随之传递给UMAT子程序的相应变量，UMAT子程序计算单元材料积分点的雅可比矩阵，并更新应力张量和状态变量，最后将这些变量的更新值通过接口返回主程序。

ABAQUS主程序与UMAT子程序的交互计算过程如下：从第t_n时刻开始，ABAQUS在Δt时间增量内产生一个由外载荷产生的应变增量Δε，UMAT子程序通过给定的本构方程为主程序提供新的柯西应力张量σ（t_n＋Δt），如果计算的应力－应变结果收敛，那么ABAQUS主程序继续计算第t_n＋1步，并根据上一步的收敛情况来选取下一步增量步长。雅可比矩阵即DDSDDE的精度影响程序的收敛速度，但是并不影响计算结果的准确性。

以下节研究的Burgers蠕变模型为例，某一材料单元积分点上ABAQUS主程序与UMAT子程序协同工作过程，如图4.12所示。

图4.12　ABAQUS主程序与UMAT子程序协同工作过程

第一步：平衡时刻t_n，ABAQUS主程序提供给UMAT子程序总应变ε（t_n）和应力σ（t_n），同时ABAQUS主程序自动产生一个时间增量Δt。

第二步：ABAQUS主程序产生一个总应变增量Δε（t_n），调用UMAT子程序计算Δε（t_n）对应的蠕应变增量Δε_c（t_n）及应力增量Δσ（t_n），此过程是一个非线性迭代求解的过程。迭代收敛后更新应力张量σ（t_n＋1）＝σ（t_n）＋Δσ（t_n）。

第三步：UMAT子程序将更新后的应力σ（t_n＋1）返回给ABAQUS主程序，同时ABAQUS主程序自动更新总应变ε（t_n＋1）＝ε（t_n）＋Δε（t_n）。

第四步：ABAQUS主程序进行最大迭代次数的检查（控制最大迭代次数是出于计算成本的考虑），如果迭代次数超过了限定的最大迭代次数（ABAQUS默认的最大迭代次数为16次），那么ABAQUS主程序就会自动减小时间增量，然后返回到第一步重新进行迭代平衡；然后ABAQUS主程序将σ（t_n＋1）代入系统平衡方程进行平衡判断，若满足系统平衡标准，结束本增量步的迭代，进入下一增量步。若不满足系统平衡，此时ABAQUS主程序将放弃本次的应力、应变更新，回复到增量步初始时刻的值，并继续进行本增量步的下一次平衡迭代，直到满足系统平衡标准为止。