Abaqus 选择二阶单元计算接触的疑问
“ Abaqus使用二阶六面体单元计算接触应力，为何出现应力波动？”
本来打算今天更新下一篇 旋转量 的。虽然阅读量不高，但也是我自己的学习过程。但世事难料，今天遇到了实在比较无法理解的问题。这个问题我暂时还没得出答案，先把症状放出来与各位探讨。问题的产生
一个比较简单的赫兹接触问题，四分之一模型，圆柱与平面接触。使用一阶六面体/四面体单元，接触面附近的von-Mises应力云图如图所示。
可以看出，最大应力并非出现在接触面表面，而是出现在内部一点的地方。这和赫兹接触的解析解相一致。但下面一点的地方，由于一阶四面体属于常应变单元，应力等值线画出来就不是很光滑。
为了提高计算精度，很自然的可以有h方法和p方法。Abaqus文档和各类有限元书籍都建议不要使用一阶四面体单元，所以我仅把单元类型改为二阶。六面体单元C3D20，四面体单元C3D10，均不带任何后缀（未使用减缩积分等技术）。
奇怪的事情发生了，算出来von-Mises应力结果是这个样子的：
很明显，可以猜测是由于应力从高斯点处外推，在节点处应力过大导致的。
而且不仅仅六面体，就连下面的四面体单元也存在节点应力过大的问题，见这里的云图：应力结果后处理
咱虽然用了这么久的Abaqus，但实在是天资愚钝找了大半天没找到怎么在Abaqus里显示单元平均应力。那就用HyperView吧~ 在HyperView里面打开这两个项目的odb结果文件，Averaging method选择None，显示单元平均应力：
从这个单元平均来看，应力一点也不光滑 而且并不能让人满意。单元平均应力属于不可用的状态。而如果把云图的显示平均方式改为Simple，使用一阶单元的结果能够画出比较漂亮、连续的等值线云图。再把离散云图改为连续（取消勾选 Discrete color），出来的结果也是挺唬人的。
至于HyperView和Abaqus用相同的odb结果文件，后处理做应力平滑后给出的最大应力数值略有区别这件事，我们也都能猜到是应力平均化算法导致的，暂且搁置不论。而二阶单元，本来理论上应力计算的精度应该是显著高于一阶尤其是后面的四面体单元的，在使用了应力平均，甚至使用连续云图的时候，结果仍然很丑。
四面体部分那星星点点的深红色，丑得就像少女脸上的痘印。而且那些应力集中的区域很明显和单元划分有关，是物理上不准确的。不然怎么会一阶的结果只有200左右，二阶单元的结果却飙升到480？二阶单元提高精度不是这么提的。
简化模型计算对比
作为一位自认还算有点动手能力的选手，我自然不会坐以待毙。
我画了一个简化的圆柱-平面接触模型，来比较使用不同单元计算出的应力结果。
C3D8R（使用减缩积分的 8节点六面体单元）
C3D8S（8节点六面体单元，增强的表面应力显示）
C3D20（20节点二阶六面体单元）
C3D20R（20节点二阶六面体单元，减缩积分）
C3D20RH（20节点二阶六面体单元，减缩积分，混合单元列式）
还比较了默认接触算法（无后缀）、罚函数算法（_pene）和增广拉格朗日接触算法（_lagr），接触刚度均为默认。
结果嘛，我们只显示四分之一的圆柱体模型，看看接触面附近的应力（由于是粗网格，所以应力肯定不精确，仅做个试验）：
（咱先不谈应力应该是多少，这几个结果应力差的都有点大。但就看这儿，C3D20单元在节点处也出现了应力集中，和刚才出现的问题一样。）
至于C3D20RH单元，和不同的接触公式，云图几乎和C3D20R没区别。
应力应该是多少，这个单元密度算的肯定不准。但是抛开那个不谈，一阶C3D8单元给出的应力结果却明显要更光滑。而C3D8S单元，看起来在结果光滑的同时，计算精度也更高了不少。
反观本应更准确的20节点二阶单元，无论是否使用减缩积分，是否使用杂交公式，都同样出现了在接触面的节点处应力奇异的现象。全 军 覆 没。
Abaqus/CAE 在GUI界面上可选的二阶六面体应力单元就这么三种。帮助文档
多年养成的习惯，有困难找文档。网上搜的结果大概率没有帮助文档里来得靠谱。
Abaqus文档关于单元的部分，页面在 Elements - Continuum Elements - General-purpose continuum elements - Solid(continuum) elements。这一页比较详细地介绍了所有固体（连续介质）单元，以及各种单元技术的选择和比较。首先是减缩积分。文档中提到，二阶单元的减缩积分通常比全积分更准确。同时，二阶减缩积分单元不存在沙漏问题。
Second-order reduced-integration elements in Abaqus/Standard generally yield more accurate results than the corresponding fully integrated elements. However, for first-order elements the accuracy achieved with full versus reduced integration is largely dependent on the nature of the problem.
对于体积自锁问题，当材料为接近完全不可压缩时，二阶完全积分单元会出现体积自锁。而一阶完全积分单元由于在体积项上使用了选择性的减缩积分，因此几乎不会出现体积自锁。
当压力像棋盘一样变化，在积分点之间出现显著跳跃时，表示发生了体积自锁。
Volumetric locking occurs in fully integrated elements when the material behavior is (almost) incompressible. Spurious pressure stresses develop at the integration points, causing an element to behave too stiffly for deformations that should cause no volume changes. If materials are almost incompressible (elastic-plastic materials for which the plastic strains are incompressible), second-order, fully integrated elements start to develop volumetric locking when the plastic strains are on the order of the elastic strains. However, the first-order, fully integrated quadrilaterals and hexahedra use selectively reduced integration (reduced integration on the volumetric terms). Therefore, these elements do not lock with almost incompressible materials. Reduced-integration, second-order elements develop volumetric locking for almost incompressible materials only after significant straining occurs. In this case, volumetric locking is often accompanied by a mode that looks like hourglassing. Frequently, this problem can be avoided by refining the mesh in regions of large plastic strain.If volumetric locking is suspected, check the pressure stress at the integration points (printed output). If the pressure values show a checkerboard pattern, changing significantly from one integration point to the next, volumetric locking is occurring. Choosing a quilt-style contour plot in the Visualization module of Abaqus/CAE will show the effect.改进的表面应力可视化单元技术，避免了应力分量从积分点外推到节点产生的误差。它采用了昂贵的27积分点方案，使单元表面显示的应力更准确。同时由于它在单元节点上有积分点，因此它在单元节点处输出的应力是没有误差的。
Improved surface stress visualization bricks：
The C3D8S and C3D8HS linear brick elements have been developed to provide a superior stress visualization on the element surface by avoiding errors due to the extrapolation of stress components from the integration points to the nodes. These elements are available only in Abaqus/Standard. The C3D8S and C3D8HS elements have the same degrees of freedom and use the same element linear interpolation as C3D8 and C3D8H, respectively. These elements use a 27-point integration scheme consisting of 8 integration points at the elements' nodes, 12 integration points on the elements' edges, 6 integration points on the elements' sides, and one integration point inside the element. To reduce the size of the output database, you can request element output at the nodes. Because these elements have integration points at the nodes, there is no error associated with extrapolating integration point output variables to the nodes.
那么，根据上述文献调研结果，我也可以猜测之前出现的问题是因为发生了体积自锁。因为确实在一个单元内不同积分点的应力能差出2倍就很离谱。但是，我还是有两个问题想不通：1）体积自锁应该只发生在几乎不可压缩材料中。而我使用的就是结构钢材料，泊松比0.3，怎么会自锁呢？ 2）根据我后面简化模型的计算，即使是C3D20R减缩积分单元，在节点处依然出现了应力奇异的现象。文档中不是说二阶减缩积分单元没问题的吗？赫兹接触
在Abaqus文档的Benchmarks验证案例中，有一个 The Hertz contact problem，赫兹接触。因为赫兹给出的球面接触解析解实在是很著名。
但这一页的Benchmark案例文档中，我有两处不明。
首先 是三维圆柱体分析中，文档提到有几个可能发生接触的单元采用了C3D27单元。这个单元我是真的在GUI当中没找到。难不成需要我手动去inp文件里输入剩下的7个节点号不成？
其次，我在赫兹接触的benchmarks文档中读到了下面这样一段。他说有两个人在1980年证明，接触算法 会导致使用二阶单元时出现表面压力分布的震荡。如果我遇到的问题可以被归类为表面压力的震荡……其实也说得通。但是文档里说，这两位学者已经提出了Simpson积分规则来解决压力震荡问题，而且已经被Abaqus软件采用了啊。为什么我今天还是会遇到这个问题呢？
这一页验证案例文档，仅Abaqus/Standard就给了多达56个inp文件（简直丧心病狂），详细地比较了各种可能的情况，还使用了子结构分析。如果后续我还有时间和进一步的好奇心的话，我就挑几个官方给的inp运行一下试试看。看这个问题官方是怎么解决的。还有这个C3D27单元究竟是个啥。
后记
接触问题……和我去年最初写公 众 号时候遇到的梁问题一样，看起来、听起来都属于有限元中比较简单 比较基础的问题之一。但实际做一遍就发现，真的没那么简单啊。
今天的推送，本来想赶在12点之前完成，但实在是力有不逮，没能赶上。我不是一个遇到问题喜欢求助别人的人，我觉得我遇到的这种问题……可能一般工程师也研究不到我这么深，也解决不了。能三言两语解决我问题的工程师……恐怕大概率我也付不起他的咨询费。
即使我导师就是给研究生开有限元课的，但我想老师对这个具体问题的Abaqus实现也应该不是很熟悉。最终这个问题，要么加密网格继续使用一阶单元算出个结果，要么还是得靠我自己找到答案。 
如果有哪位老师能够告诉我，我使用的二阶单元为何会出现这样的现象，该如何解决，那就更好不过了，善哉善哉。 
