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

==概述==

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

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

iSolver介绍视频：

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

==第17篇：几何非线性的物理含义 ==

几何非线性，即考虑大变形对于结构平衡位置的影响。爆炸冲击、冲压成型、大型结构件的弯曲等都含有几何非线性问题，几何非线性也是现代结构有限元商业软件的必备发展方向。在Abaqus中只要简单的在Step中勾选NL Geom这个开关就行。

K.J.Bathe教授1979年提出的几何非线性理论也是目前应用于有限元分析最广泛的几何非线性力学。但要想在自主结构有限元程序中编程加入几何非线性的理论，远远不是Abaqus或者Ansys表面上的看起来只要加一个NLGeom=on/off这个开发那么容易。同时，要让自研程序的几何非线性做到和商软结果接近远比线性或者材料非线性难，我们在iSolver编写几何非线性的过程中也发现，除了刚度矩阵的修改，增量迭代法的自动步长选取，收敛准则等都有极大的影响。从本章开始，将介绍几何非线性的一些理论和Abaqus的实现方式，同时通过iSolver的程序验证Abaqus的实现方式。本章将介绍几何非线性的简单的物理含义，并通过几何非线性的悬臂梁Abaqus和iSolver的小应变情况的结果，从直观上理解几何非线性和线性的差异。配合本章的视频解说和操作演示可看下方：

http://www.jishulink.com/college/video/c12884 20理论系列文章17-几何非线性的物理含义

2.1 几何非线性的物理含义

2.1.1 从线性到几何非线性

一个物体从初始状态A由于受到外部载荷运动，如果现在已知了另一种状态B的位移，那么其它的任意状态C的位移怎么求？

如果能直接从B的位移乘以一个常量就得到C，那么这个系统就是线性系统。譬如下面的800mm的悬臂梁问题，在Abaqus中用线性计算，载荷F和位移u是直线关系。

载荷1N的时候Abaqus计算得到最大位移时1.177mm。

那么载荷是1000N时是多少？显然，不用计算也知道就是1177mm。 但1177这个值明显有问题，已经超过了梁的长度，按生活经验判断这个梁估计都断了或者极端扭曲了，所以这种情况需要用几何非线性来计算。

2.1.2 几何非线性简单物理含义

虚功原理如下：

在物理上可解释能量守恒原理，即在某一个时刻点，假定在外力作用下有个虚拟的位移，那么外力在虚拟位移下做的虚功=内部应变能的变化相同。

为了更好的理解上面的物理解释，如果我们把当成真实的位移，那么外层加上对时间的积分，可以理解为外力在虚拟位移下做的虚功=内部应变能的一段小时间内对应变能的积分：

那么就和我们高中时学过的小球受重力作用后势能、动能相互转换是一样的原理。

举个简单的有限元例子，譬如下方一个减缩积分S4R单元的右端受X方向两个力F：

得到的位移为U（相当于虚位移），那么力F做的功是W=2F*U。

增加的应变能为S*E*V在对所有时刻点t的积分，S、E、V都是当前时刻的值。

应变E的取法有很多种，采用真实应变，那么E取为位移U(t)和长度L(t)的比值，按虚位移的定义，虚位移必然相对原始长度比较小，也就是L(t)=L0+U(t)可以用L0代替，E=U(t)/L0，如果是线性系统，U(t)=U*t/总时间，积分很容易计算出来，得到应变能V=S*U*b*h，因为内力和外力平衡，减缩积分S4R面积内的所有点的应力和中心点一样，所以S=F/截面积=2F/(b*h)，此时V=2F*U=W。

如果是非线性系统，那么应变就没法简单的用E=U(t)/L0，W随t的变化就是个非线性过程。每个时刻点可以求出一个斜率，这个斜率最终会形成当前时刻点的刚度矩阵。

如果是对当前时刻的体积积分，那么对W求导就很困难，因为V也是与时间有关的，可以选择一个不变的初始构型V0，此时应力和应变也需要做相应的变化，我们假定分别变为了S和E。

也就是刚度矩阵将分为两块，上式的前面一部分依然是以前的BDB形式，只不过B换成了当前时刻的应变位移矩阵，而后面新增项一般称为几何刚度阵，在Abaqus中称为初始应力矩阵（initial stress stiffness）。

2.1.3 几何非线性的计算机求解方式

理论上受力曲线是一条光滑曲线，计算机没法求解曲线上每个时刻点的结果，只能求解部分有限间隔点的结果。非线性问题不是一条直线，所以需要多次迭代才能实现，而不再考虑u=F/K这种一次性就能计算的简单问题。非线性问题求解有多种方法主要分为以下几类：增量法、迭代法、增量迭代法和弧长法。

具体的理论和Abaqus实现过程可参考我们以前系列文章：第四篇：非线性问题的求解。介绍Abaqus在非线性分析中采用的数值计算的求解方法。

http://www.jishulink.com/content/post/360565

2.2 悬臂梁的几何非线性算例

此次验证，依然使用自研求解器iSolver与Abaqus计算结果对比的方式。

2.2.1 算例说明

我们采用上面提到的悬臂梁的例子，我们采用壳单元来模拟整个梁。

参数如下：

尺寸：X方向长度L=800，Y方向宽度b=20，Z方向厚度h=20。

材料：Young’s Modulus E=11000。

点载荷：F=1000。

网格：在X方向划分20份，Y方向划分8份（为了避免Abaqus和iSolver中都会出现的沙漏现象，厚度方向尽量多划分网格）。

Step中设置NLGeom打开，并且步长=0.25。

• 我们发现就算不设置这个步长，采用自动步长，Abaqus也依然会自动调整为0.25的步长，为了和iSolver对比，我们都设置0.25，避免因为收敛判据的不同造成的结果差异。

如果线性情况下，上面已经得到了就是1177，现在几何非线性情况，我们采用小应变的几何非线性。

• 注：这个问题虽然位移很大，但依然还是个小应变问题，小应变单元的计算速度远远高于大应变，所以我们只用小应变单元来模拟。

2.2.2 Abaqus结果

Abaqus中设置单元为S4R5。

得到最后一个增量步的位移为628.4：

同时可以发现应变为0.03，和1相比是小量，所以可以用小应变来模拟。

2.2.3 iSolver结果

在运行结束后，得到的位移为627.9，和Abaqus相差0.08%。

2.3 遗留问题

Abaqus应力判据为：

其中，

为最大应力残差；

为误差因子，默认为5e-3。

q为平均时间应力，即.msg中的Average Force值，在Abaqus的帮助文档中有具体的计算公式，但我们按照公式第一个增量步的结果和Abaqus完全一致，后面的后面的增量步差异很大，这是导致我们的迭代次数和Abaqus没法完全一致的原因，有知道q第二增量步怎么计算的大神还请不吝赐教。

==总结==

本文介绍了几何非线性的简单的物理含义，并通过几何非线性的悬臂梁Abaqus和iSolver的小应变情况的结果，从直观上理解几何非线性和线性的差异。

视频解说及算例操作演示如下：

http://www.jishulink.com/college/video/c12884 20理论系列文章17-几何非线性的物理含义

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

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