测试功能：无摩擦接触；

参考文献：

[1] K.L. Johnson, Contact Mechanics, Journal of tribology, 1985, 108(4):464.

[2] 某商业软件 Benchmarks Guide.

问题描述

几何：

半径10，厚度1。

材料：

弹性模量 E=3e7，泊松比v=0.3。

边界条件：

壳约束全部自由度，实体约束z方向自由度（平面应变），x=0表面约束Ux=0。

荷载：

上表面沿y负方向施加位移Uy=-0.4。

接触：

壳面选为主面，圆柱面选为从面，无摩擦接触。

单元：

Hex8，Tet10，Hex20。

参考解

由赫兹公式得出，见参考文献[1]。

结果对比

三种单元得到的最大等效应力分别为：

  某商业软件手册中的最大Mises应力可能并不是网格收敛的解（因其云图中的网格显示较为粗糙，其加密网格后应力可能会相应的降低），笔者用Simdroid在网格相对粗糙时得到的等效应力与某商业软件更接近，随着网格的逐渐加密，最终得到表中Hex8和Tet10单元的结果。而Hex20单元存在典型的二阶serendipity单元接触问题（因其缺乏面中点导致），应力结果不可信，而某商业软件采用了变节点C3D27单元，包含面中点，所以不存在该问题。

变形后的接触压力和x坐标曲线如下图：

  作为对比，某商业软件手册中关于本问题的结果为：

结果云图

Hex8单元的等效应力云图为：

Tet10单元的等效应力云图为：

Hex20单元的等效应力云图为：

某商业软件手册[2]中关于本问题的云图为：

结论

  各个单元的接触压力在靠近接触中心处的计算结果均与Hertz接触公式结果相一致，在远离接触中心处略有偏差。一阶六面体单元Hex8和二阶四面体单元Tet10的等效应力云图结果基本与某商业软件结果一致，Hex20单元的应力分布存在较大的偏差。导致该问题的原因主要为接触力在Hex20单元的表面上分配时采用的serendipity单元二阶形函数，会使得边中点处出现负数（即接触压力为负），从而出现局部的应力集中。