ANSYS Workbench 接触（二）

每次学完了基础操作以及相关理论之后总要做个小例子，在整个机械专业学习中似乎没有一门课专门讲解接触力学，想到也就是机械设计里面有过一个赫兹接触公式。在面对一个如此陌生的问题，如果不寻找一个有背景的例子，做完了心里也是不踏实的。赫兹接触有许多大神做过，可那毕竟是别人做的，现在我想做自己的赫兹接触。

赫兹接触

赫兹接触说的是德国科学家海因里希·鲁道夫·赫兹提出的关于弹性体接触的理论公式，关于物体之间相互作用的基本方程【百度百科】。每个公式定律都有其使用条件，在使用前务必确认使用范围，如果要延伸或者跨领域引用也应该给出合理依据。

基本假设

• 弹性体。这说的是材料的性质，为我们选择材料给与参考

• 小应变。排除了几何非线性中的大应变一类

• 光滑表面接触且不允许彼此穿透。这里提的是接触对象的表面质量，不允许穿透是实际物理现象，理想接触状态。

观察上面三个基本假设可以发现，均是对实际问题的理想化。弹性体与小应变是容易理解的，材料的本构以及几何变形影响结构刚度变化，第三条说明接触对象的表面质量是可能影响接触质量的。

典型接触

最常见的是下面两种情况：球接触与圆柱接触，基本描述如下

上图给出圆柱接触时基本力学表达式。常识告诉我，初期接触时两个圆柱之间是一条直线，当受到外力作用时，接触处应力是无穷大的。但是对于可变形体而言，其逐渐变形导致接触变成一个小而狭长的面。从上图的侧面观察到，接触压力呈现椭圆分布状态。

上图给出圆柱接触时同样适用的情况：圆柱与一个半径无穷大的平面接触、圆柱与一个圆柱形凹槽接触。

上图是两个球相互接触，接触初期是一个点，随着对象的变形，接触区域变成一个圆，接触压力亦是呈现椭圆分布。

适用情况道理类似，不赘述了。下面开始做一个球与大平面相接触的问题，具体按照如下展开。

Workbench模拟接触

使用DM建立模型，考虑到载荷与边界沿着中心轴对称，因此使用轴对称模型降低计算开销，且能够使结果满足轴对称效果。前面提过，边界载荷是对称的，其分析结果是应该含有对称的。但是使用对称分析是否适合一切结构问题，需要具体问题具体讨论。

由于与球接触的物体极大，以及小应变行为，故将大的柱体考虑为刚体，进一步缩小计算问题。

DM建模

在DM的XY平面里面建立如下的面体，且所有的几何体都落在X正半轴所在的一侧。

这里有一个问题，图中的红色框对应的选项，在图标的左下角有一个颜色不同的图标，右击查看可以看到提示，因为采用的是开放线条拉伸，这里是为了切分上面的几何体，无需理会。上半球半径为8mm，下侧由于要考虑为刚体，尺寸自己定一下，没有太严格要求。

设置材料数据仅需要弹性模量和泊松比，其余无需，两个物体采用一样的材料参数。

分析流程框图

自己搭建即可，至于是否需要采用什么样的流程图，自己确定就好。不过这里推荐下面这种方法，注意设置分析类型为2D，以前详细说过，这里就不凑字数了

设置轴对称

进入mechanical里，点击Geometry，设置详细栏如下：

设置接触

将默认的绑定接触（Bonded）修改为无摩擦接触（Frictionless），并且设置接触行为为“非对称接触”，因为打算把其中一个看成刚体，两者刚度差异巨大，如下图所示：

设置接触算法和刚度更新，算法为增广拉格朗日，刚度为每次迭代更新刚度

网格划分

设置单元为低阶

因为接触是由于压力作用产生的，其直接力学效应为变形。带有中间节点的单元，在变形较大的情况下可能导致中间节点滑移，从而影响网格质量参数：雅可比。如果网格过度扭曲可能导致无法收敛，而低阶单元则没有这种缺陷。

设置下侧柱状体的边长一个单元，因为它考虑位刚体，设置整个为一个单元，其具有很大刚性，等效为刚体。

设置上侧变形球体面的单元尺寸

之所以切分设置网格大小在于缩小计算资源，没办法，电脑硬件跟不上。所以仅对下方的进行了控制，并且要确保接触处的网格具有较好的一致性（网格要整齐）。为了使软件能够捕捉到接触对，建立良好接触，接触的网格不应该过于粗大，因为接触面需要探测几个下伏单元深度。

上侧顶点处需要施加集中力，如果划分的网格过于细，可能导致极其明显的应力奇异，而这种可能导致我们偏离了分析意图，所以为了减小此种不利情景出现，建议将上一侧的单元划分的粗一些，单元的粗细直接体现着单元的刚度。因为没有对上面进行尺寸控制，软件自动给出较粗的网格。

得到划分结果

分析设置

设置分析的载荷步

求解设置

关闭弱弹簧，确保大变形也是关闭

输出控制

分析数据管理设置

边界条件

设置下侧柱形体为固定，上侧半圆直径线X轴方向限制移动，半圆顶点受到竖直向下的集中力188.5N

查看分析结果

2D分析结果

边界上在Y向施加的是188.5N，软件得出的Y向支反力为187.76N，误差0.39%

上方基于节点加载，不可避免地出现应力高于周边的情况，下侧接触处存在挤压也是应该存在极大的应力。

初始接触的地方在竖直方向上应该是没有位移的，数值上分析结果是合理的。节点施加处应该是竖直方向位移最大的地方，上图亦是合理。

查看接触状态

初始接触处存在着滑移接触，附近是近场接触（是数值上的接触，在Pinball以内），更远处是远场接触，即不在Pinball里面。

接触压力随着接触状态发生变化，其具体数值分布可以观察下图：

自定义输出Y向应力，取绝对值，便于比较。因为实际是受压，颗粒当时负值，这里去掉符号。

可变形球体在竖直方向上最大接触应力86.086MPa

以上是以2D形式查看了基本后处理，因为是轴对称模型，这里简单看下3D状态下的扩展显示。如下图所示：

为了与理论计算做简单对比，通过球与圆柱理论公式知道，需要知晓接触半径，这里通过获取接触面积来间接得到接触半径。

据说新版本已经可以直接输出接触面积了，不知道是不是可以，本版本18.0没有找到办法，只能在网上搜集了，获得如下基本的命令流。首先需要定义一个Name Selection，对象是上半球网格控制的那条边，它在命令流里面作为一个变量：cnt_press_edges。在Solution下侧插入命令流，查看结果即可

命令流如下：

resumeset,last allsel
*get,num_el,elem,0,count *get,start_id,elem,0,num,min *dim,cnt_area,array,num_el
cmsel,s,cnt_press_edges esln,s,1 esel,r,ename,,172
*vget,elem_sel,elem,start_id,esel
allsel
*vabs,1*vmask,elem_sel *vmask,elem_sel *vscfun,elem_tot_area,sum,elem_area *vmask,elem_sel etable,carea,nmisc,27 *vmask,elem_sel *vget,cnt_area,elem,start_id,etab,carea *vmask,elem_sel *vscfun,tot_area,sum,cnt_area
my_contact_area = tot_area

获得接触面积大小为：3.237mm^2，因为接触区域形状是一个圆，可计算接触半径为1.0153mm，理论接触半径为1.00964mm，误差0.56%。

根据公式计算理论接触应力大小为88.29MPa，分析计算得到的是86.086MPa，误差约2.5%。而如果根据接触工具直接看接触压力 ，大小为88.03，这个相比较更近，但是理论计算最大值为中心处的，所以需要注意接触压力是一个总的值，不仅仅含有竖直方向的。

以上就是球与无限大平面接触的简单示例，简单走下流程，感受下接触分析。接触虽然是一件很复杂的事情，但是这个经典问题收敛确实很容易的。按照上面的设置，大概10秒就可以收敛。做出来了，感觉是不是掌握了呢。其实这里面还是存在的一些问题的。对于经典问题，已经知道其力学表现形式，应该根据力学响应形式精确的设计网格分布，如数量、位置、路径、大小、形状等，才在网格这个点上有一个基本认识。网格到目前为止，还没有一本较好的参考资料，简单的处理也就不用这里总结了，说什么梯度以及关心处做好网格也都是不具备太大实际价值，因为这样的话几乎已为人熟知，但是具体操作却不仅仅是所谓的局部精细化布置。这里只是做一个粗糙的过程，还不清楚网格是如何影响接触反应的。这里存在着不小的误差，可误差究竟是如何产生的呢，难道只是网格一个原因么？还有两个圆柱相接触，并且还有各向应力计算没有对比，是不是按照上面粗糙的过程也能满足的呢一定的精度呢？鉴于问题太多，只能慢慢学，慢慢总结了。有兴趣可以自己做一做，找了一下应力的计算公式如下：

此为球接触时的主应力和剪应力计算

此为圆柱接触时主应力与剪切应力计算

诸多问题以后再讨论学习，关于赫兹接触这一节总结到此。

注：仅记录学习FEM的一个过程，表达的是个人观点与认识，欢迎一起讨论学习。有疑问可以私，本号没有留言功能，无法互动。本人小白一枚，正在努力的路上。