有限元理论基础及Abaqus内部实现方式研究系列27： Abaqus内部计算和显示的应变

SnowWave02 浏览：2820 评论：25 收藏：5

（原创，转载请注明出处）

==概述==

有限元理论基础及Abaqus内部实现方式研究系列27： Abaqus内部计算和显示的应变的图3 本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过

（1）基础理论

（2）商软操作

（3）自编程序

三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。

有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

有限元理论基础及Abaqus内部实现方式研究系列27： Abaqus内部计算和显示的应变的图4

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元求解器，通过自研求解器和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。

自主结构有限元求解器iSolver介绍视频：

有限元理论基础及Abaqus内部实现方式研究系列27： Abaqus内部计算和显示的应变的图5

http://www.jishulink.com/college/video/c12884

==第27篇：Abaqus内部计算和显示的应变 ==

应变度量在线性情况下没有差异，但在几何非线性情况下有不同的度量方式。系列文章18：几何非线性的应变中，我们提到了几何非线性常用的三种应变：工程应变、Green应变、对数应变，其中工程应变和Green应变分别用于线性和小应变情况，对数应变用于几何大应变情况，那么在Abaqus的几何大应变中，是否真的就简单的采用了对数应变呢？Abaqus后处理显示的应变和计算的应变是一样的吗？我们将举壳单元和体单元两个例子通过iSolver和Abaqus的结果对比来验证几何大应变Abauqs内部计算用的应变和后处理显示的应变。前面有朋友也问几何非线性时iSolver计算用的应变是什么？趁此文章我们也将证明iSolver计算时不同的单元用不同的应变，每种应变的选取方式和Abaqus一致。

1.1 对数应变的计算方法

当存在几何大变形的情况，变形比较大时，此时终止时刻位移x-初始时刻X已经和X比拟了。

对数应变取从初始时刻到最终时刻点这段时间的累加，与材料点在整个时间段内的变形路径，是一种积分行为。

（1）在一维情况下，设拉伸率为r：

这个等式其实隐含了应变是沿x方向的，如果在三维空间可以认为有一个主方向n=（1,0,0）’，因为应变是二维张量，那么上式就是：

（2）在三维情况下，存在三个主方向n1,n2,n3，拉伸率分别为r1，r2，r3：

剩下的问题就是怎么求这三个主方向，这个主要是求一个与变形F相关的3X3矩阵的特征值和特征向量，具体可参考其它论文书，其实和Abaqus后处理中显示应变时有三个Principal的值是一样的求法。

1.2 变形率积分的计算方法

由上面计算方法发现每次计算对数应变都需要求一个特征值和特征向量，在数值计算中，特征问题的求解耗时相对较多，且计算相对复杂（一般人都是认为计算复杂才采用别的应变，个人不太认可），而实际许多非线性材料中，都有这样一个规律，就是弹性应变都相对较小，譬如典型的钢材料，杨氏模量为2.1e11Pa，屈服应力为235Mpa，那么达到屈服时的应变为235e6/2.1e11=1e-3，同时，典型的应力应变本构曲线如下图，在塑性段譬如C点的弹性应变和屈服应变差异并不大。

因此，Abaqus中假定内置的所有材料都满足弹性应变相对较小，此时，理论可以证明，对数应变可以简单的取为变形率D的积分：

上式无法直接得到数学表达式，但在有限元中，可以通过增量形式累加。

1.3 调试Abaqus内部计算应变的方法

由于对数应变和试验最接近，因此Abaqus后处理中的E都是用对数应变来显示的，Abaqus为了进一步提示对数应变，直接在后处理中如果选了NLGeom=On，应变的显示从E变味了LE，但Abaqus几何非线性实际计算的应变并不完全一致。

Abaqus真正计算的变量度量可以通过它的子程序的输入参数获取。在Abaqus中，增量步即代表时刻点，可以查看增量点时刻的子程序输入来猜测Abaqus的内部量描述方式。UMAT子程序中，在材料本构函数中要利用应变增量和当前应力等物理量更新应力，查看UMAT等子程序的接口：

可知其中STRAN和DSTRAN分别表示当前增量步最后时刻的应变全量和增量。具体的介绍也可参考下面视频讲解：

http://www.jishulink.com/college/video/c13034

变形率D在一维上代表对数应变的导数，但三维上并不是对数应变的导数，这是有很大区别的，同时，可以利用iSolver分别采用上述两种应变度量和Abaqus子程序接口的结果比对来确认Abaqus计算的应变是哪种度量。所以下面我们将找一个体单元和一个壳单元的例子来验证到底Abaqus计算和显示的应变是什么。

1.4 体单元的例子

1.4.1 算例介绍

体单元算例参数如下：

尺寸：5X1X0.1。

材料：Young’s Modulus 1e8， Poisson Ratio 0.3。

左侧四个节点固支。

右侧四个节点约束位移为5,1,1。

划分为一个壳单元C3D8R。

几何非线性开关NLGeom=On，且控制只迭代一次。

1.4.2 Abaqus的应变

Abaqus中采用壳的UMAT子程序进行计算。

（1）显示应变：Abaqus计算完毕后得到导入结果，在后处理中查看，应变E11=0.6819,E22=5.621e-3如下：

（2）计算应变：Abaqus中采用UMAT子程序，利用我们的子程序调试插件DUS调试UMAT，在Visual Studio中查看dStran的值，发现在计算完应变后，进入UMAT时，显示如下，可得E11=0.6667， E22=0：

1.4.3 iSolver的应变

猜测体单元中，abaqus显示的是对数应变，而计算还是用的变形率。为了证明这点，在iSolver以同样方式建模，只是iSolver中采用自带材料进行计算，材料参数和UMAT的输入完全一致。

计算中，iSolver也采用变形率计算方式，得到的应变显示E11=0.6667， E22=0，如下，可发现和Abaqus完全一致。

1.5 壳单元的例子

1.5.1 算例介绍

壳单元算例参数如下：

尺寸：5X1，厚度0.1。

材料：Young’s Modulus 1e8， Poisson Ratio 0.3。

左侧两个节点固支。

右侧两个节点每个加集中力6.73128E+006，x方向。

划分为一个壳单元S4R。

几何非线性开关NLGeom=On，且控制只迭代一次。

1.5.2 Abaqus的应变

Abaqus中采用壳的UMAT子程序进行计算。

（1）显示应变：Abaqus计算完毕后得到导入结果，在后处理中查看，应变E11=8.528e-1,E22=-5.173e-1如下：

（2）计算应变：Abaqus中采用UMAT子程序，利用我们的子程序调试插件DUS调试UMAT，在Visual Studio中查看dStran的值，发现在计算完应变后，进入UMAT时，E11=8.528e-1,E22=-5.173e-1，调试如下：

可以发现壳单元Abaqus的计算应变和显示应变一样，猜测都是对数应变。

1.5.3 iSolver的应变

iSolver中采用自带材料进行计算，材料参数和UMAT的输入完全一致。

为了计算和Abaqus完全一致，iSolver也采用对数应变计算方式，得到的应变显示如下，可发现和Abaqus完全一致。

==总结==

由上可以看到，在实际计算中，对体单元，Abaqus和iSolver都采用变形率积分方式来计算应变，对壳单元，Abaqus和iSolver都采用对数应变。一般理论书都认为Abaqus是因为对数应变计算复杂才采用别的应变，但个人认为应该不是这个原因，因为Abaqus对体单元为了显示对数应变，依然重新计算了一遍，说明Abaqus体单元采用变形率是有其它原因的，具体什么原因我也没研究清楚，欢迎探讨。

如果有任何其它疑问或者项目合作意向，也欢迎联系我们：

snowwave02 From www.jishulink.com

email: snowwave02@qq.com

以往的系列文章：