经典线性粘弹性本构模型（ansys模拟）
1前言
发展到今天，有限元软件使用的单元，本构模型以及相关求解技术越来越复杂，能模拟的问题也越来越宽广。这对于软件本身来说无疑是件好事儿，但对于使用者来说却是把双刃剑，理解其中原理的人如鱼得水，不理解的很可能陷入各种云图曲线而难以抉择。个人暂时也是处于两者之间，因此常常回顾一些经典的案例以提升分析素养。今天就如何使用ansys模拟经典的线性粘弹性本构模型进行说明，作为向非线性粘弹性模型以及粘塑性模型的一个过渡。
2分析问题
参考书籍：《Elasticand Inelastic Stress Analysis》，具体如下图所示
图 1固体力学参考书籍
书中在第六章中推导了经典的线性粘弹性本构模型：
Maxwell模型
Kelvin模型
三参数标准固体模型
三参数标准流体模型
四参数的博格斯流体模型，
描述了两个基本实验
蠕变试验
松弛试验
而本文正是使用ansys模拟蠕变试验和松弛试验来研究书中本章所述各模型的力学性质。
3稍打理论
具体分析之前，我们先对上述模型有个基本的了解：
3.1 Maxwell模型
图 2 Maxwell模型示意图
Maxwell模型是一个弹簧和一个理想粘壶模型（后文简称为粘壶）串联，根据串联模型的特点我们可以得到如下方程：
根据上述式子我们可以得到Maxwell本构模型为：
3.2 Kelvin模型
图 3 Kelvin模型示意图
3.3 其它经典模型
图 4其它经典粘弹性模型
上述模型从左到右依次是：
标准固体模型：弹簧+Kelvin模型
标准流体模型：粘壶+Klevin模型
博格斯流体模型：Kelvin模型+Maxwell模型
这些模型的具体推导可以参照教程，这里不讲太多理论。
4问题分析
通过上述模型，大致可以判断出线性粘弹性模型仅包含两个基本要素：线性弹簧，线性粘壶（所谓线性就是弹簧的弹性系数以及粘壶的粘性系数保持常数）。对应到ansys中的单元，能直接找到弹簧阻尼单元combine14，我们先看下弹簧阻尼模型：
图 5弹簧阻尼单元说明
看起来弹簧阻尼单元（2D轴向弹簧阻尼）和Kelvin模型一样，但是有点区别得注意，材料本构表述的是应力应变之间的关系，而弹簧阻尼表述的是力和位移之间的关系，其实个人理解起来两者是一样，因为K=EA，如果想用弹簧的K代替前面的E，则只需将A看成单位1即可，对于应变我们也可以将长度理解成单位1，这样在弹簧阻尼模型中，力和应力，位移和应变实际是一个等效的概念（后面输出曲线就直接使用应力和应变的说法）。
对于分析类型，由于引入了粘壶，因此需要使用瞬态动力学进行分析以考虑粘性效果，否则粘壶将失效。
5Maxwell模型的蠕变以及松弛实验
选择单元类型为combine14，更改关键字K3选择为2D轴向弹簧，并且为弹簧单元创建两个实常数，一个仅包含弹性系数k=1，另一个仅包含阻尼系数c=1。直接使用节点（0,0,0），（0,1,0），（0,2,0）创建对应的单元，得到如下Maxwell模型（D0为阻尼，K0为弹簧）：
图 6 Maxwell模型的有限元模型
5.1 蠕变试验（creep test）
按照蠕变实验要求需要对整体施加一个瞬时的总应力，随后保持总应力不变，研究模型总应变与时间的关系。为了让蠕变的性质能体现的更加明显，这里我采用三个载荷步去模拟：
载荷步1：底端节点固定，对上端节点在1e-3的时间（较短时间）内施加大小为1的力载荷，载荷形式为斜坡加载
载荷步2：外载保持不变，分析2s
载荷步3：卸载，分析1s
注意：整个分析过程采用瞬态分析模拟，并且需要打开时间积分效应；虽然变形很大，但是由于各元素在全局坐标系下的刚度并不因构型改变而变化，因此并不属于大变形，不必打开瞬态大变形开关）
得到Maxwell模型的变形云图以及应变-时间曲线图如下：
7 Maxwell模型蠕变试验动图
图 8 Maxwell模型典型蠕变曲线
上图为典型的Maxwell模型的蠕变曲线
初始弹性变形阶段：可以看到随着瞬时应力的施加，弹簧产生了瞬时弹性应变，这时粘壶刚度无限大，不产生应变（载荷步1）
蠕变变形阶段：随着应力的持续作用，即使大保持不变，但是粘壶正在不断产生应变，这里称为蠕变应变（载荷步2）
卸载阶段：卸载之后，弹性应变立即消失，但是蠕变应变不可恢复（载荷步3）
5.2 松弛试验（relaxation test）
与蠕变实验不同，松弛试验是施加一个瞬时的总应变，随后保持应变不变，研究模型应力与时间的关系。因此按照要求，我们分两个载荷步去模拟松弛实验：
载荷步1：底端节点固定，对上端节点在1e-3的时间（较短时间）内施加大小为1的位移载荷，载荷形式为斜坡载荷
载荷步2：外载保持不变，分析10s
得到Maxwell模型的变形云图以及应力-时间曲线图如下：
图 9 Maxwell模型松弛试验动图
图 10 Maxwell模型典型松弛曲线
上图为典型的Maxwell模型的松弛曲线：
初始应力阶段：由于瞬时应变的施加，弹簧产生瞬时的应力
松弛阶段：随着时间的累积，弹簧的弹性应变逐渐转变为粘壶的应变，体系整理的应力水平不断降低，对于Maxwell模型会直接降为0，意味着所有的弹性应变都会转换为粘壶的不可恢复应变。
6博格斯模型的蠕变以及松弛试验
上面使用Maxwell模型作为引例，简单介绍了下蠕变实验以及松弛试验，下面就复杂点的博格斯模型为例，说明下复杂的线性粘弹性现象如何去模拟（个人认为是比较直观的建模方法）。
回顾第三部分的博格斯模型，我们知道是将一个Kelvin模型与一个Maxwell模型串联起来的，而且一个弹簧阻尼单元本身相当于一个Kelvin模型，但是为了扩展可视性，准备将弹簧与阻尼拆成两个，使用如下方法进行建模：
图 11博格斯模型的有限元模型
其中自由度耦合可以看成连接上下部分的枢纽，至于为什么不将弹簧阻尼单元合并，前面也有说明，主要是为了分析的可拓展性和可视性，这样建模对于更加复杂的模型也能处理，而且能很清晰看到结构组成。
6.1 蠕变试验
对上述模型进行第5部分所述相同的操作，我们可以得到如下结果：
图 12博格斯模型蠕变试验动图
图 13博格斯模型典型蠕变曲线
可以看到，博格斯模型与Maxwell模型的区别主要有以下两点：
Maxwell模型蠕变阶段呈现线性，而博格斯模型蠕变阶段呈现非线性
Maxwell模型卸载后弹性应变回复，应变不再变化，而博格斯模型卸载后，弹性应变回复，但是后续仍然发生蠕变应变的变化。
其实不管是上述的非线性还是后续的蠕变变化都是因为Kelvin模型的特点导致，感兴趣的读者可以研究下他的性质，这里不作说明。
6.2 松弛试验
同样，我们对博格斯模型进行第5部分的松弛试验，可以得到以下结果：
图 14博格斯模型松弛试验动图
图 15博格斯模型典型松弛曲线
图 16耦合点位移曲线
由博格斯模型的松弛试验，我们得到了和Maxwell模型一致的松弛曲线，但是两者松弛速率是有区别的，主要体现在前期Maxwell模型松弛速率低于博格斯模型，后期松弛速率又高于博格斯模型。
这里附上耦合点处的位移曲线是为了让大家感觉下位于模型中一点的运动情况，可以看到耦合点的位移从初期到后期呈现先迅速上升到最高点，后迅速下降，随后缓慢下降，这个过程脑补下动图，其实可以感觉到应变在组件内流动，也是弹性应变，蠕变应变在进行相互转化，进行后续粘塑性学习可能体会更深刻，这里就简单点一点。
7时间
最后，专门就时间这个概念进行下说明，因为以前也会困惑，使用材料库中的粘性本构进行分析时，算出的时间到底是1s，1h，1年的蠕变曲线？貌似在分析中都没有体现时间单位这个概念，那么我怎么判断自己的时间单位？以Maxwell模型为例，具有两个基本参数，就是弹性系数和阻尼系数，其中弹性系数毫无疑问与时间无关，但是阻尼系数=力/速度，所以其单位为N*S/m，很明显，如果给的阻尼系数是建立在s单位制下的，那么我的分析结果也是s，如果给的想以h作为分析单位，那么本例中的阻尼系数都应该改乘以3600，这样得出的结果也是在h单位下的，那么对应到软件中的材料本构，我们也可以得到类似的结论，因此在给定粘性本构的时候，一定要确定其中的参数所建立的单位制。
8小结
本文只是就经典粘弹性模型中的一部分进行了详细说明，但是根据上述方法很容易对更加复杂的粘弹性本构进行模拟，希望通过文章的阅读，读者能体会到粘弹性的一些性质，这也是理解更加复杂的具有粘性性质本构的基础。
