有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变
（原创，转载请注明出处）
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图1]()
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图2]()
==概述==
本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式。有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。
一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元求解器，通过自研求解器和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。
iSolver介绍视频：
http://www.jishulink.com/college/video/c12884
==第18篇：几何非线性的应变 ==
上一章我们介绍了几何非线性的物理含义，这章接着介绍几何非线性中重要的物理量：应变。物体在受到外力作用下会产生一定的变形，变形的程度称应变。理解应变的物理含义对有限元结果的判别和试验的对应有重要意义。线性情况下的应变非常容易计算，但如果存在几何非线性，就不是所有人都知道Abaqus等商软结果中的应变E代表的是什么了。本章首先从位移、变形和应变的区别说起，然后再具体介绍一下几何非线性下的应变是如何度量的。
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图8]()
1.1 位移 Displacement
一个物体从初始状态C0受到外部载荷运动到另一个状态C，由于材料不能凭空消失，那么C0下的任意的一个物质点a（设坐标为X）必然在C下有个对应点A（设坐标为x）。位移表示两个位置坐标相减，为u=x-X。
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图11]()
1.2 变形 Deformation
位移无法表示出变形情况，譬如大黄蜂变形为汽车时，身体360度转动，位移很大，但没有任何的变形。
变形必须要由两个点之间材料纤维的拉伸才能决定。对初始构想C0任意的两点a、b，变形后设为为AB，那么AB矢量和ab之间的关系可以由
dx=F*dX
来表示。
其中F定义为deformation gradient matrix。
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图16]()
1.3 应变 Strain
显然F与X和x的坐标系相关，dX和dx都是在特定三维空间笛卡尔坐标系下的矢量值。它是个张量，有9个量。F在理论上可以表示物体受力后的变形情况，但实际中试验没法测出这么多分量，因为拉伸试验等测试的几何量只能是位移，长度等有限的几个标量，这些基本量是无法得到F的全部分量的。所以只能寻找新的量来表示物体变形情况。
一种很自然的想法就是取模，也就是取F的绝对值，称为拉伸率strech ratio：
r=||dx||/||dX||
橡胶等材料，变形比较大，r可以达到2以上，那么在试验上很容易测量，但对金属等变形比较小的情况，一倍的力和两倍的力得到的r只是1.0001和1.0002，差异非常小，依然很难标定。
所以，为了测量和计算方便，人为假定了另一个表示物体变形情况的变量，就是应变，因为是人为设定的，所以应变的取法有很多种，主要考虑两个条件的便捷性：
（条件a）在试验上可以简单的表示为位移，长度等基本量，且工况不同时，差异较大，能容易反应出物体加载情况。
（条件b）在理论上可以方便地和位移、应力一起形成最终正确的求解方程。
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图21]()
1.4 应变的度量选择
为了简单起见，我们只针对一维情况来说明应变度量选择，应变的一般公式可以表示为r的一个函数：
Strain=f(r)
由于应变是反应物体变形的量，那么当不变形时，也就是lambda=1时，显然Strain=0。
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图22]()
1.4.1 工程应变
当变形特别小时，也就是r~=1时，假定应变和位移是线性的，那么取伸长量和原始长度的比值。
也就是：
Strain=u/X=(x-X)/X =r-1
x和X分别为最终和最初的长度。一倍的力和两倍的力得到的应变变为了0.0001和0.0002，差两倍关系，很容易区分出来。
此时，显然条件a很容易满足，同时，在理论上将应变和位移看成是线性关系去求刚度阵，受力后得到的位移也容易计算，且结果和试验一致，因此，条件b也满足。
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图25]()
1.4.2 真实应变
当存在几何大变形的情况，变形比较大时，此时x-X已经和X比拟了。
此时工程应变中的应变和位移是线性关系这种假定求出的位移已经和实际偏差较大了，因此，需要采用新的应变度量方式。一种选择是应变取从初始时刻到最终时刻点这段时间的累加，显然，这种选择就与材料点在整个时间段内的变形路径有关。
问题：那么从上图的X到最终的x，我们取哪个路径呢？
答：取真实的路径。
实际情况是什么路径，就应该取什么路径，此时在真实路径下得到的应变就是真实应变。当然真实路径不一定总是直接从X直线到x点的，譬如下方
(1) 如果受力是FA，那么的确是X直线到x点A路径；
(2) 但如果是受力FB，那么可能是路径B，只不过试验上总是采用FA来做拉伸实验罢了。
在路径A下，可以得到
也就是真实应变和拉伸率的对数相等，这也是为什么真实应变也称为对数或者自然应变的原因。
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图32]()
1.4.3 Green应变
真实应变在一维下很容易计算，但是在三维情况下，Strain和F相关，上面的积分就没法直接求了，只能得到一个类似
这种表达式求积分，时间的度量由增量来表示，真实应变由应变增量累加得到。因此，在计算机中每次求真实应变必须花费大量时间。
注：实际中有限元的几何非线性求应变增量需要计算局部模态，而模态分析的求解时间成本很大。
一般情况只能老老实实积分计算，但对大位移大转动小应变的特殊情况，发现应变如果取为Green应变，计算精度不受影响，同时，三维的应变也可以简单的表示出来。所以对小应变情况，可以采用Green应变，一维形式如下：
Strain=0.5（r^2-1）
![有限元理论基础及Abaqus内部实现方式研究系列18： 几何非线性的应变的图37]()
1.5 壳的应变度量
如果是壳单元，那么用的是局部坐标系，假设为L坐标系，初始构型下存在坐标系变换dX=F2*dL。类似的，几何大变形中任意一个构形也有类似的变换。此时变形量F或者真实应变就和L坐标系相关了。
具体的也可看一下第十四篇：壳的应力方向或者下面配合iSolver求解器的视频解说
https://www.jishulink.com/college/video/c12884 20.14 理论系列文章14：壳的应力方向
==总结==
本章首先从位移、变形和应变的区别说起，然后通过一维的简单例子具体介绍了几何非线性下的应变的度量方式。
至于线性的应变在商软或者自主软件中在单元中实现，可以参考下面的视频：
https://www.jishulink.com/college/video/c14948
深入浅出有限元：基础理论->Abaqus操作->matlab编程
如果有任何其它疑问或者项目合作意向，也欢迎联系我们：
snowwave02 From www.jishulink.com
email: snowwave02@qq.com
以往的系列文章：
第一篇：S4壳单元刚度矩阵研究。介绍Abaqus的S4刚度矩阵在普通厚壳理论上的修正。
http://www.jishulink.com/content/post/338859
第二篇：S4壳单元质量矩阵研究。介绍Abaqus的S4和Nastran的Quad4单元的质量矩阵。
http://www.jishulink.com/content/post/343905
第三篇：S4壳单元的剪切自锁和沙漏控制。介绍Abaqus的S4单元如何来消除剪切自锁以及S4R如何来抑制沙漏的。
http://www.jishulink.com/content/post/350865
第四篇：非线性问题的求解。介绍Abaqus在非线性分析中采用的数值计算的求解方法。
http://www.jishulink.com/content/post/360565
第五篇：单元正确性验证。介绍有限元单元正确性的验证方法，通过多个实例比较自研结构求解器程序iSolver与Abaqus的分析结果，从而说明整个正确性验证的过程和iSolver结果的正确性。
https://www.jishulink.com/content/post/373743
第六篇：General梁单元的刚度矩阵。介绍梁单元的基础理论和Abaqus中General梁单元的刚度矩阵的修正方式，采用这些修正方式可以得到和Abaqus梁单元完全一致的刚度矩阵。
https://www.jishulink.com/content/post/403932
第七篇：C3D8六面体单元的刚度矩阵。介绍六面体单元的基础理论和Abaqus中C3D8R六面体单元的刚度矩阵的修正方式，采用这些修正方式可以得到和Abaqus六面体单元完全一致的刚度矩阵。
https://www.jishulink.com/content/post/430177
第八篇：UMAT用户子程序开发步骤。介绍基于Fortran和Matlab两种方式的Abaqus的UMAT的开发步骤，对比发现开发步骤基本相同，同时采用Matlab更加高效和灵活。
https://www.jishulink.com/content/post/432848
第九篇：编写线性UMAT Step By Step。介绍基于Matlab线性零基础，从零开始Step by Step的UMAT的编写和调试方法，帮助初学者UMAT入门。
http://www.jishulink.com/content/post/440874
第十篇：耦合约束（Coupling constraints）的研究。介绍Abaqus中耦合约束的原理，并使用两个简单算例加以验证。
https://www.jishulink.com/content/post/531029
第十一篇：自主CAE开发实战经验第一阶段总结。介绍了iSolver开发以来的阶段性总结，从整体角度上介绍一下自主CAE的一些实战经验，包括开发时间预估、框架设计、编程语言选择、测试、未来发展方向等。
http://www.jishulink.com/content/post/532475
第十二篇：几何梁单元的刚度矩阵。研究了Abaqus中几何梁的B31单元的刚度矩阵的求解方式，以L梁为例，介绍General梁用到的面积、惯性矩、扭转常数等参数在几何梁中是如何通过几何形状求得的，根据这些参数，可以得到和Abaqus完全一致的刚度矩阵，从而对只有几何梁组成的任意模型一般都能得到Abaqus完全一致的分析结果，并用一个简单的算例验证了该想法。
http://www.jishulink.com/content/post/534362
第十三篇：显式和隐式的区别。介绍了显式和隐式的特点，并给出一个数学算例，分别利用前向欧拉和后向欧拉求解，以求直观表现显式和隐式在求解过程中的差异，以及增量步长对求解结果的影响。
http://www.jishulink.com/content/post/537154
第十四篇：壳的应力方向。简单介绍了一下数学上张量和Abaqus中壳的应力方向，并说明Abaqus这么选取的意义，最后通过自编程序iSolver来验证壳的应力方向的正确性。
https://www.jishulink.com/content/post/1189260
第十五篇：壳的剪切应力。介绍了壳单元中实际的和板壳近似理论中的剪切应力，也简单猜测了一下Abaqus的内部实现流程，最后通过一个算例来验算Abaqus中的真实的剪切应力。
https://www.jishulink.com/content/post/1189260
第十六篇：Part、Instance与Assembly。介绍了Part、Instance与Assembly三者之间的关系，分析了Instance的网格形成原理，并猜测Abaqus的内部组装实现流程，随后针对某手机整机多part算例，通过自编程序iSolver的结果比对验证我们的猜想。
https://www.jishulink.com/content/post/1195061
第十七篇：几何非线性的物理含义。介绍了几何非线性的简单的物理含义，并通过几何非线性的悬臂梁Abaqus和iSolver的小应变情况的结果，从直观上理解几何非线性和线性的差异。
