接触问题有限元Matlab编程求解(上)|接触问题（点-线接触） | Matlab源码 | 理论文本 |接触非线性

suoge223

已于 2024-03-18 10:22:40 修改

阅读量2.9k

点赞数 20

分类专栏：有限元编程从入门到精通文章标签：人工智能

于 2024-02-22 10:30:13 首次发布

本文链接：https://blog.csdn.net/u010542847/article/details/136223916

版权

有限元编程从入门到精通专栏收录该内容

36 篇文章

订阅专栏

专栏导读

作者简介：工学博士，高级工程师，专注于工业软件算法研究
本文已收录于专栏：《有限元编程从入门到精通》本专栏旨在提供 1.以案例的形式讲解各类有限元问题的程序实现，并提供所有案例完整源码；2.单元类型包含：杆单元，梁单元，平面三角形单元，薄板单元，厚板单元，壳单元，四/六面体实体单元，金字塔单元等；3.物理场问题涉及：力学、传热学、电磁学及多物理场耦合等问题的稳态（静力学）和瞬态（动力学）求解。专栏旨在帮助有志于有限元工业软件开发的小伙伴，快速上手有限元编程，在案例中成长，摆脱按部就班填鸭式教学。
【所有专栏文章均提供对应视频课+源码】，文本教程+视频教程+源码联系，三向强化学习！学习地址： https://www.bilibili.com/video/BV1kP4y1d7Zo/
欢迎订阅专栏，订阅用户可私聊进入有限元编程交流群（知识交流、问题解答），并获赠丰厚的有限元相关学习资料（教材、源码、视频课）
专栏订阅地址：有限元编程从入门到精通_suoge223的博客-CSDN博客

文章目录

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

前言

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

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

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

【课程完整案例源码链接】：接触有限元编程 | Matlab源码 | 理论文本

本次课程主要围绕接触问题的有限元Matlab编程求解进行介绍，主要分为两个部分，一是接触理论有限元求解的基本原理和算法，二是利用Matlab对接触问题有限元求解算法进行实现，并将计算结果与商业有限元计算结果进行对比，验证了程序算法的可靠性。今天我们重点给大家介绍接触问题有限元Matlab编程求解第一部分。

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

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

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

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

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

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

（1）

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

（2）

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

图2 接触对的定义（图中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对接触问题有限元求解算法进行实现。本文以图1所示的两物体接触的平面应力问题为对象介绍接触问题的有限元Matlab编程求解过程。图1中上部块体上表面受均布荷载作用，左图为Matlab计算结果，有图为abaqus模拟结果，经比较，二者应力分布形式基本一致，但因采用算法的不同计算结果略有差异。本文的接触问题的有限元算法在计算精度方面还有待提高，但作为接触问题有限元编程入门来讲还是有一定借鉴意义，姑且抛砖引玉吧，欢迎大家批评指正！

图1 两物体接触的平面应力问题Matlab有限元编程求解

图2 两物体接触的平面应力问题Abaqus求解

本程序基于四节点平面应力单元的代码进行开发，四节点平面应力单元相关的理论介绍和代码实现可参考文章《四节点八节点四边形单元悬臂梁的Matlab有限元编程》或仿真秀视频教程《Matlab有限元编程从入门到精通20讲》中对应的四节点四边形单元的课程。

1、主程序main.m

Main.m是主程序，部分代码展示如下。由于接触问题是两个独立的物体相互接触的过程，所以主程序对于两物体的节点单元信息、边界条件以及各自的弹性刚度矩阵都是分开定义。其中Geometry()返回节点单元信息(P,E)、荷载施加节点(q0_node)和边界约束施加节点(Boundary)，通过Patch=1或2分别定义上下物体的上述信息；然后主程序通过while循环，分100个荷载步逐步施加荷载，并计算每个荷载步对应的节点位移d，然后通过d_final=d_final+d与上一步计算的节点位移d叠加。节点位移的计算也是通过Kd=F求得，而这里的K即本篇part1中提到的切线刚度矩阵，包含了两物体各自的弹性刚度和接触刚度两部分。两物体各自的弹性刚度矩阵通过Linear_Elasticity()函数得到，接触刚度矩阵通过Contact_Formualtion()。Linear_Elasticity()函数可参考文章《四节点八节点四边形单元悬臂梁的Matlab有限元编程》或仿真秀视频教程《Matlab有限元编程从入门到精通20讲》中对应的四节点四边形单元的课程。本篇重点讲一下，接触刚度矩阵的编程求解方法。

2、接触刚度矩阵求解程序Contact_Formualtion()

接触刚度矩阵计算是求解接触问题的核心，具体算法如下:

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

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

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

在实际问题中，并非所有这些接触面上的点都发生接触力学行为，只是定义一个可能发生接触的节点集合，在计算中只在此集合中进行搜索判断是否发生接触。接触算法的核心就是搜所可能发生接触的物体确定实际被激活的接触面。一旦激活接触，就必须建立基于穿透函数 (5) 的局部接触约束条件，进而产生接触刚度。接触刚度矩阵求解函数Contact_Formualtion()首先通过接触面几何定义函数Contact_Geometry()返回主面和从面的节点信息和对应的自由度的全局索引值，索引信息方便被激活的接触对的刚度矩阵组装至整体接触刚度矩阵中。然后程序对所有从节点和主面线段（segment）的组合进行搜索，通过cntelm2d()函数来判断是否形成接触对，同时判断接触对是否被激活从而贡献刚度，并返回对应的接触对刚度矩阵和接触力。接着将接触对刚度矩阵组装至全局接触刚度矩阵中。因此，cntelm2d()函数对于整体刚度矩阵的计算至关重要，接下来一小节就介绍接触对刚度矩阵的计算程序。

3、接触对刚度矩阵求解程序cntelm2d()

接触对刚度矩阵可根据part1中的公式19进行计算，具体的实现程序如下所示。程序首先求得segement线段的单位切向量和法向量，根据公式3计算出从节点x1到主面segement的距离gNs，程序中用GAPN表示；接下来就根据gNs的正负来判断是否发生穿透，如果为负值，说明发生node to segement接触，根据公式11计算从节点x1对应的主面上投影接触点的参数坐标ξ，程序中用alpha表示。然后通过判断alpha是否在[0,1]之间来判断从节点是否落在主线段segement上，如果落在主线段上，则根据公式1计算Ns，程序中用NN表示，然后根据公式19计算接触对刚度矩阵。这里需要补充介绍接触力的求解方法。对于从节点在接触点ξ处的接触力pN根据公式13给出，该接触力分配到主线段segement上的两节点力为

（20）

将接触对的接触力卸载一起可以得到整体接触对的接触力向量为

(21)

参考文献：

目前我已经录制了对应的视频教程发布在《Matlab有限元编程从入门到精通20讲》

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

本课程为matlab有限元编程专题课，课程主要以案例的形式进行讲解，中间会穿插案例中所涉及到的有限元基本理论，案例不局限于力学问题的有限元求解，还会涉及传热学、电学等问题的有限元求解。所涉及的单元有杆单元，梁单元，平面三角形单元，薄板单元，厚板单元，四面体实体单元等等，力学问题有静力学问题，也有动力学问题，还会涉及材料非线性、几何非线性、接触非线性等三大类非线性问题。

希望大家关注笔者仿真秀专栏，点击文章末尾阅读原文即可关注，大家学习Matlab编程过程中，有任何疑问都可以联系仿真秀客服进入仿真秀Matlab有限元编程用户交流群，群里还会分发更多Matlab编程学习资料。