黏弹性本构：从一维模型到三维的UMAT数值实现
0. 开篇
（其实这里标题应该叫“引言”？）但我不想按论文的风格写的那么严肃。
出于最近的科研需要，小喵完整学习了一下关于黏弹性本构的知识。小喵获取的这些知识散落在许多来源，很多课本里的推导要么不完整，要么甚至有错误。所以我学的也是磕磕绊绊。或许之前有老师或者教材完整讲过，但至少我在调研学习的时候没找到。本着造福后人的精神，希望其他对黏弹性相关理论感兴趣的朋友，在搜索到、读过这篇文章以后，学习过程能稍微轻松一点。
这篇文章，和过去小喵写的许多推送一样，也是一篇我自己的学习笔记。略有不同的是，之前学习笔记性质的推送，很多与仿真软件操作有关的内容我是边学边写整理出来的。而这一篇，由于大部分是理论和公式，我是在自己完全看懂以后，再梳理出正确的逻辑一气呵成写下来的。这样对读者更友好。
这篇文章只在开头有一个简单的黏弹性Abaqus仿真案例。更多还是以公式推导为主，从一维讲到三维，一直讲到线性黏弹性本构的UMAT实现。我使用了mdnice墨滴编辑器的Markdown语法和LaTeX公式，让文章的公式排版看起来更清爽。
顺带一提，幸好最新版本的Word同时支持Unicode和LaTeX格式的公式，我可以在Word和OneNote里用简洁的Unicode语法敲公式，然后把它转成LaTeX放进Markdown。有些我在OneNote里推过的公式，就不需要再用LaTeX语法重新敲一遍了。
小喵是很菜的博士生，在推导公式方面更是没什么天分。这篇文章的公式推导我会尽可能讲的慢一点，详细一点，以匹配我这种水平选手的理解能力。希望对这个主题感兴趣的读者也都能够读懂。
正文目录：
“黏”还是“粘”？
Abaqus黏弹性案例
一维黏弹性模型
Abaqus中的黏弹性参数
微分与积分形式的桥梁：Boltzmann叠加原理
一维到三维的桥梁：黏弹性本构向三维的扩展
黏弹性本构的UMAT实现与对比验证
书籍资料勘误及参考文献
在正文内容安排的顺序上，
首先用一个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内置黏弹性材料，有：
而对于(4-8)形式的UMAT，其参数的对应关系为：
注意到这里 对应的是长期模量。想要获得瞬时模量，需要：
有了这些关系，我们就可以进行对比验证了。
验证案例取
对应的Abaqus内置材料参数为：
| E_inst | v | k_i | g_i | tau_i |
|---|---|---|---|---|
| 200 | 0.3 | 0.5 | 0.5 | 1 |
UMAT参数为：
| lambda | G | lambda1 | G1 | nu1 |
|---|---|---|---|---|
| 57.69 | 38.46 | 115.385 | 76.923 | 1 |
使用尺寸为1*1*1的单个单元模型作拉伸验证，内置材料和UMAT材料在静态分析步位移同为 5.293e-2，黏性分析步位移为 7.537e-2，结果一致，验证通过。
附 - 结合了线弹性和黏弹性的UMAT代码：
C=====================================================C
C 公 众号 《CAE知识地图》 C
C 作者 毕小喵 C
C=====================================================C
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),
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),
5 JSTEP(4)
DIMENSION DSDEL(NTENS,NTENS),DSDEC(NTENS,NTENS),STREND(NTENS),
1 PSTRAN(3),ASTRAN(3,3),barCp(3,3),Q(6,6),
2 DSTRES(6),D(3,3)
C
C
Elambda = PROPS(1)
Emu = PROPS(2)
Elambda1 = PROPS(3)
Emu1 = PROPS(4)
Enu1 = PROPS(5)
C
C determine step
C
IF(JSTEP(2).NE.21) Then
ElamInst = Elambda1*Enu1
EmuInst = Emu1*Enu1
EBULK3 = 3*ElamInst+2*EmuInst
EG2 = 2*EmuInst
EG = EmuInst
EG3 = 3*EmuInst
ELAM = ElamInst
C
DSDEL=0.0
Do 10,K1=1,NDI
Do 20,K2=1,NDI
DSDEL(K1,K2)=ELAM
20 continue
DSDEL(K1,K1)=ELAM+EG2
10 continue
c
Do K1=NDI+1,NTENS
DSDEL(K1,K1)=EG
ENDDO
c
Do K1=1,NTENS
STREND(K1)=STRAN(K1)+DSTRAN(K1)
ENDDO
c
DDSDDE=DSDEL
STRESS=MATMUL(DDSDDE,STREND)
c
ELSE
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 + Enu1
TERM1I = 1./TERM1
TERM2 = (.5*DTIME*Elambda+Elambda1)*TERM1I*DEV
TERM3 = (DTIME*Emu+2.*Emu1)*TERM1I
C
DO K1=1,NDI
DSTRES(K1) = TERM2+TERM3*DSTRAN(K1)
1 +DTIME*TERM1I*(Elambda*EV
2 +2.*Emu*STRAN(K1)-STRESS(K1))
STRESS(K1) = STRESS(K1) + DSTRES(K1)
END DO
C
TERM2 = (.5*DTIME*Emu + Emu1)*TERM1I
I1 = NDI
DO K1=1,NSHR
I1 = I1+1
DSTRES(I1) = TERM2*DSTRAN(I1)+
1 DTIME*TERM1I*(Emu*STRAN(I1)-STRESS(I1))
STRESS(I1) = STRESS(I1)+DSTRES(I1)
END DO
C
C CREATE NEW JACOBIAN
C
TERM2 = (DTIME*(.5*Elambda+Emu)+Elambda1+
1 2.*Emu1)*TERM1I
TERM3 = (.5*DTIME*Elambda+Elambda1)*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 = (.5*DTIME*Emu+Emu1)*TERM1I
I1 = NDI
DO K1=1,NSHR
I1 = I1+1
DDSDDE(I1,I1) = TERM2
END DO
C
C TOTAL CHANGE IN SPECIFIC ENERGY
C
TDE = 0.
DO K1=1,NTENS
TDE = TDE + (STRESS(K1)-.5*DSTRES(K1))*DSTRAN(K1)
END DO
C
C CHANGE IN SPECIFIC ELASTIC STRAIN ENERGY
C
TERM1 = Elambda + 2.*Emu
DO K1=1,NDI
D(K1,K1) = TERM1
END DO
DO K1=2,NDI
N2 = K1-1
DO K2=1,N2
D(K1,K2) = Elambda
D(K2,K1) = Elambda
END DO
END DO
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
I1 = NDI
DO K1=1,NSHR
I1 = I1+1
DEE = DEE + Emu*(STRAN(I1)+.5*DSTRAN(I1))*DSTRAN(I1)
END DO
SSE = SSE + DEE
SCD = SCD + TDE - DEE
ENDIF
RETURN
END
8. 书籍资料勘误及参考文献
首先是前文提到的，阚前华老师的这本《非线性本构关系在ABAQUS中的实现》。
在这本书的37页，给出了一维Kelvin形式标准线性模型对应的三维形式，并在下方注明了参数之间的对应关系。然而，书中 (3-11)式是错的，漏掉了应力率的迹 这一项。按照书中 (3-11)的形式是无法推出(3-12)的。
这本书后面从(3-13)开始的公式都是正确的，这里的笔误无伤大雅。这些参数之间的对应关系仍很有价值。
——————————
然后，是一份在网上传的比较广泛的《各向同性粘弹性umat理论讲解.pdf》，它还配套了一段黏弹性UMAT子程序。
原PDF文档是英文。如作者所说，这是一个Kelvin-Voigt体的黏弹性本构。事实上，它的参数除了 和 ，就只有一个黏度 。
这个本构的三维形式完整写出来，和(4-8)的形式相比，彻底缺少了 这一项。如果仅用来学习参考……倒也不是不行，但反正不太好用。想要用广义Maxwell模型模拟它，只需要把黏壶对应的那根弹簧的刚度设置的非常大即可。但反过来，它是无法模拟广义Maxwell模型的。从参数的数量上就不对。
本文撰写过程中，其余的参考资料还包括：
[1] 维 基百科-Viscoelasticity. https://en.wiki pedia.org/wiki/Viscoelasticity
[2] 冯元桢教授的教材。我参考的是很厚的英文书 Classical and computational solid mechanics，但读者去翻翻更简单的中文版，《连续介质力学初级教程》其实也行。
[3] 来自University of Glasgow, UK 的 René de Borst教授的一本英文书，NON-LINEAR FINITE ELEMENT ANALYSIS OF SOLIDS AND STRUCTURES - SECOND EDITION. 这本书的第8章，Time-dependent Material Models里面有Boltzmann叠加原理，和许多更深入的推导。
[4] B站\仿真秀平台上 九千CAE 老师的黏弹性理论视频和Abaqus-Fortran子程序直播视频。
[5] CAE中学生公 众号 图惜 大佬 写的推送 Ansys Workbench工程应用之——结构非线性（中）：材料非线性（4）粘弹篇
——————————
本文相关的资料文件，有：
Abaqus验证案例，黏弹性垫圈的inp文件 visco_cae.inp； 线弹性-黏弹性结合的UMAT子程序及验证inp: OneElem_UMAT.inp； 冯元桢教授的英文书 PDF； Borst教授的英文书 PDF； 阚前华老师的《非线性本构关系在ABAQUS中的实现》 PDF（转自CAE中学生）； 这篇推送的markdown源代码，包含全部LaTeX格式公式。
