产品
==概述==
本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式。有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。
一方面我们查阅Abaqus软件手册得到修正方法的说明，另一方面我们自己编程实现简单的结构有限元求解器，通过自研求解器和Abaqus的结果比较结合理论手册如同管中窥豹一般来研究Abaqus的修正方法，从而猜测商用有限元软件的内部计算方法。在研究的同时，准备将自己的研究成果记录下来写成一个系列文章，希望对那些不仅仅满足使用软件，而想了解软件内部实现方法甚至是做自己的软件的朋友有些帮助。
进一步，对体来说，商用软件对体的修正方式相对梁更少，因为梁对轴向拉伸、轴向扭转、横向弯曲、横向剪切等的刚度都需要不同的修正，但体只需要一个方向的修正就行，只要掌握了体的修正方式，那么自编程序的体单元组成的任意结构分析都可以得到商业软件完全一致没有误差的结果。
在Abaqus的Standard和Explicit求解中默认的体单元都为C3D8R，C表示体单元，3D表示三维单元，8表示八节点，R表示Reduced减缩积分，也就是八节点一次线性六面体单元，该单元在Abaqus中是计算耗时最少的六面体单元，而六面体的计算精度本身就比四面体和锲形单元高，所以C3D8R在体单元实际应用中最广泛，我们本章以此为研究对象。
本文先讨论一般的六面体单元的基本理论和C3D8R的修正方式，然后在自编有限元程序iSolver实现同样的修正方式，最后验证iSolver的结果和Abaqus完全一致，从而证明Abaqus的内部修正和我们设想的一致。
体单元刚度矩阵的理论相对梁和壳来说要简单很多。按照有限元的一般的等参理论：
1）单元的节点上的坐标只要网格划分完后就已经确定了，那么单元内的任意一点的坐标可以表示为节点上的坐标值。
i=1到节点N，在六面体单元中一次和二次单元节点个数一般为8、20，四面体分别为4、10。20节点六面体的另外12个节点正好位于12条边的中点，10节点四面体类似。
Ni表示形函数，类似于一个比重，所求的x点约靠近xi点必然这个比重越大，最大显然应该是x点位于xi点上，此时比重就是1，其它的形函数为0。
2）有了坐标的表达式，由等参的定义，位移的表示和坐标完全一致，那么就可以得到单元内任意一点位移的表示：
3）位移对坐标的微分就得到应变：
4）由输入的材料得到应力与应变的关系，即本构关系
5）体的刚度矩阵就是：
一个[-1,1]内连续的积分此时就转换为了离散的积分，积分点的数目由gi的个数决定。最简单的情况就是取一个积分点，就是C3D8R单元，表示减缩积分。Abaqus的C3D8R除了减缩积分外，还做了下面的修正。
只要是减缩积分都会存在沙漏，C3D8R也不例外。
单元受到弯曲作用时，正常情况下形状应该出现如下图所示：
上方材料纤维受拉，下方纤维受压，导致上方是拉力，下方是压缩力，这样使得在左右边界力矩保持平衡，如下图所示。
对线性减缩单元，将出现如下情况：
由于单元拉伸应力只计算中心点P的值，而中心点P在X方向的纤维L1的长度没有变化，使得单元拉伸应力为0，也就是和外力矩M无法平衡，导致在M作用下单元扭曲的非常厉害。
一个单元扭曲的形状如上图，当存在多个这样的单元时，会形成如下的形状。上下两个单元就类似一个沙漏。
显然，沙漏只存在于弯曲现象中，对拉伸没有此现象。
从刚度的数学表达式也可以解释这个沙漏现象。
在减缩积分时，只有一个积分点，那么上述值就是
上述应变与位移的关系B0矩阵为中心点P的值，存在多个项是相同的，使得刚度阵K也存在多个相同值的元素，导致刚度矩阵秩的缺失，举个简单的例子，就容易导致类似下方的刚度矩阵。
2X2的矩阵满秩的话秩是2，但上面这样的矩阵秩只有1，在加了载荷后是无解的，只不过在计算机中由于精度误差，计算中的刚度矩阵将会变为类似下方的矩阵：
这时只要加一个非常小的载荷将导致解非常大。这也是当存在沙漏时，单元变形明显变大的原因。
沙漏控制现在最常用的方法是增加一个人工的沙漏刚度，譬如：
[0.001, 0
0, 0.001]
增加这个沙漏刚度单元刚度矩阵将变为下方的矩阵：
[ 1.001 ，1
2， 2.001]
这个时候解依然存在，不一定是正确解，但肯定比由于计算机精度造成的误差小很多了。
那么具体如何增加沙漏刚度和大小应该取多少合适呢？
这个沙漏控制应该只对沙漏模式（对应的位移函数设为Γ）起作用，但对其它模式（对应的位移函数为O）不起作用，也就是在没有沙漏现象的正常情况这个刚度应该不起作用，这点通过增加一个归一化的γ矩阵，使得γ*Γ=1，同时对其它模式γ*O=0来实现。
在一次减缩积分的六面体单元中，以x方向为例，存在四个沙漏模式
对应的沙漏位移函数为：
后面就是怎么取这个γ的算法了，与Γ和应变应力矩阵Bi有关，实际值为：
和普通的刚度矩阵一样，沙漏刚度的值为
V是体积，Dhg为沙漏的本构关系，Factor来控制沙漏刚度的大小。
Abaqus也是通过增加一个沙漏刚度来控制C3D8R的沙漏：
即在Element Type中可以设置Hourglass Control的方法，Standard中default就是stiffness的控制方法，只有Explicit其它几项才有作用，同时设置了一个缩放因子Scaling factors来控制沙漏刚度的幅度。
这个沙漏刚度在物理上没有任何的物理意义，是人为加进去的，所以由它计算出来的任何量是没有意义的，这也是abaqus后处理中ALLAE这个变量中的AE表示Artificial Strain Energy的含义。因为没有意义，想要结果正确，那么这个人工的应变能相对总的内能ALLIE就需要越小越好。一般认为ALLAE/ALLIE的比值必须控制在5%-10%以下，才认为计算结果可信。
沙漏的因子Factor很多时候需要修改，它的取值受下面两者影响：
1）有效抑制沙漏模式，这点要求这个因子越大越好。
2）沙漏刚度会导致产生人工的沙漏应变能，该能量在实际中并不存在，这点要求这个因子越小越好。
上面两点本身就是矛盾的两点，所以很多时候这个因子根本不存在，当出现沙漏的时候，修改沙漏因子并不能真的解决沙漏问题。
在前面章节S4R的沙漏控制也提到了，通过算法来控制减缩单元譬如C3D8R/S4R的沙漏，这个沙漏因子很多情况下没法同时满足不增加太多的人工刚度又控制了沙漏，但Abaqus还是默认在C3D8R和S4R中增加了一个小的沙漏因子，既然这个小的沙漏因子对存在的沙漏无法得到有效控制，那么意义在哪里？
可能是笔者有误，对Abaqus了解的还不够深入，也希望跟大家一起探讨，只能猜测一种可能是这种沙漏对Explicit分析及其重要，毕竟Abaqus的Explicit对沙漏的控制要比Standard有更多的选择，另一种可能情况是不在于这个因子的大小，而只在于是否有这个因子，如果没有沙漏，那么只要涉及到刚度矩阵求逆的操作都直接失败，而存在这个增加的人工沙漏后可以让刚度矩阵的秩增加，从而所有矩阵求逆的操作都能顺利进行。以壳单元S4R为例，譬如下面的这个例子：
计算一个四周简支板的模态分析，采用减缩壳单元S4R。
第一种情况：如果没有沙漏因子，Abaqus就会计算失败。
采用iSolver自研程序计算，可得到下面的前10阶模态，和理论值相比，可以发现多了前三项不存在的伪模态。猜测这可能也是Abaqus计算失败提示的内部意义。
第二种情况：当加入沙漏因子后。
Abaqus和iSolver都可以计算出前10阶模态精确值，说明沙漏的存在使得模态分析求解顺利进行。
Abaqus的修正方式仅仅是猜测，下面我们在自研求解器iSolver中对八节点六面体单元按照上面猜测的修正方法，和理论和Abaqus比较从而验证我们的猜测。
模型：该算例取之Abaqus Verification Manual 6.12-1的1.5.2 Patch test for three-dimensional solid elements，是Abaqus用来证明体单元的分片试验。inp文件为：Job-PatchTest-Solid-7Hexa.inp，我们只取第一个step。
模型由7个不规则的斜六面体组成，内部的四个节点任意，结果都是一样的。在Abaqus中建模和网格划分，采用C3D8R单元：
输入材料和载荷如下：
S11和S12的结果和理论完全一致，证明iSolver的六面体单元通过基本的分片试验。
模型：该算例取之MacNeal论文[2]，为当初验证Nastran体单元的一个例子。
inp文件为Job-CurvedBeam-Z.inp。
模型为一个1/4的半圆形的梁，一端固支，另一端受面外拉力1N，分析受力点在面外的位移。
几何：内径4.12，外径4.32，厚度h=0.1。
材料：杨氏模量1e7，泊松比0.25。
网格：整体划分0.05大小的结构化网格，同时单独设置厚度方向为4个网格，C3D8R单元。
边界和载荷：一端固支，另一端受面外拉力1N。
由论文上可得理论上U3最大位移为0.5022。
Abaqus分析后得到U3最大值=0.5332，相比理论值为1.06：
Abaqus分析后得到U3最大值=0.5332，相比理论值为1.06，和Abaqus完全一致：
Abaqus和iSolver的C3D8R的位移结果都比理论值偏大一些，是因为弯曲载荷下存在沙漏作用，沙漏控制只是减弱了位移的幅度，但刚度还是比实际的要偏软。如果采用C3D8完全积分，可得到下面的位移，明显比理论值偏小，此时没有沙漏现象，但存在剪切自锁使得单元偏刚。
对任意六面体单元组成的结构，iSolver的结果都和Abaqus完全一致，具体可看演示录像：
https://www.jishulink.com/college/video/c12884
本文首先简单介绍了一下8节点六面体单元的基本理论，分析了刚度矩阵的来源，并研究了Abaqus中应用最广泛的C3D8R单元的刚度矩阵的修正方式，自编程序如果采用这些修正方式可以得到和Abaqus完全一致的分析结果，并用多个算例验证了该想法。
Abaqus的C3D8R六面体单元的刚度矩阵在八节点线性单元理论基础上的修正如下表：
| 项次 | 问题 | 修正情况 | 说明 | |
| 修正 | 不修正 | |||
| 1 | 积分 | √ | 采用减缩积分 | |
| 2 | 沙漏现象 | √ | 增加人工的沙漏刚度对沙漏进行控制 |
有兴趣的可以自行下载iSolver进行验证
免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删
155-2731-8020
