当前位置：服务支持 > 软件文章 > ABAQUS UMAT子程序趣味教程（二）：线弹性本构编写

ABAQUS UMAT子程序趣味教程（二）：线弹性本构编写

阅读数 8

各位小伙伴，大家好，我是艾米雷斯。

说实话，在写这个不正经教程的第二讲时，艾米还是有那么些许自我感动的。毕竟以艾米这么懒的性格，本来应该在拖了四五天后，甚至一两周后才开始写这篇稿子，这次竟然第二天就想写了，感觉我大学给心仪的女孩子发短信也没这么积极啊。

拖延症，这是人类行为学中的无解问题。艾米回想起博三的时候，买了一本书，叫做《战胜拖拉》，然而，我硬是到博士毕业时都没翻开那本书，这已经不是战胜拖拉了，而是被拖拉给教做人了。

闲话就先扯到这里，还是开始我们今天的正题，那个不正经教程的第二讲。毕竟要填一下上一讲给自己挖的坑，否则以后坑挖多了，容易崴脚。

想必小伙伴们都注意到了那个题目中的那根大粗棒子了，没错，今天，我们就是用它来辅助讲解我们的线弹性本构子程序的编写。这根棒子的尺寸是艾米拍屁股想出来的，截面是100毫米的正方形，长度是1500毫米。

你问我为啥用这个尺寸？这个嘛，拍屁股想出来的东西应该不用解释吧。

ABAQUS的具体建模过程艾米就不复述了，普通青年可以用图形用户界面（GUI）的方式建立模型，二 * 青年可以用Python自动化脚本来建立这个模型。

这里顺带提一句，UMAT不正经教程写到中途，艾米可能开一个“ABAQUS参数化建模中二青年训练营”的教程，愿意陪艾米一起中二的小伙伴们，届时可以一同前往围观。不过不要期待过深哦，说不定由于艾米太懒，到时候你们发现连这个知乎号都注销了，以逃避更新，嘿嘿。毕竟，处理“扔过墙帽子”的最完美方案，就是当这顶帽子不存在，重新买一顶。

虽然ABAQUS中的建模不是本教程的重点，不过，艾米还是要说一下，想要在ABAQUS中使用UMAT技术，有几个地方是需要处理，那就是，得让ABAQUS知道：你使用了UMAT技术。总结来说，就是对你建好的模型，做两个地方的处理：

****************************************************************************************

材料模块 (Property) 的处理
任务模块 (Job) 的处理

****************************************************************************************

材料模块的处理在以下框框中进行：

没错，就是你选择ABAQUS自带材料的地方。User material 在中部的那个 General 选项卡中可以找到。这里要注意，在用户材料定义中，一般是每行代表一个材料属性，这个材料属性代表啥，是你自己定的哦。

譬如，艾米这里就定义了两个材料属性，弹性模量206000MPa，和泊松比0.3 （很明显这是钢材的常规材料性能）。上面那个Depavar是求解过程中你自己定义的储存变量，目前还用不到，用到了艾米再和大家说，你用的时候也可以删除。

后面我们会涉及到塑性，往往小伙伴需要输入自己获得的材料性能数据，譬如，应力应变关系曲线，有十组点。这里无法输入两列数，所以你可以定义2~11行对应了十组点的应力，而12~21行对应应变。

这里艾米要提一句哦，定义在哪一行，甚至怎么定义（譬如奇数行是应力，偶数行是应变），都是可以的，为什么呢？因为这里ABAQUS只是让你打包自己的数据。

有个形象的例子：

你毕业了，需要将自己学校寝室书架上的书搬回家，并且家里的摆放顺序与学校寝室相同（重度强迫症鉴定完毕！），当把书放进搬家公司箱子里的时候，你当然可以随意放，只要记得哪本书在哪里就行了，回家后，你再按照你的记忆将书都拿出来，恢复顺序即可。当然，艾米一般都是把这些书都卖了。

扯远了，总结一下就是，这些数据怎么放都是小伙伴们自己决定，到时候在程序里面怎么用，你自己心里有数就行。

任务模块的处理在以下框框中进行：

这一步就是告诉ABAQUS，你编写的.for文件在什么地方。

处理好这两个地方，下面我们的主要关注点，就是我们的UMAT子程序了。

这里提一句哦，艾米下面用的公式啥的，都参考了下面这本书，陈惠发老先生的《弹性与塑性力学》：

不知道小伙伴们有没有这样一个错觉，为啥我看懂这些教材了，但是编写UMAT程序的时候，还是看不懂别人编的程序？

能看懂艾米就呵呵了，要知道每个人的理解都是不一样的，这会直接导致他在生成刚度矩阵时的思路都不一样。这也是艾米在看别人教程时的感受：

为了有一个统一的思路，艾米决定采用陈惠发老先生的书。嗯，艾米是不会告诉你，真实原因是艾米懒得再去读其他教材。

下面我们就开始我们第一个UMAT子程序之旅吧。

还记不记得，上次我们提到了UMAT的架构中，艾米有一块的东西是没有解释的，就是下面这些东西：

****************************************************************************************

CHARACTER*8 CMNAME

DIMENSION STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS,NTENS),

1 DDSDDT(NTENS), DRPLDE(NTENS), STRAN(NTENS), DSTRAN(NTENS),

2 PREDEF(1), DPRED(1), PROPS(NPROPS), COORDS(3), DROT(3,3),

3 DFGRD0(3,3), DFGRD1(3,3)

****************************************************************************************

有没有发现，什么 STRESS 啦， STATEV 啦，都是最上面 UMAT 子程序括号里出现的。

小伙伴们不要觉得恐慌哦，其实，这是ABAQUS在定义一些储存的变量。

我们知道，对于子程序最上面的声明， SUBROUTINE UMAT(STRESS,STATEV, ... )，括号里面的这些变量实际上是 UMAT 的子程序自变量。但是ABAQUS不会自动在计算机中建立一些区域，存储这些变量的值，供你直接使用。这类似于函数的形式参数（形参）。

所以这个地方，讲白了，就是ABAQUS让你自己定义下这些变量的类型。举个例子：

STRESS(NTENS)，其实就是定义了一个长度为NTENS的数组，STRESS是变量名。NTENS是啥呢？它是正应力和剪应力的数量之和，这个ABAQUS会根据单元的类型自动传输进来，譬如实体单元，NTENS的数值就是6。

由于STRESS是数组，所以必须用 DIMENSION 来声明，那么 DDSDDE(NTENS,NTENS) 就很容易理解了，当 NTENS = 6时，它就是一个6行6列的数组。

实际上，你不一定需要 DIMENSION STRESS(NTENS) 这样来定义，如果你确定只使用实体单元，DIMENSION STRESS(6) 也是可以的。

这里提一下数组的初始化，一般你定义了一个数组，不作操作的话，数组的元素都是0。

但是，根据艾米的操作经验，这里有一个小坑：

如果你自己定义了一个数组，譬如 AAA(NTENS,NTENS)，如果你不对它做任何操作，按道理，它应该时时刻刻每个元素都是零，但是你会发现，实际情况是，这样自己定义的数组，运行到后面，ABAQUS写入一些很奇怪的数值。

艾米也不知道为啥，但测试后有两种解决方案：

****************************************************************************************

这么定义： DIMENSION AAA(6, 6)，也就是直接填入数字
每次使用前，都将该变量初始化为零，AAA = AAA * 0.D0

****************************************************************************************

肯定有小伙伴要问了，“0.D0”这是啥？

不慌，这是FORTRAN 的数值表达方式，D 代表Double，表示这是一个双精度实数，然后采用类似于科学计数法，0.D0，就是0，因为是 0.0 * 10^0 ; 1.0D-6就好理解了，就是 1.0 * 10^(-6)，也就是0.000001。

其实以上这些变量定义照抄就行，相当于ABAQUS让你做个罐子，这样在给你“圣水”的时候，你有东西接。

这里再多插一句，如果小伙伴们要定义自己的数组，请一定要放在ABAQUS这些语句之后，主程序语句开始之前，一旦你开始了一些变量的赋值，计算之后，ABAQUS就不允许你再利用DIMENSION定义你的数组，不然会报错。

下面就开始顺序编写我们的UMAT子程序。

首先是将ABAQUS传进来的数值放进我们自己定义的变量中，如下所示：

****************************************************************************************

EMOD = PROPS(1)

ENU = PROPS(2)

AG = EMOD / 2.0D0 / (1.0D0+ENU)

AK = EMOD / 3.0D0 / (1.0D0-2.0D0*ENU)

****************************************************************************************

注意哦，Fortran里面有一个 I-N 准则，也就是你定义变量的时候，如果开头字母位于 I 和 N 之间，包括：I J K L M N，此时，Fortran会自动认为这些是整型变量，会自动修改这个变量的类型为整型。有小伙伴说，可以用 IMPLICIT NONE这个命令，但艾米的实践经验表明，在固定格式情况下加入这个语句，有时候会玄学报错，所以建议不加，直接记住这个规则就行。

你看艾米定义的变量就是 E 和 A 开头的，机智吧。

给没有编程经验的小伙伴科普一个概念，一个整型变量和整型变量的除法运算，其结果是一个整型变量，譬如 1/2 = 0，没错，确实是0；所以想要获得0.5的正确计算结果，你至少得让其中一个数，1 或者 2 为实型变量。这就是为毛艾米要在此处和大家唠嗑 I-N 准则的原因。

EMOD = PROPS(1) 这个命令就是将我们在Property模块输入的第一行的数值读进来，这里就是206000，是弹性模量。ABAQUS实际上，就是将我们输入的那个数字，传进了它自己预设的一个数组 PROPS，你按照自己之前设计的顺序提取就行。

这里艾米再插播一个小知识：

****************************************************************************************

临时的数值型的变量可以在程序的任何地方定义

****************************************************************************************

什么意思呢？

譬如，我想定义一个变量 A，这个变量不是数组，只是一个数，那么你可以在程序的任何地方如下定义，而不必非得在程序开始前，和 DIMENSION 一样进行定义：

A = 1.0D0

还有什么好处呢？

你的程序在任何时刻计算出来的数值，都可以放进A，而且你不用提前告知程序，这个A是实型的还是整型的，你输入的是啥，A的类型就是啥。这种功能在循环中，你定义一个临时的变量时，极为有用。

ENU = PROPS(2) 这个也就好理解了，就是读取到了第二个数，泊松比 0.3。

有小伙伴要问，可不可以不定义啥ENU，直接使用 PROPS，当然可以，但你这不是吃饱了撑的么，你后面再回头检查程序的时候，PROPS(2) 说不定就忘了这是啥意思了，懂？

AG = EMOD / 2.0D0 / (1.0D0+ENU) 和 AK = EMOD / 3.0D0 / (1.0D0-2.0D0*ENU) 这两条语句是做什么用的呢？

答案是，用来定义线弹性本构的两个材料参数。

这里再插播一个对 UMAT 本质工作的介绍，我们编写UMAT，目标是做啥呢？其实主要目标就是两个：

****************************************************************************************

更新 UMAT 中的应力数组 STRESS
更新 UMAT 中的应力应变矩阵 DDSDDE

****************************************************************************************

UMAT中，STRESS 数组一共有6个分量（这里我们用实体单元），排列顺序为：

$\left[\sigma_{11},\sigma_{22},\sigma_{33},\sigma_{12},\sigma_{13},\sigma_{23} \right]$ \left[\sigma_{11},\sigma_{22},\sigma_{33},\sigma_{12},\sigma_{13},\sigma_{23} \right]

分别是三个方向的正应力和三个方向的剪应力。

还有一个应变增量数组 DSTRAN，里面的排列顺序是：

$\left[d\varepsilon_{11},d\varepsilon_{22},d\varepsilon_{33},d\gamma_{12},d\gamma_{13},d\gamma_{23} \right]$ \left[d\varepsilon_{11},d\varepsilon_{22},d\varepsilon_{33},d\gamma_{12},d\gamma_{13},d\gamma_{23} \right]

这里有一个坑哦：

$d\gamma_{12}=2d\varepsilon_{12}, d\gamma_{13}=2d\varepsilon_{13},d\gamma_{23}=2d\varepsilon_{23}$ d\gamma_{12}=2d\varepsilon_{12}, d\gamma_{13}=2d\varepsilon_{13},d\gamma_{23}=2d\varepsilon_{23}

讲白了就是，DSTRAN 里面用的是工程剪应变，然而到了 VUMAT子程序后，又不用了，这不是玩艾米呢。不管怎么样，小伙伴们要注意这个坑哦。

那DDSDDE是啥，其实就是：

$\left[d\sigma_{11},d\sigma_{22},d\sigma_{33},d\sigma_{12},d\sigma_{13},d\sigma_{23} \right]$ \left[d\sigma_{11},d\sigma_{22},d\sigma_{33},d\sigma_{12},d\sigma_{13},d\sigma_{23} \right] 和 $\left[d\varepsilon_{11},d\varepsilon_{22},d\varepsilon_{33},d\gamma_{12},d\gamma_{13},d\gamma_{23} \right]$ \left[d\varepsilon_{11},d\varepsilon_{22},d\varepsilon_{33},d\gamma_{12},d\gamma_{13},d\gamma_{23} \right] 之间的变换矩阵。

注意了，这里艾米在应力符号之间也加了"d"，代表应力增量。

UMAT 就是要你提供这两个数组的定义，所谓定义，就是这里面元素都取啥值。

变量 AK 和 AG 就是为定义 DDSDDE 这个数组的元素而准备的，我们接着往下说。

是不是看到这里有点想关掉页面，先收藏,下次再看？毕竟开始涉及公式了

然而艾米只想说：

怎么定义，我们就要用到陈惠发老先生书里的弹性理论了

众所周知，各向同性线弹性材料的独立材料常数只有两个，最源头的我们称为 Lame 系数，当然你也可以定义其他两个系数，这两个系数和两个Lame系数之间存在函数关系即可。艾米这里用的就是剪切模量 G，和体积模量 K，也就是变量 AG 和 AK。

为啥艾米用这两个系数？不要问，问就是量子力学。

当然，小伙伴看到那个“众所周知”也不要虎躯一震，你不知道很正常，艾米写“众所周知”主要是让自己看起来“很专业”，实际内心也慌得一匹。

那按照陈老先生的理论， $\left[d\sigma_{11},d\sigma_{22},d\sigma_{33},d\sigma_{12},d\sigma_{13},d\sigma_{23} \right]$ \left[d\sigma_{11},d\sigma_{22},d\sigma_{33},d\sigma_{12},d\sigma_{13},d\sigma_{23} \right] 与 $\left[d\varepsilon_{11},d\varepsilon_{22},d\varepsilon_{33},d\gamma_{12},d\gamma_{13},d\gamma_{23} \right]$ \left[d\varepsilon_{11},d\varepsilon_{22},d\varepsilon_{33},d\gamma_{12},d\gamma_{13},d\gamma_{23} \right]之间的 DDSDDE 该怎么写呢？如下所示：