浅谈ANSYS Workbench接触设置

浅谈ANSYS Workbench接触设置
0、引言
ANSYS中的接触可涉及位移、电压、温度、磁场等自由度，在这些接触中，涉及位移自由度的接触是比较复杂的。

本文大概介绍了ANSYS中接触求解的原理，并使用ANSYS Workbench计算了两圆柱接触和轮齿接触的接触应力并与赫兹公式进行了对比，最后给使用ANSYS Workbench求解接触时提供了一些建议。

鉴于作者水平有限，难免会存在一些错误，希望广大读者批评指正。

1、ANSYS接触公式理论
接触处理往往是复杂的。

可能的话推荐使用程序默认的设置。

因为现实接触体之间不会相互穿透，程序必须在两个表面之间建立一种关系，在分析中阻止彼此穿透。

程序阻止相互穿透的行为被称之为强制“接触兼容性”。

图1 接触穿透示意图
为了在接触界面上强制执行兼容性，Workbench Mechanical通常提供了几个接触公式。

这些公式定义了使用的求解方法。

图2 接触算法设置界面
•纯罚函数法
•增广拉格朗日法
•常规拉格朗日法
•多点约束（MPC）法
•梁(beam)
如果穿透在一个接触容差（FTOLN*下层单元的深度）范围内，接触兼容性则是满足的。

接触深度是一个接触对中每个接触单元深度的平均值。

如果程序检测到任意穿透大于这个容差，全局求解仍然认为是不收敛的，即使残余力和位移增量达到了收敛准则。

图3 下层单元深度示意图
2、纯罚函数法和增广拉格朗日法接触公式
对于非线性实体接触面，可使用纯罚函数公式或者增广拉格朗日法公式。

这两个都是基于罚函数接触公式：
F Normal=K Normal*X Penetration
有限接触力F Normal，是接触刚度K Normal的函数。

接触刚度越高，接触穿透X Penetration越小，如下图说明：
图4 接触刚度与接触穿透的示意图
理想的，对于一个无限大的接触刚度K Normal，可以获得一个0穿透。

在基于罚函数方法下这在数值上是不可能的，但是，如果只要X Penetration足够小或者可以忽略，则认为求解结果是精确的。

纯罚函数法和增广拉格朗日法之间最大的区别是增广拉格朗日法补充了接触力(pressure）的计算：
纯罚函数公式：F Normal=K Normal*X Penetration；
增广拉格朗日公式：F Normal=K Normal*X Penetration+λ；
因为有额外项λ，所以相较于纯罚函数法，增广拉格朗日法对接触刚度K Normal的大小不太敏感。

3、常规拉格朗日接触公式
常规拉格朗日法公式增加了一个额外的自由度(接触压力)来满足接触兼容性。

而不是将接触力分解成接触刚度和穿透值，而是接触力（接触压力）被作为一个额外的自由度来求解。

F Normal=DOF。

•用压力自由度来强制0或接近0的穿透；
•不需要一个法向接触刚度（0弹性滑移）;
•需要直接求解器，会增加计算的要求。

4、常规拉格朗日震颤
震颤问题经常发生在常规拉格朗日法中。

如果不允许0穿透（下图左半部分），接触状态就在开放或关闭（一个阶跃函数）之间变化。

这通常使得收敛更加困难，因为接触点会在一个开放或关闭状态之间振动，这被称之为“震颤”。

如果允许一些轻微的穿透（下图右半部分），就会使得求解容易收敛，因为接触不再是一个阶跃变化。

图5 接触状态示意图
5、多点约束(MPC)接触公式
对于面之间是确定的绑定或者未分离接触类型，多点约束公式是可用的。

多点约束内部增加约束方程来将接触表面之间的位移“绑定”起来。

这种方法不再是基于罚函数法或者基于拉格朗日乘子法。

是一个直接且有效的方法来联系绑定接触区域的表面。

基于多点约束法的绑定接触支持大变形效应。

6、梁接触公式
梁接触公式仅对于绑定类型是可用的。

这个公式是使用无质量的线性梁单元来将接触拓扑“缝合“在一起。

7、公式的比较
8、增广拉格朗算法
增广拉格朗日算法是罚函数算法的一系列更新迭代，去寻找拉格朗日乘子（即，接触牵引）。

增广拉格朗日算法比罚函数法通常会带来更好的矩阵状态并且对接触刚度系数的大小不太敏感。

然而，在一些分析中，增广拉格朗日算法需要额外的迭代，特别是变形的网格变的更加扭曲，接触压力定义如下：
ε—兼容性容差（以定义实常数指令R 输入的FTOLN ） λi —在第i 次迭代时的拉格朗日乘子分量。

拉格朗日乘子分量λi 是局部计算的（对于每个单元）并且是迭代的。

通过上式可知，接触压力与接触刚度（这里仅指法向接触刚度）和接触穿透有关，默认的接触刚度与材料参数、单元大小和用户定义的穿透容差有关。

如果用接触刚度因子（默认是1.0，范围是0.1-10）来调整接触刚度的大小，在材料和单元大小一定的情况下，接触刚度只与穿透容差有关，默认下接触刚度与最终的穿透容差成反比，穿透容差越小，接触刚度越大。

接触穿透容差也可以由因子和实际值来定义（ANSYS 中正值为因子，其范围一般是小于0.2，默认是0.1，负值表示实际容差值，长度单位）。

这给接触压力的调整指明了方向，网格一定情况下，若想增大接触压力，可以调小接触穿透容差或者增加接触刚度；其他因素一定下，可以调整网格，使接触附着单元深度减小。

ANSYS 官方文件中给出了ANSYS 计算的接触压力与赫兹公式计算的赫兹接触应力值的比较，使用的赫兹公式为：
σH =√
F n (1ρ1±1ρ2
)
πL(1−μ12E 1+1−μ2
2E 2
)
σH —接触应力（Mpa ）
F n —法向接触力（N ） L —接触线长度（mm ）
ρ1—圆柱1接触半径（mm ）； ρ2—圆柱2接触半径（mm ）； E 1—圆柱1弹性模量（Mpa ）； E 2—圆柱2弹性模量（Mpa ）； μ1 —圆柱1泊松比； μ2 —圆柱2泊松比；
π—圆周率。

这是两个圆柱接触的接触应力计算公式。

9、验证部分
9.1 两圆柱接触分析
但是需要注意，使用线性单元在关键区域不要使用这种单元的退化形式。

ANSYS在验证接触应力时也是使用的线性单元，所以本文使用线性单元solid185。

图6 接触设置界面
在验证的过程中通过加密网格，接触应力趋于一个稳定的值，但是与赫兹公式计算的结果相差甚远，因此调整网格并不能来解决这个问题，通过改变接触刚度和接触穿透容差得到了接触应力与赫兹公式计算的结果相一致。

为了避免偶然性，只改变材料的参数(弹性模量)得到一些列结果，表3，这些结果与赫兹公式计算的结果一致，说明接触的设置有效。

从赫
兹接触公式中可以看到，接触应力与弹性模量的关系是非线性的，而ANSYS计算不同弹性模量的接触压力与赫兹公式计算的接触应力一致，两种结果算法不同，说明若针对某一接触问题，只要得到合适的接触刚度和接触容差，那么求解接触问题是完全可行可信的。

图7 验证模型网格图
图8 弹性模量为2E2Mpa时的接触应力云图
图9 等效应力云图
图10 位移矢量和
9.2 轮齿接触应力
半径，α是压力角，本文使用的是渐开线圆柱直齿轮，压力角是20°，通过赫兹接触公式计算分度圆处的接触应力为1776.15Mpa，ANSYS计算的分度圆处接触应力云图如下：
图11 轮齿接触压力云图
图12 齿轮分度圆处啮合时的等效应力
10、建议
○1ANSYS Workbench中默认的接触公式：对于刚体接触是罚函数接触算法，对于其他类型的是增广拉格朗日接触算法；
○2ANSYS官方资料中说明，接触刚度越大结果越精确，但是会导致收敛困难，刚度越小收敛良好但是会影响精确度，当接触刚度太大时，如1e16，则机器的精度不能够保证，可能的话缩小模型力的单位，通常接触刚度的单位是force/length^3，基于力的模型时接触刚度单位是force/length；
○3初始刚度计算公式为：K contact=E/H/FTOLN*2，K contact为初始接触刚度，E为弹性模量，H为接触下层单元的平均深度，FTOLN为穿透容差因子。

○4大多数情况不需要调整接触刚度，但是推荐使用刚度更新，当使用接触刚度更新时，接触刚度因子的大小对最终的接触压力大小影响不大，每个迭代接触刚度的更新基于当前下层实体单元的平均应力和穿透容差而决定。

○5从接触公式来看，接触力一定的情况下，接触刚度和接触穿透量是互相影响的，使用增广拉格朗日接触算法时，系统将接触穿透作为一个收敛指标，接触刚度的大小是未知的，所以在通过调整接触刚度因子对接触结果无影响时，可以考虑调整穿透容差，也可以考虑调整单元的尺寸，虽然一般而言，网格越多结果越精确，但是无限的减小网格尺寸是不可能的。

○6使用常规拉格朗日接触算法不需要接触刚度，并且需要使用直接求解器，这种情况下的接触结果靠网格来保证，但是在网格大扭曲的情况下，增广拉格朗日接触算法比常规拉格朗日算法要好。