Ansys Workbench工程应用之——结构非线性（下）：状态非线性（3）间隙接触
存在初始间隙而在加载后间隙闭合的情况都可以统称为间隙接触问题，在两个几何体分离时，接触是没有刚度的，在接触瞬间，刚度突变，所以往往会出现不收敛问题，本文主要讨论的便是接触稳定性问题。
1 概述
线性接触只在初始检测接触状态，而非线性接触在每个子步更新接触状态。所以存在初始间隙的模型计算接触碰撞过程中，程序将在每一子步检测接触的两个对象是否已经闭合。如果第一子步就发生刚体无阻尼的无限距离平移，便会发生接触失败和刚体平移穿透现象，这也是载荷加载接触失败的原因，这将在下文详解。
此处需要注意，在间隙接触分析中应谨慎使用几何修正，因为它们都是采用数字补偿的方法调整间隙，可能与我们所需的计算要求不同，反而使计算结果不准确。比如下图轴孔配合，左图使用调整接触，相当于在轴与孔之间的间隙处创建了刚性区域，右图未使用调整间隙，在轴受横向力的工况下，明显右图更符号实际情况。
计算前推荐（不是必须）检查模型间隙是否在接触球半径范围内，因为对于非线性接触，只有间隙在接触球半径内，程序才认定两面之间为近场(接近打开)状态，而不是远分离状态，计算效率更高。
2 加载方式的影响
间隙接触问题计算方法有控制位移加载法、控制载荷加载法、显式动力学方法等，其中控制载荷加载法不易收敛，必须设置接触稳定因数。显式动力学不在本文讨论范围内。
控制位移加载法比控制载荷加载法更容易收敛，因为在分子步后，每一子步加载位移是确定的，不会发生刚体不受控的无限平移。接触的两个对象总会在某个子步或某两个子步之间发生闭合，使接触建立，所以控制位移加载往往无需设置接触的稳定阻尼也能收敛。
位移控制的缺点是只有在知道施加什么位移时才适用。
实例1 位移加载下的球面碰撞问题。
Step1 建模。
建立静态结构算例，创建100×100×5mm的平板，在平板上方创建R15的球体，两者间隙为1mm。
Step2 接触设置。
设置球面与平面的无摩擦接触，球面为接触面，平面为目标面，接触行为不对称，其余默认。
设置球体网格方法为六面体主导，尺寸3mm。平板在厚度方向设置2层单元，其余默认。
Step3 边界条件。
打开大变形选项。
固定平板的4各端面，给球体施加向下10mm的位移约束，水平2个方向位移为0。
Step4 求解与后处理。
求解后总变形如下。
平面变形如下，最大变形为8.9mm。
在结果中插入接触工具——渗透，如下图，最大渗透为0.14mm，对比平面变形比为0.14/8.9=1.5%，比值较小，完全可接受。
在结果中插入探针——力反应，分别设置固定约束处的反力、定向位移的反力，可以发现它们是平衡的，说明全局反力满足平衡。
注意：本例还可使用法向拉格朗日算法，计算结果渗透接近0。
2.2 控制载荷加载
控制载荷加载法不容易收敛性，特别是间隙较大的情况，所以控制载荷加载法一般只用于小间隙闭合分析，且需要使用稳定性控制。
若不使用稳定阻尼，哪怕第一子步所分配的力再小，由于没有抵抗力，刚体将在发生不可控的无限平移，发生刚体平移穿透现象。
稳定性阻尼系数(Stabilization DampingFactor)在接触设置的高级选项中。稳定性阻尼应用于法向，只对非线性接触有效，接触稳定性阻尼添加粘性阻尼牵引力，该力与虚拟的速度成正比，但是与刚体运动反向。
接触稳定性控制类似于非线性控制选项中的稳定性，向模型中引入了伪能量。这个方法可以增强收敛，但是当引入的阻尼力太大时，会影响求解精度。所以一般不要超过0.2，建议设置为0.01~0.1之间。
实例2 载荷加载下的球面碰撞问题。
Step1 分析前设置。
复 制实例1算例，接触算法改为程序控制，其余均不变。
固定平板四周，给球体施加竖直向下的1000N的力，约束球体只能在竖直方向移动，水平的两个方向位移为0。
Step2 初次求解。
求解后报错不收敛。
Step3 稳定性控制。
在接触中将稳定因素设置为0.1，此系数一般需要反复调试，取能收敛的较小值，使稳定伪能量对求解结果影响最小。
设置分步如下。由于碰撞接触是强烈非线性问题，建议将牛顿法改为非对称，有利于收敛。
Step4 求解与后处理。
求解后平面变形如下图，最大值为0.028mm。
接触渗透如下图，最大值为0.0015mm，相对于平面变形为5%，属于较小水平。
支反力验证如下图，固定约束处的支反力为1000N，系统处于平衡状态。
3 间隙接触收敛控制
除了上文提到的稳定性控制以外，还可以通过以下方法控制间隙接触问题的收敛。
3.1.1 接触行为
非对称接触也被称为“一次性接触”，这通常是模拟表面到表面接触的最有效方法。而对称接触可以在接触表面之间生成两组接触对，这称为对称接触。与非对称接触算法相比，对称接触算法在更多表面位置强制实施接触约束条件。
显然，对称接触的效率低于非对称接触。然而，许多分析将需要使用它（通常是为了减少渗透）。需要对称接触的特定情况包括：
（1）接触面和目标面之间的材料相近（刚度相近），形状相近（其实不相近也行）。
（2）两个表面都有非常粗糙的网格。
所以当出现因网格粗大而导致的穿透与不收敛问题，读者可以尝试使用对称接触行为，以实例2为例，采用默认的程序控制算法 接触稳定 默认网格情况下，使用非对称接触经过50次迭代收敛，使用对称接触经过29次迭代收敛。
3.1.2 法向拉格朗日算法
在接触求解中，很多时候的不收敛是由于穿透过大引起，此时可在高级——公式化中设置拉格朗日算法，可得到接近于0的穿透量，有利于收敛。
法向拉格朗日只支持接触非对称行为，当接触行为设置为程序控制、自动不对称、对称时，内部计算器都会启用自动不对称算法。
法向拉格朗日只支持直接求解器（即稀疏矩阵求解器），不支持迭代求解器。程序控制将自动调用直接求解器。当用户将求解器强制设置为迭代求解器后，程序仍然调用直接求解器并在求解信息中发出警告。
法向拉格朗日不支持基于积分点的探测方法，所以探测方法中是没有积分点选项的。
3.1.3 接触时步控制
接触时步控制只支持非线性接触，默认不启用。它的作用是当接触发生突变时，允许时间步自动二分，以便减小载荷（或位移），促进收敛，副作用是会消耗更多时间和计算资源。
自动平分即自动二分，默认情况下，接触状态改变是不会影响时间步长，选择此项后，在每个子步骤结束时审查接触行为，以确定是否发生过度渗透或接触状态的剧烈变化。如果是这样，则使用二分的时间增量重新计算子步。
预测影响即预测碰撞冲击，在上述选项基础上，增加预测检测接触行为变化所需的最小时间增量。如果接触中有冲击，或者遇到间隙接触不收敛，建议使用此选项。
第三个选项是使用影响约束，用于没有能量损失的弹性体碰撞问题。
3.2 非线性控制中的相关设置
3.2.1 牛顿选项
如果遇到接触摩擦问题不收敛，特别是摩擦系数＞0.2时容易发生不收敛，可以将牛顿选项改为非对称，比如以下齿轮接触啮合问题，
3.2.2 线搜索
因为线性搜索的主要作用是程序在检测到硬化效应（例如接触突变时），就在位移增量上自动乘以0~1之间的系数，相当于给位移减速，特别有助于间隙接触问题的收敛。用户只需要保持程序控制（或手动强制打开）即可，不要手动关闭分析设置——非线性控制——线搜索功能。
例如以下2D接触问题，自动时步均设置如下，线搜索使用程序控制（或手动强制打开）与关闭线性搜索的对比，前者明显更易收敛。
3.2.3 非线性稳定性
在分析设置——非线性控制中，有稳定性选项，这已经在屈曲专题中详细介绍过了。它也可以用于间隙接触的收敛控制，而且与接触中的稳定性原理相同，但是需要注意，必须激活第一个子步，因为在载荷加载的间隙接触问题中，产生刚体平移的正是第一个子步。需要注意的是，在接触中设置了稳定性，就不要在非线性中设置稳定性，即只能选择在其中一处设置稳定性，否则可能引入过大的稳定伪能量。
如果遇到不收敛问题，读者可在求解方案信息——Newton-Raphsion残余力总设置为1，计算后生成不收敛的最后一子步的残余力，读者可根据云图显示分析是否网格过于稀疏（应力残余发生在网格中间区域），是否穿透过大等。
附 文献【1】中相关实例
实例3 载荷加载下的接触非线性。
燕尾槽广泛用于机床导轨设计中，与梯形导轨配合使用,由于槽内需要储存润滑脂，所以槽与导轨的配合为极小的间隙配合。下面以一个简单的燕尾槽及梯形导轨为例，说明间隙配合接触分析过程。
Step1 建模。
建立静态结构分析项目，点击A3（几何结构），将分析类型设置为2D。
以XY平面建立二维矩形，单位均为mm，矩形长为36，高为20，在矩形内掏梯形燕尾槽，槽上宽12，下宽30，梯形低角60°，圆角为2.5mm。以槽尺寸向内偏移0.25mm建立一梯形导轨。
在Mechanical中在Geometry的细节菜单中将2D Behavior定义为平面应变，这是因为导轨及燕尾槽都很长，只取一截面分析,所以为平面应变。其余默认设置。
Step2 网格与接触。
网格仅定义全局尺寸2mm，质量——误差限值设置为强力机械，其余默认设置。
设置接触：轨道的5条边为接触侧，燕尾槽的5条边为目标侧，接触类型为无摩擦，行为为不对称，接触方法为法向拉格朗日（可得到近似于0的穿透量）。稳定阻尼因子为0.04（反复测试得到的较小值）。
插入接触工具并计算初始接触，确认接触状态不是远打开。
Step3 边界条件。
分析设置中，设置1步，初始子步10、最小子步10、最大子步500，打开大变形。
边界条件：选择燕尾槽的下边加载固定约束，选择导轨的上边加载力,大小为10N，方向沿Y轴方向。
Step4 计算与后处理。
计算顺利收敛。
查看变形云图，中变形如下，最大位移为0.5mm，这其中主要是轨道的刚体平移。
燕尾槽变形如下图，最大变形仅为0.002mm，位于燕尾槽上端，说明计算结果是合理的。
验证穿透量：插入接触工具——渗透，结果显示穿透为0。
验证支反力：插入固定约束的支反力，求得支反力为10N，说明系统是平衡的，结果可信。
实例4 位移加载下的材料非线性 接触非线性。
接触分析过程中，往往伴随着材料非线性特征，这两种非线性结合在一起，极易不收敛。下面以一个2D拉拔型材的例子说明这类接触分析方法。
Step1 建模。
建立静态结构分析项目，点击A3（几何结构），将分析类型设置为2D。
建立一拉拔模及型材，单位均为mm。拉拔模为5边形，总高30，左侧边长高15，距中心线8，上边宽11，下边宽6；型材左侧与Y轴对齐，总高28，宽9，与模具重合高度为8，圆角为2.5。建立2D模型的目的是为了计算快捷，让读者更方便地了解计算收敛的方法。
Step2 材料设置。
双击A2（工程数据）进入工程材料库。空白处右击——工程数据源。调用WB默认材料库内的一般非线性材料——铝合金NL。如果没有，可以关注后私信我。
在Mechanical中，在Geometry的细节菜单中将2D Behavior定义为轴对称。
修改模型名称分别为模具和型材，定义模具的材料为结构钢，定义型材的材料为铝合金NL。
Step3 网格与接触。
按默认设置生成网格。
设置接触：选择型材的上边线及圆角线两条线为接触边，模具的斜线和左边线两条线为目标边，接触类型为摩擦、 摩擦系数0.1，其余设置均默认。
Step4 边界条件。
对模具整个面加载完全约束，选中型材的下边线定义位移约束，X向为0， Y向沿Y轴正方向30mm。
分析设置：设置分步，初始子步50，最小子步10，最大子步500。打开大变形开关。输出节点力设置为是。
Step5 首次求解。
在书上此处不能收敛，但是随着ANSYS版本不断更新，现在的版本在非线性方面已经强化了不少，所以能收敛。插入总变形结构如下。
注意看云图，箭头所指位置型材已经发生了畸变，这是由于网格太稀疏造成。椭圆标识位置发生了穿透，是因为型材的右边线没有与型材边线建立接触。
Step6再次求解。
网格细化：对圆角进行网格细化，对圆角边线划分8分，再次划分网格。
接触修改：将接触边改为型材的3上边、圆角、右边共3条边。
再次求解，总变形如下图。
渗透验证：在结果总插入接触工具——渗透，最大渗透仅为0.008mm（如果使用法向拉格朗日算法，可得到近似于0的渗透结果），属于很小水平，完全可以接受。
摩擦应力：在结果总插入接触工具——摩擦应力，整个过程的摩擦应力最大约为100MPa。
接触反力：在结果总插入摩擦接触的反力，注意此操作前必须在分析设置——输出设置中打开节点力输出。
位移约束反力：在结果总插入位移约束的反力，在0.95左右达到最大值，Y向此时反力为76000N左右，可以作为拉拔型材所需的力。
扩展显示如下图设置。
写在最后
对于间隙接触问题，由于接触初始有间隙，从而产生刚体平移或接触刚度突变导致收敛困难。使用位移加载更容易收敛，在使用载荷加载时必须使用一些的收敛奇淫技巧才能收敛。
本期解读了间隙接触，下期将详细解读过盈配合，敬请期待。
由于图惜实践经验实在有限，文中也难免纰漏百出，敬请批评指正。
[1]《Ansys Workbench有限元分析实例详解》——周炬、苏金英
[2]《ANSYS Mechanical连接功能系列课程》——安世亚太苏睿
[3] ANSYS 2022帮助文件
