有限元理论基础及Abaqus内部实现方式研究系列7：C3D8R六面体单元的刚度矩阵

==概述==

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

一方面我们查阅Abaqus软件手册得到修正方法的说明，另一方面我们自己编程实现简单的结构有限元求解器，通过自研求解器和Abaqus的结果比较结合理论手册如同管中窥豹一般来研究Abaqus的修正方法，从而猜测商用有限元软件的内部计算方法。在研究的同时，准备将自己的研究成果记录下来写成一个系列文章，希望对那些不仅仅满足使用软件，而想了解软件内部实现方法甚至是做自己的软件的朋友有些帮助。由于水平有限，里面可能有许多错误，欢迎交流讨论。

==第七篇：C3D8R六面体单元的刚度矩阵==

前面几篇文章都是介绍了梁和壳的结构，其实梁和壳只是体的两种特殊形式，一般情况下体才是一个结构正常的形状，而且当计算机足够强大时，仅仅只要用体单元就能求解结构问题，这个时候甚至可以不用做几何简化和中面提取等辅助操作了，大大减少了前处理的时间。

前面文章提到梁相对壳来说，商业软件的修正方式相对较少，如果自己编程序，采用这些修正方式可以得到和商业软件完全一致的梁单元刚度矩阵，如果刚度矩阵完全一致，那么对任何的梁的算例都可以得到和商业软件完全一致的结果了。进一步，对体来说，商用软件对体的修正方式相对梁更少，因为梁对轴向拉伸、轴向扭转、横向弯曲、横向剪切等的刚度都需要不同的修正，但体只需要一个方向的修正就行，只要掌握了体的修正方式，那么自编程序的体单元组成的任意结构分析都可以得到商业软件完全一致没有误差的结果。

在Abaqus的Standard和Explicit求解中默认的体单元都为C3D8R，C表示体单元，3D表示三维单元，8表示八节点，R表示Reduced减缩积分，也就是八节点一次线性六面体单元，该单元在Abaqus中是计算耗时最少的六面体单元，而六面体的计算精度本身就比四面体和锲形单元高，所以C3D8R在体单元实际应用中最广泛，我们本章以此为研究对象。

本文先讨论一般的六面体单元的基本理论和C3D8R的修正方式，然后在自编有限元程序iSolver实现同样的修正方式，最后验证iSolver的结果和Abaqus完全一致，从而证明Abaqus的内部修正和我们设想的一致。具体验证过程也可以参考我们的演示录像：

https://www.jishulink.com/college/video/c12884 章节3

1.1 六面体单元刚度矩阵的基本理论

体单元刚度矩阵的理论相对梁和壳来说要简单很多。按照有限元的一般的等参理论：

（1）单元的节点上的坐标只要网格划分完后就已经确定了，那么单元内的任意一点的坐标可以表示为节点上的坐标值。

i=1到节点N，在六面体单元中一次和二次单元节点个数一般为8、20，四面体分别为4、10。20节点六面体的另外12个节点正好位于12条边的中点，10节点四面体类似。

Ni表示形函数，类似于一个比重，所求的x点约靠近xi点必然这个比重越大，最大显然应该是x点位于xi点上，此时比重就是1，其它的形函数为0。

（2）有了坐标的表达式，由等参的定义，位移的表示和坐标完全一致，那么就可以得到单元内任意一点位移的表示：

（3）位移对坐标的微分就得到应变：

• 注意Abaqus的剪切应变顺序为12，13，23，和这边的B矩阵顺序不同。

（4）由输入的材料得到应力与应变的关系，即本构关系

（5）体的刚度矩阵就是：

1.2 Abaqus的修正方式

一个[-1,1]内连续的积分此时就转换为了离散的积分，积分点的数目由gi的个数决定。最简单的情况就是取一个积分点，就是C3D8R单元，表示减缩积分。Abaqus的C3D8R除了减缩积分外，还做了下面的修正。

只要是减缩积分都会存在沙漏，C3D8R也不例外。

1.2.1 沙漏产生的原因

单元受到弯曲作用时，正常情况下形状应该出现如下图所示：

上方材料纤维受拉，下方纤维受压，导致上方是拉力，下方是压缩力，这样使得在左右边界力矩保持平衡，如下图所示。

对线性减缩单元，将出现如下情况：

由于单元拉伸应力只计算中心点P的值，而中心点P在X方向的纤维L1的长度没有变化，使得单元拉伸应力为0，也就是和外力矩M无法平衡，导致在M作用下单元扭曲的非常厉害。

一个单元扭曲的形状如上图，当存在多个这样的单元时，会形成如下的形状。上下两个单元就类似一个沙漏。

显然，沙漏只存在于弯曲现象中，对拉伸没有此现象。

从刚度的数学表达式也可以解释这个沙漏现象。

在减缩积分时，只有一个积分点，那么上述值就是

上述应变与位移的关系B0矩阵为中心点P的值，存在多个项是相同的，使得刚度阵K也存在多个相同值的元素，导致刚度矩阵秩的缺失，举个简单的例子，就容易导致类似下方的刚度矩阵。

    

2X2的矩阵满秩的话秩是2，但上面这样的矩阵秩只有1，在加了载荷后是无解的，只不过在计算机中由于精度误差，计算中的刚度矩阵将会变为类似下方的矩阵：

这时只要加一个非常小的载荷将导致解非常大。这也是当存在沙漏时，单元变形明显变大的原因。

1.2.2 控制沙漏的算法

沙漏控制现在最常用的方法是增加一个人工的沙漏刚度，譬如：

                                                                        [0.001, 0

                                                                        0, 0.001]

增加这个沙漏刚度单元刚度矩阵将变为下方的矩阵：

                                                                        [ 1.001 ，1

                                                                        2， 2.001]

这个时候解依然存在，不一定是正确解，但肯定比由于计算机精度造成的误差小很多了。

那么具体如何增加沙漏刚度和大小应该取多少合适呢？

这个沙漏控制应该只对沙漏模式（对应的位移函数设为Γ）起作用，但对其它模式（对应的位移函数为O）不起作用，也就是在没有沙漏现象的正常情况这个刚度应该不起作用，这点通过增加一个归一化的γ矩阵，使得γ*Γ=1，同时对其它模式γ*O=0来实现。

在一次减缩积分的六面体单元中，以x方向为例，存在四个沙漏模式

对应的沙漏位移函数为：

后面就是怎么取这个γ的算法了，与Γ和应变应力矩阵Bi有关，实际值为：

和普通的刚度矩阵一样，沙漏刚度的值为

V是体积，Dhg为沙漏的本构关系，Factor来控制沙漏刚度的大小。

1.2.3 Abaqus的沙漏控制方法

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）              沙漏刚度会导致产生人工的沙漏应变能，该能量在实际中并不存在，这点要求这个因子越小越好。

上面两点本身就是矛盾的两点，所以很多时候这个因子根本不存在，当出现沙漏的时候，修改沙漏因子并不能真的解决沙漏问题。

 

1.2.4 控制沙漏的实际意义

在前面章节S4R的沙漏控制也提到了，通过算法来控制减缩单元譬如C3D8R/S4R的沙漏，这个沙漏因子很多情况下没法同时满足不增加太多的人工刚度又控制了沙漏，但Abaqus还是默认在C3D8R和S4R中增加了一个小的沙漏因子，既然这个小的沙漏因子对存在的沙漏无法得到有效控制，那么意义在哪里？

可能是笔者有误，对Abaqus了解的还不够深入，也希望跟大家一起探讨，只能猜测一种可能是这种沙漏对Explicit分析及其重要，毕竟Abaqus的Explicit对沙漏的控制要比Standard有更多的选择，另一种可能情况是不在于这个因子的大小，而只在于是否有这个因子，如果没有沙漏，那么只要涉及到刚度矩阵求逆的操作都直接失败，而存在这个增加的人工沙漏后可以让刚度矩阵的秩增加，从而所有矩阵求逆的操作都能顺利进行。以壳单元S4R为例，譬如下面的这个例子：

计算一个四周简支板的模态分析，采用减缩壳单元S4R。

第一种情况：如果没有沙漏因子，Abaqus就会计算失败。

采用iSolver自研程序计算，可得到下面的前10阶模态，和理论值相比，可以发现多了前三项不存在的伪模态。猜测这可能也是Abaqus计算失败提示的内部意义。

第二种情况：当加入沙漏因子后。

Abaqus和iSolver都可以计算出前10阶模态精确值，说明沙漏的存在使得模态分析求解顺利进行。

1.3 算例验证

Abaqus的修正方式仅仅是猜测，下面我们在自研求解器iSolver中对八节点六面体单元按照上面猜测的修正方法，和理论和Abaqus比较从而验证我们的猜测。

1.3.1 分片试验

1.3.1.1 模型描述

模型：该算例取之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单元：

输入材料和载荷如下：

1.3.1.2 分析结果

1.3.1.2.1 理论结果

1.3.1.2.2 iSolver结果

S11和S12的结果和理论完全一致，证明iSolver的六面体单元通过基本的分片试验。

     

1.3.2 弯曲梁的算例

1.3.2.1 模型描述

模型：该算例取之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。

1.3.2.2 分析结果

1.3.2.2.1 理论结果

由论文上可得理论上U3最大位移为0.5022。

1.3.2.2.2 Abaqus结果

Abaqus分析后得到U3最大值=0.5332，相比理论值为1.06：

1.3.2.2.3 iSolver结果

Abaqus分析后得到U3最大值=0.5332，相比理论值为1.06，和Abaqus完全一致：

1.3.2.2.4 结果分析

Abaqus和iSolver的C3D8R的位移结果都比理论值偏大一些，是因为弯曲载荷下存在沙漏作用，沙漏控制只是减弱了位移的幅度，但刚度还是比实际的要偏软。如果采用C3D8完全积分，可得到下面的位移，明显比理论值偏小，此时没有沙漏现象，但存在剪切自锁使得单元偏刚。

1.3.3 任意六面体组成的结构

对任意六面体单元组成的结构，iSolver的结果都和Abaqus完全一致，具体可看演示录像：

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

==总结==

本文首先简单介绍了一下8节点六面体单元的基本理论，分析了刚度矩阵的来源，并研究了Abaqus中应用最广泛的C3D8R单元的刚度矩阵的修正方式，自编程序如果采用这些修正方式可以得到和Abaqus完全一致的分析结果，并用多个算例验证了该想法。

Abaqus的C3D8R六面体单元的刚度矩阵在八节点线性单元理论基础上的修正如下表：

项次	问题	修正情况	说明
修正	不修正
1	积分	√		采用减缩积分
2	沙漏现象	√		增加人工的沙漏刚度对沙漏进行控制

有兴趣的可以自行下载iSolver进行验证，因为看不到Abaqus的源代码，上述C3D8R的修正方式也仅是猜测，如果你在使用iSolver测试其它的由C3D8R组成的模型结果时发现和Abaqus结果不一致，欢迎联系我们。

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

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