计算不收敛？接触问题有限元Matlab编程求解攻略(上)

一、写在文前

本次课程主要围绕接触问题的有限元Matlab编程求解进行介绍，接触问题与我们之前讲过的几何非线性问题同属非线性有限元的范畴，但相比材料非线性和几何非线性来说他更难求解，其非线性更强，相信好多用商业有限元软件做接触分析的小伙伴最常遇到的问题就是不收敛了，为啥不收敛就是因为接触分析有很强的非线性。

首先接触面的几何特征是未知的，两个物体的哪些节点会相互接触，形成的接触面什么样子，这些都是未知的，都需要反复的试错迭代才能找到；其次接触面的力学本构关系也有很强的非线性，不要简单认为接触就是两物体之间相互挤压这么简单，这只是法向上的接触力，此外还包括切向上的接触力，也就是摩擦力，摩擦力学行为本身就是有很强的非线性，我们用高中的物理学知识去想一下摩擦力，首先如果小于最大静摩擦力物体不会发生相对，如果大于最大静摩擦力物体表面会发生相对滑动，可以看出摩擦力相对位移来说是一个突变的过程，而且最大摩擦力的大小与法向的压力还相互耦合，再考虑的因素多一些还包括接触面的粗糙度，环境温度湿度等等。当然我们这节课只是接触问题有限元编程求解的入门课，不会讲这么复杂。

这次课我们主要针对图1所示的两矩形块相互接触的平面应力问题进行介绍，主要分为两个部分，一是接触理论有限元求解的基本原理和算法，二是利用Matlab对接触问题有限元求解算法进行实现，并将计算结果与商业有限元计算结果进行对比，验证了程序算法的可靠性。

一、 接触问题有限元求解原理及算法

本文以图1所示的两物体接触的平面应力问题为对象介绍接触问题的有限元Matlab编程求解过程。接触问题的有限元求解需要确定运动学关系和接触界面的本构方程、变分方程及其在接触区域内有限元的离散化。以下内容分别从这四方面进行讲解。

图1 两物体接触的平面应力问题

1、接触问题的运动学关系

本节总结了定义接触条件所需的运动关系，通过运动学关系来描述合适发生接触。接触条件可分为法向接触（法向穿透）和摩擦接触（摩擦滑移）。法向接触公式需要距离（穿透）函数来定义。摩擦接触需要相对切向相对位移来定义。本文只围绕法向接触问题基于小变形假定（有限变形）进行探讨。对于距离（穿透）函数的公式，必须区分接触的表面。在本文中，一个表面定义为从表面，另一个表面定义为主表面。通常接触分析的目的是找到接触点和对应的接触点上的接触力。在有限元分析中，位移和力知道其中一个便可通过平衡方程求得另一个。但在接触分析中，两接触表面是否接触并不直接判断，如果接触，接触点坐标和接触力均为未知，这就是接触问题的最大挑战。这就需要通过试错迭代的过程去寻找接触点，一旦形成接触对，就将接触约束条件施加上去。

所以接触分析问题第一步就是要找到接触点，即找到接触对，具体来讲就是找到与从面已知点x1对应的主面上的未知接触点x2（主从面的定义对于本文要讨论的问题是任意的），如图2所示。通俗地讲，就是找到主面上距离x1最近的点。为此，在自然坐标系（参数坐标系）下定义主面上的未知接触点。这里的自然坐标系可参考第四节内容加以理解。主面上的点x2与从面上的固定点x1的最小距离由下式确定

（1）

根据数学上的正交投影理论，上式可等价为

（2）

式中向量与向量相互垂直，为主面上的切向量，二者关系如图所示。这样就可以确定主面上对应从面节点的投影接触点x2即为距离从面固定点x1最近的点，x1与x2就形成了一个接触对。

接触点对确定之后，下一步就是通过计算x1和x2两点之间的距离来判断是否发生接触从而施加接触约束。两点距离通过下式定义

（3）

从上式可以看出，gn>0，不接触；gn=0，完美接触；gn<0，穿透。

2、接触面本构方程

根据描述接触界面处的力学行为所需的精度，文献中对接触区的本构方程建模有多种不同方法。当切向接触力学行为通过摩擦本构方程描述，法向接触力学行为通过微观力学激励或者由非穿透运动学约束来描述。由于本文只考虑法向接触，因此，上述接触约束条件即抑制穿透的约束条件为

（4）

在这个公式中，法向接触力是反作用力。将接触力考虑进去之后可以得到所谓的 Kuhn-Tucker-Karush 条件

（5）

其中pN是接触力，即在接触的情况下，由约束g = 0条件引起的反力。

3. 弱形式

引入上述公式（4）接触条件不等式之后的变分不等式如下式所示

（6）

可能大家对上述公式乍一看并不理解，对于几何线性问题，上述方程可改写为如下变分形式

（6’）

上述变分方程可以理解为虚功原理，正是由于存在接触约束，因此为内力虚功大于外力虚功。

求解变分不等式的算法种类繁多，其中大部分源于优化理论，最流行的算法就是惩罚因子法和拉格朗日乘数法。本文只讨论惩罚因子法。

当将有限元方法与惩罚或拉格朗日乘子方法一起使用时，不等式约束被公式化为激活和非激活的接触约束，这些约束在求解过程中是动态变化的。因此，上述变分不等式（公式6）)可改写为包含激活约束的变分方程

（7）

对于两个物体之间的接触，可通过惩罚因子方法来确定被激活的接触区域。公式（7）中的接触贡献项（contact contributions）可通过下述弱形式的惩罚项引入法向的接触约束gN = 0。

（8）

式中即为惩罚因子，惩罚因子是一个很大的数，理论上讲，时趋于真实约束，但是惩罚因子取值过大会导致数值求解问题。将罚系数取为一个较大的数，一般可以取刚度矩阵中最大元素的1e6倍。

4、有限元离散

在对接触区域做有限元离散时，单元节点在接触界面处不匹配，如图所示。在接触算法实现时，必须区分接触表面上的点哪些是接触的点，哪些是不接触点。为此，要定义所有可能发生接触的节点的集 合，这一过程对应我们在商业有限元软件中定义接触面的操作，接触面包括主面和从面，是可能发生接触的一对面。在实际问题中，并非所有这些接触面上的点都发生接触力学行为，只是定义一个可能发生接触的节点集 合，在计算中只在此集 合中进行搜索判断是否发生接触。因此为了节约计算资源，尽可能减少集 合中的元素，但又不能太少导致发生不合理的穿透行为。

如下图3所示，可能发生接触的节点集 合定义为Jc，在计算中的某一时刻其被激活的节点定义为JA。。在后续矩阵公式中，仅包括被激活的节点。

图3 接触问题的有限元离散

（1）Node-to-segment（NTS）离散

对二维问题，使用Node-to-segment（NTS）方法离散，这里node代表节点，segment可以理解为描述物体边界的两个相邻节点的连线。如图所示，假设从节点x1s，与x21，x22定义的主线段接触。

图4 NTS离散方法（图中1，2上标有误，应互换）

通过参数坐标 ξ 和主面节点的线性插值来定义接触点，

（9）

对应的切向量为

（10）

对切向量进行归一化，得到 a21 = ̄ a21 /l，其中主单元的长度 l =‖ x22 − x2 ‖。根据单位切线向量a21，可推导出主线段的单位法向量为n2 = e3 × a21，e3表示面外单位向量。

根据公式（2）和公式（3），可求得与x1s对应的主面上投影接触点的参数坐标ξ 和从面节点x1s到主面的距离gNs，分别如下式所示。

（11）

上述距离函数的变分可表示为

（12）

η 表示试函数或虚拟位移。上述变分表达式在可以在弱形式 (7-8) 中使用，其中积分应当在实际被激活的接触面范围内进行计算。

（13）

式中表示所有被激活的接触面的个数。上式中的节点法向力 PNs 由 PNs = pNs*as 在从节点处给出。对于惩罚方法，给出，其中等于拉格朗日乘子法中的拉格朗日乘子，代表接触力的大小。在二维情况下，接触单元的面积as等于 （对应主节点间的距离），但是如果一个主段（segment）上有多个从节点时，如图5所示，asi可通过下式计算得到

图5 主段上有多个从节点与之接触

对于公式(12)表示的距离函数变分，写成矩阵的形式为

（14）

（15）

（16）

这样一个从面节点的接触贡献的弱形式为 η*Gs ，其中对应的单元残差Gs表示为：

（17）

（2）算法

接触算法的核心就是搜所可能发生接触的物体确定实际被激活的接触面。一旦激活接触，就必须建立基于穿透函数 (5) 的局部接触约束条件。处于被激活接触状态的节点总数用nc表示。弱形式的全局的矩阵表达式为

（18）

其中Gv表示公式（7）中两物体贡献的弱形式，Gsc表示接触贡献的弱形式，具体由公式（17）定义。

具体算法如下：

(1)初始化v1=0;

(2)对所有节点i进行迭代循环求解直至收敛

a.搜索是否发生接触：当gNsi ≤ 0 → 激活节点接触条件

b.计算新的位移增量

c. 收敛性检查：满足则结束迭代，否则继续迭代。

就是我们熟知的切线刚度矩阵。该矩阵包含了两物体各自的刚度以及接触贡献的刚度，在加载过程中，每时刻的接触状态（接触面，接触力）都在不断变化，因此存在非线性，需通过牛顿迭代法来确定当前时刻的切线刚度矩阵。

（3）接触残余力的线性化

求解上述接触非线性的算法是牛顿迭代法，这就需要计算切线刚度矩阵。法向接触问题的正切矩阵来自公式 (13) 中的项 δgNs PNs。在第一项δgNs (14-16)的线性化中，还必须考虑 ̄ξ 对当前位移的依赖性以及法向量的n1的变化。根据可以得出切线刚度矩阵的表达式为（参考P404 5.65,对比11.70），

考虑到gNs是一个极小值，因此，与gNs相乘的项均和省略

（19）

其中，Ns的定义见公式16。

需要注意的是，本文只考虑法向接触，且为两个平面接触，因此在接触之后的整个过程中，接触刚度不会发生变化，因此无需通过牛顿迭代法求解切向刚度矩阵。大家重点关注以下几点：(1) 接触搜索算法，(2) 判断是否施加接触约束，(3) 接触刚度矩阵的计算，其中前面两点最终服务于第三点。因此接触刚度矩阵计算是求解接触问题的核心，具体算法如下:

1）搜索接触：在预定义的主面和从面之间遍历所有可能的node-to-segement组合；

2）判断是否施加接触约束，条件一：node与segement的距离gN<0 ；条件二: node在segement上投影点的参数坐标满足 ；

3）根据公式19计算刚度矩阵。

三、Matlab编程实现接触问题的有限元求解

由于篇幅过长，笔者将在下一篇文章给大家介绍。目前我已经录制了对应的视频教程发布在《Matlab有限元编程从入门到精通20讲》后面进行加餐，欢迎识别下方二维码试看。

我的Matlab有限元编程精品课

Matlab有限元编程从入门到精通20讲：快速获得各典型有限元案例的Matlab代码