黏弹性本构：从一维模型到三维的UMAT数值实现

0. 开篇

（其实这里标题应该叫“引言”？）但我不想按论文的风格写的那么严肃。

出于最近的科研需要，小喵完整学习了一下关于黏弹性本构的知识。小喵获取的这些知识散落在许多来源，很多课本里的推导要么不完整，要么甚至有错误。所以我学的也是磕磕绊绊。或许之前有老师或者教材完整讲过，但至少我在调研学习的时候没找到。本着造福后人的精神，希望其他对黏弹性相关理论感兴趣的朋友，在搜索到、读过这篇文章以后，学习过程能稍微轻松一点。

这篇文章，和过去小喵写的许多推送一样，也是一篇我自己的学习笔记。略有不同的是，之前学习笔记性质的推送，很多与仿真软件操作有关的内容我是边学边写整理出来的。而这一篇，由于大部分是理论和公式，我是在自己完全看懂以后，再梳理出正确的逻辑一气呵成写下来的。这样对读者更友好。

这篇文章只在开头有一个简单的黏弹性Abaqus仿真案例。更多还是以公式推导为主，从一维讲到三维，一直讲到线性黏弹性本构的UMAT实现。我使用了mdnice墨滴编辑器的Markdown语法和LaTeX公式，让文章的公式排版看起来更清爽。

顺带一提，幸好最新版本的Word同时支持Unicode和LaTeX格式的公式，我可以在Word和OneNote里用简洁的Unicode语法敲公式，然后把它转成LaTeX放进Markdown。有些我在OneNote里推过的公式，就不需要再用LaTeX语法重新敲一遍了。

小喵是很菜的博士生，在推导公式方面更是没什么天分。这篇文章的公式推导我会尽可能讲的慢一点，详细一点，以匹配我这种水平选手的理解能力。希望对这个主题感兴趣的读者也都能够读懂。

正文目录：

1. “黏”还是“粘”？

2. Abaqus黏弹性案例

3. 一维黏弹性模型

4. Abaqus中的黏弹性参数

5. 微分与积分形式的桥梁：Boltzmann叠加原理

6. 一维到三维的桥梁：黏弹性本构向三维的扩展

7. 黏弹性本构的UMAT实现与对比验证

8. 书籍资料勘误及参考文献

在正文内容安排的顺序上，

首先用一个Abaqus案例来引入。不会介绍每一步具体操作流程，只是让读者对黏弹性材料初步有一个感性认知。

接下来，分别介绍 （1）经典的一维黏弹性模型，（2）Abaqus中黏弹性参数对应的公式，和 （3）Abaqus文档中黏弹性UMAT的控制方程。到这一步，读者会发现这三者的公式看上去长得似乎很不一样。

那么后续，就会通过Boltzmann叠加原理架起形式（1）和（2）之间的桥梁；再通过一系列推导，将一维的黏弹性模型扩展到三维，将形式（1）和（3）统一起来。

最后，介绍在Abaqus中，使用UMAT实现黏弹性本构模型的公式并给出源代码。与Abaqus内置的本构模型作对比验证，从而将形式（2）和（3）也统一起来。

全文一万多字，还挺长的。如果对这个主题的相关推导感兴趣，建议收藏慢慢看。

1. “黏”还是“粘”？

在正文开始之前，先小小的咬文嚼字一下。在现代汉语字典中，“黏”字一般读nián，“粘”字一般读zhān。当粘字读nián时，字典里的标记是【同“黏”】。所以在类似“黏稠、黏糊、黏性、黏液、黏土”等词的写法上，正体字应该是“黏”。

在《大辞海（材料科学卷）》中，相关描述使用的字都是“黏”。

 

 

综上，根据字典，实际写作“黏弹性”和“粘弹性”都是正确的。在Abaqus软件界面上翻译为“粘弹性”，但在这篇推送里，我还是统一使用“黏弹性”的字形。

那么，让我们正式开始。

2. Abaqus黏弹性案例

所谓黏弹性材料，简单来说就是：给它加一个恒定不变的力，随着时间的推移它的位移会逐渐增大；给它一个恒定的变形，那随着时间推移反力会逐渐减小。

在进行理论推导之前，先用一个仿真案例来带大家直观地认识一下黏弹性行为。

案例来自Abaqus文档中的示例问题（在帮助文档的这个位置能找到官方inp文件）：

 

Abaqus | Example Problems | Static Stress/Displacement Analyses | Static and quasi-static stress analyses | Transient loading of a viscoelastic bushing

官方给的inp文件使用了一个小技巧，用了*NGEN 和*NCOPY 两个Abaqus/CAE界面不支持的关键字，生成了一圈圆环形的平面应变垫圈网格。这样的inp文件直接求解可以正确执行，但导入Abaqus/CAE环境以后会有点小问题。导入HyperMesh甚至直接不能正确生成网格。小喵自己重新做了一版inp，下载方式见文末。

这些都不重要。前处理我们只看三件事：1）材料定义；2）分析步；3）边界条件，后处理再看个结果云图动画。其中，关于黏弹性，最重要的就是材料定义。我们看一看对黏弹性材料做仿真 需要提供什么样的参数。

材料定义方面，在Abaqus里，只允许在线弹性和超弹性材料下面叠加黏弹性材料属性。这里由于变形较大，材料本身定义为超弹性，应变势能为多项式形式，模量时间尺度为瞬态。

 

 

分析步方面，这个案例定义了四个分析步。两个 静力,通用 分析步用来施加载荷，两个 黏性 分析步用来模拟材料蠕变。Abaqus的粘性分析步里，蠕变/膨胀/粘弹性应变错误容差 指定为0.0005，粘弹性积分方式为显式/隐式。

 

载荷方面，圆环中心定义一个参考点与圆环内圈耦合。Step-1施加向右50000的集中力，Step-3施加10000的弯矩。圆环外圈保持固定。Abaqus本身没有单位制，作为一个演示案例，我这里也就不提具体单位了。

 

计算结果放一张动图出来。注意，两段位移突变都来自 静力通用 分析步；比较流畅的动画都来自 黏性 分析步。在这两段动画发生的时间里，载荷都没有变化，位移变化全部来自材料的黏弹性行为。

这里有必要强调一下，Abaqus内置的黏弹性本构，在 静力通用 分析步中是不会发生黏弹性行为的。静力分析步的时间没有物理意义，无论静力分析步的时间是0.001还是1，计算结果都一样。只有黏性分析步的时间才会让材料发生黏弹性的蠕变。（而Abaqus官方提供的黏弹性UMAT示例子程序里则无此特征，即使是静力通用分析步，也会有黏性行为。当然也可以手动修改，后文会讨论）

我们在这个案例中看到，Abaqus中定义材料黏弹性行为，需要提供三种参数：

可以写很多行，即     可以取1,2,3,...

有花花绿绿云图的案例到此结束。后面就基本都是公式推导啦。其实也就是大学本科难度再稍微扩展一丢丢。我会尽量把每一步写的清楚一些，读者认真跟着读下来的话，相信是一定可以读懂的。

在后文中，我们先介绍一维的黏弹性简化力学模型，再介绍Abaqus中这三个参数对应的公式。

3. 一维黏弹性模型

让我们从两种最简单的一维概念模型开始。它们分别是Maxwell模型和Kelvin模型，后者有时也称为Kelvin-Voigt模型。

 

注意到模型里面有一个黏壶，或者也可以叫它阻尼器。阻尼器提供的反力与变形的速度成正比。使用应力和应变来描述的话，设其阻尼为     ，则有：

（沿袭牛顿的记法，变量头顶加点表示对时间求导。V表示visco黏性，E表示elastic弹性）

最简单的Maxwell模型和Kelvin模型，分别是将一根弹簧与一个阻尼器串联和并联。在串联时，左右两部分应力相等，总应变等于二者叠加：

因此，对于Maxwell模型，其控制微分方程为：

在弹簧和阻尼器并联时，上下两部分应变相等，应力则为二者叠加：

因此，对于Kelvin模型，其控制微分方程为：

对黏弹性模型，最常见的两种测试分别是蠕变和松弛。即：1）瞬间施加一个恒定的力之后保持不变，观察其位移随时间的变化——蠕变；和 2）瞬间施加一个恒定的位移之后保持不变，观察其反力随时间的变化——松弛。

根据维 基百科，Maxwell模型可用于描述软固体，如靠近熔化温度的热塑性聚合物、新混凝土（忽略其老化），以及温度接近熔点的各类金属；Kelvin模型可用于负载不太大的有机聚合物、橡胶和木材。

这里由于篇幅原因先不展开详细推导，相信读者很容易直观想象这两种模型的蠕变和松弛行为以及其问题：Maxwell模型在蠕变条件下，位移会无限增大、松弛条件下反力最终会降到零；Kelvin模型在蠕变条件下初始位移为零、在松弛条件下初始反力几乎无限大，后续反力几乎不变。

为了解决这两种模型的问题，提出了标准线性模型，也称为Zener模型。（冯元桢教授的书中把弹簧阻尼器并联的模型称为Voigt模型，而把标准线性模型称为Kelvin模型……无所谓了，反正画出示意图不至于引起混乱就好）

Zener模型有两种描述方式。分别是1）在Maxwell模型上并联一个弹簧，和2)在Kelvin模型上串联一个弹簧。其控制方程分别为：

 

抄写下来。对于Maxwell形式的标准线性模型，控制方程为：

对于Kelvin形式的标准线性模型，控制方程为：

这两个控制方程的推导，还是需要稍微花一番功夫的。感兴趣的读者可以暂停自己试着推导一下。我们这里以Kelvin形式为例，简单做一下解释。

对于Kelvin形式的标准线性模型 ，左侧有：

右侧有：

左右(3-5)和(3-6)式联立：

(3-7)式里还有右侧的应变     ，这不行。得去掉：

整理可得上面的公式(3-4)：

用类似思路就可以得到Maxwell形式的控制方程(3-3)。

这两种形式虽然表达式看起来不太一样，但控制方程里都包含了应力、应力率、应变、应变率     这四项。所以其实在本质上，这两个模型是等价的。

除此之外，还有Burgers模型如下图所示。这里的控制方程居然出现了应力和应变的二阶导，说实话小喵也没太理解这是为什么，不过这个模型不太重要。

 

以及广义Maxwell模型，这是最通用的模型，也是各类有限元软件使用的模型。理论上可以用它模拟出前面所有的简化模型。

 

请读者先对这些一维黏弹性模型留个印象。我们继续看Abaqus里面黏弹性参数对应的公式长什么样。

4. Abaqus中的黏弹性参数

其实，不仅仅在Abaqus里黏弹性参数是这个公式。几乎所有商用有限元软件里，定义黏弹性材料的方式都大同小异。

本节的推导和符号参考了Abaqus帮助文档。

先回顾一下各向同性线弹性理论。

 

我们知道，对于线弹性模型来说，本构关系用张量形式可以写为：

 

 

   一起被称为Lamé常数。其中的     同时也是固体力学中的剪切模量     。    里面，    是哑指标表示求和，    在i=j时为1，否则为0。

如果我们把应变张量分解成球张量和偏张量:

 

 

 

则本构关系可以写为：

 

 

其中，    为体积模量，    为剪切模量。(4-2)这种形式是将体积变形和剪切变形区分开来的写法。弹性常数之间的关系，找个最全的表格放在这：

 

回顾结束。在Abaqus文档中介绍，三维空间中，广义Maxwell本构模型的公式可以写为：

很显然这里     .

(4-3)这个形式与线弹性本构的(4-2)式很像。

其中，    和     是Prony级数形式的剪切松弛模量和体积松弛模量，它们都是时间的函数。以剪切松弛模量为例，可以有许多种写法。例如：

这里，    就是长期剪切模量，后面指数里的     是对应第i项的剪切弛豫时间，含义是让第i项的黏性反力减小到初始的    倍所需的时间。式(4-4)的变体可以写为（为了清晰，后面e指数函数写为exp):

这里小写的     就是和瞬时剪切模量     相比，第i项的相对剪切模量。

类似地，体积松弛模量也可以写为（我这里就偷懒只写出一种形式了）：

在Abaqus中，认为     。这个其实没关系，想要让剪切黏弹性行为和体积黏弹性行为相互独立的话，把对应项的    或者    单独设为0就好了。

对应前面第2节的内容，在Abaqus里定义材料的黏弹性行为，需要提供的参数就是公式(4-6)和(4-7)里面的     。

在Abaqus中，定义线弹性和超弹性本构时，都有一个模量时间尺度选项，可选长期或瞬态，就分别对应了上面黏弹性松弛模量中的     和     。

到这为止，认真跟着推导看下来的读者应该很容易理解前面的(4-4)至(4-7)式。对那个e指数函数的理解，分别把时间为0和时间无穷大代入进去就很直观了。但是，相信大多数没有黏弹性相关知识基础的朋友都会有个疑惑：

这些公式，和前面一维的黏弹性模型，是怎么建立起联系的呢？

(4-3)(4-6)(4-7)和前面Maxwell模型的控制方程(3-3)，它们看上去长得并不像啊。

 

这还不是全部。还有第三个长得不像的东西，在Abaqus文档的UMAT章节：

 

Abaqus | User Subroutines | Abaqus/Standard User Subroutines | UMAT | Example: Simple linear viscoelastic material

里面，介绍了对应于前面Kelvin形式的标准线性模型，一维控制方程为式(3-4)，

在小应变、各向同性条件下扩展到三维，本构关系可以写成这样：

 

这公式(4-8)，虽然它长得倒也很像线弹性的式(4-1)，但怎么看上去好像和Abaqus内置的黏弹性本构(4-3)也都不太一样呢。

而且，不知道各位读者面对这一页截图，你自己能不能试一下从(3-4)如何推导出(4-8)，并且给出     这五个参数与Maxwell形式的     之间的转换关系？

反正，小喵资质驽钝，没有拉马努金那种级别的“注意力”。所以下面的段落，我们就来仔细讲一讲，上述三种形式：1）一维形式的黏弹性弹簧-阻尼器模型；2）Abaqus内置黏弹性材料的剪切和体积松弛模量形式；3）Abaqus UMAT例子中的三维黏弹性本构模型，这三者是如何相互联系和统一起来的。

5. 微分与积分形式的桥梁：Boltzmann叠加原理

首先从简单Maxwell模型开始。它的微分方程为(3-1)：

如果施加一个恒定的应变，考察其松弛行为，那么在后续的时间里，应变保持不变，    . 满足这一条件的通解形式可以写为：

式(5-1)很容易验证。这里面有两个待定常数，如果在     时刻，施加数值为     的恒定应变，那么施加时刻可以确定     ，应变数值可以确定     ，松弛应力(5-1)即可定解为：

Boltzmann叠加原理说的是，线性黏弹性系统对两个独立的应变历史的响应，等于对每个应变历史响应之和。即对于    时刻施加的应变    ，和    时刻施加的应变    ，有：

那么，对于n个激励，就有：

对式(5-3)取极限，可知对于连续的应变激励，有积分形式：

或者写成更熟悉的Riemann时间积分形式：

前面Abaqus的公式里面，弛豫时间被记为    ，被积分的时间就可以写成     . 但在这里为了不至于混淆，我们的积分时间变量写为     ，其实是一码事。

(5-4b)这个形式，还可以进一步写成

是不是就和前面的式(4-3)很像了？

对于最简单的一维Maxwell模型，松弛模量     写为

那么对于Maxwell形式的标准线性模型，对应控制方程(3-3),

其最终的松弛模量可以写为：

公式(5-7)可以用和上面类似的方式推导得到。

如果是有n个弹簧-阻尼器的广义Maxwell模型，松弛模量可以用Abaqus的风格写为：

这个(5-8)是不是就和前面的(4-5)式长得很像啦。

向三维扩展的时候，参考线弹性常数之间的换算关系，对于使用     定义的1阶广义Maxwell模型，有：

这样，微分形式的弹簧-阻尼器模型(1)就与积分形式的黏弹性定义(2)统一起来了。

还剩下从(3-4)到(4-8)，从一维到三维的扩展没有搞定。

6. 一维到三维的桥梁：黏弹性本构向三维的扩展

这部分，毕小喵在阚前华老师的《非线性本构关系在ABAQUS中的实现》书中找到了解答，但阚老师可能由于笔误，这部分有一个公式应该是写漏了一项。关于这些勘误我们留到文章末尾详细讲。这一节我们直接来推导从一维到三维的扩展。

再看一遍公式(3-4)，一维Kelvin形式 标准线性黏弹性模型的微分方程：

参考弹性力学里面广义胡克定律的推导，引入泊松比：

这个写法……反正读者能理解就好。

对于x方向的应变和应变率：

x方向正应力的贡献是：

y（z同理）方向正应力的贡献就应该是：

那么，把式(3-4)扩展到三维，最终xx方向的分量就应当写为：

用张量下标形式写为：

可以看到，式(6-3)中，包括应力和应力率的迹     ,     。现在我们要把它改写成式(4-8)的形式：

就需要在关于xx的方程里，消去     和     的yy和zz项，作为代价，引入关于应变     和应变率     的yy和zz项。这部分推导和弹性力学中变换的方式是一样的。

 

在弹性力学中，用应力表示应变：

 

 

用应变表示应力：

 

整体思路是，(为了方便，这部分推导我们不再按全局编号，只使用局部公式编号），把式(6-3)展开，关于xx的等式（就是前面的式(6-2)）我们叫做(1)，关于yy的等式叫(2)，关于zz的等式叫(3)，以此类推。那么，用

消去应力的xx项，只包含对称的yy和zz项，得到（这个公式俺就不换行了）：

再用(1)式和(4)式相加，消除应力的yy和zz项：

得到最终的xx分量：

仿照上面的推导写出其余5个方向分量，再把它们合并写成张量形式，就终于又可以给一个正式编号(6-4)了：

式(6-4)就是与式(4-8)一致的形式。其中，

这样，就将一维弹簧-阻尼器形式的标准线性黏弹性模型(1)扩展到三维，完成了和Abaqus文档中三维形式黏弹性本构(3)的对应。

7. 黏弹性本构的UMAT实现与对比验证

这一节，我们参考Abaqus官方文档里的推导，介绍黏弹性本构的实现，以及与Abaqus内置黏弹性本构材料的对比验证。这里会稍微多几行我自己写的内容，但大框架还是和官方文档保持一致。

 

Abaqus | User Subroutines | Abaqus/Standard User Subroutines | UMAT | Example: Simple linear viscoelastic material

其实，黏弹性本构推导过程本身非常简单。回顾公式(4-8)：

对于这个方程，一个简单、稳定的 时间积分格式是采用中心差分。对于任意函数     ：

将(7-1)的形式代入式(4-8)，在     时刻将(4-8)的本构展开，再整体乘以     ，合并同类项得：

这里     的写法与Abaqus文档保持一致，就是三方向正应变之和，应变张量的迹。

有了 (7-2)和(7-3)式，雅可比矩阵就可以很容易地写出：

我们知道，定义一个UMAT本构，最重要的就是写出它的DDSDDE，雅可比矩阵，也称为切线刚度阵。有了这三项，对称的雅可比矩阵可以写出来，UMAT剩下的部分就是一些标准语法，没什么难度了。官方文档后面还有一些能量密度的推导，感兴趣的网友可以去翻帮助文档。

想要验证这段UMAT的黏弹性和Abaqus内置的黏弹性是否能得到相同的结果，还有一些问题需要解决。

第一个问题是，官方提供的这段UMAT代码，是不区分分析步类型的。即使是在静力，通用分析步下，它也会按照黏弹性本构来进行时间积分。这个行为和Abaqus内置的黏弹性本构就不一致。

为了解决这一问题，在UMAT的传入参数里可以找到一个叫 JSTEP() 的变量。在UMAT那一页的官方文档中是这么写的：

这个变量的第一个分量 JSTEP(1)，存储的是分析步编号。JSTEP(2)是分析步类型。

具体分析步类型对应的编号，可以在帮助文档这一页找到：

 

Abaqus | Output | File Output Format | Results file output format | Procedure type keys

 

所以，在UMAT代码里面加一个判断，就可以设置为只要分析步类型不是21，就按线弹性计算。或者更简单的方法是，分析步编号1先按线弹性计算，后面的分析步不管是啥，再按黏弹性计算也行。

第二个问题是，(4-8)形式的黏弹性本构，其参数怎么和Abaqus内置的黏弹性参数对应。

这里可以参考第6节的推导。第6节我们推的是Kelvin形式的标准线性模型到三维的扩展，使用相同的方法，也可以得到广义Maxwell形式的三维扩展。

具体来说，对于使用     来定义的广义Maxwell模型三维形式：

长期模量为     ，黏弹性部分为     ，瞬时模量就是     .那么，对于Abaqus内置黏弹性材料，有：