非线性有限元编程 | 接触（1）
今天开始，将要更新接触方面的非线性有限元编程知识，本篇推文会给大家带来以下内容：
接触相关概念 结合单点接触问题，使用理论解析方法求出接触力 分别使用拉格朗日乘子法和罚函数法求解接触力
接触相关概念
接触是将多个部件连接在一起的主要工具，包括螺钉、螺栓、焊缝等，接触非线性可以被归类为边界非线性。
接触分析的目的是回答以下问题:
是否有两个或多个物体在接触； 如果有接触，接触的位置或区域在哪里 接触界面中有多少接触力或压力 接触后界面中是否存在相对运动
接触的两个特点：
若两个单独的物体接触，当两个物体分开时，接触力保持在零；物体接触后，接触力瞬间增加 在位移边界上，若位移已知，则可以计算出反力；在牵引边界上，若已知作用力，则可以计算出相应的位移。然而，在接触的情况下，_位移和接触力都是未知的_，除了非常有限的情况。
在求解接触问题时，可以将接触视为优化公式的约束，使势能在满足接触约束的情况下最小化，即一个物体不能穿透另一个物体。通常使用罚函数法（penalty regularization）和拉格朗日乘子法（Lagrange multiplier methods），将约束优化问题转化为无约束优化问题。
在有限元分析中，主从面（the master/slave surface）的概念常常用来描述接触问题，通过施加特定的接触约束，使得从属边界上的节点不能穿透主边界上的表面单元。
Kim 教授的书中主要从以下几个方面讨论接触非线性：
单点接触 基于变分法提出接触通用公式（柔性接触-刚性接触、柔性接触-刚性接触、线-线接触、面-面接触、摩擦等） 接触变分形式的有限元离散和数值积分 三维接触公式
单点接触问题
考虑受均布荷载作用的悬臂梁，梁的自由端受刚性块的限制，如下图所示，梁的末端与刚性块之间存在一个小间隙。均布荷载 =1 kN/m，梁长 =1 m，抗弯刚度 EI= N ，初始间隙 =1 mm。
直接法
我们先忽略刚块的存在，悬臂梁（Euler beam model）在均布荷载作用下，自由端的位移为：
得出的挠度值大于间隙 ，即可判定两者发生了接触，假定接触只发生在梁的尖端，即单点接触。由于刚块阻止了梁的弯曲，其效果可以通过施加一个力来模拟，即一个接触力（ ），使梁不能穿透刚性块。
由于梁的挠度较小，因此对两种荷载的影响采用叠加规则，梁在末端力（负号表示接触力与均布荷载方向相反）作用下的挠度曲线和末端挠度可表示为：
根据叠加原则，
计算得出接触力 为 75N。
拉格朗日乘子法
直接发为了确定接触状态，必须首先计算出解(梁的挠度)。当接触问题变得复杂（如多点接触问题）时，往往不能求得解析解。这时可以将梁与刚体之间的间隙作为未知量，加上一个附加约束，可以得到一个更系统的接触公式。
未知的接触力用 表示，它作用于梁和刚体的相反方向。为了分配一致的方向，两个接触点中的一个被认为是主接触点，另一个被认为是从接触点。
将接触力视为外加荷载，利用两个独立荷载的叠加规律，可得到梁的挠度曲线：
接触的物理要求是无穿透，接触力为正，当间隙大于零时，接触力为零，反之亦然。这些要求要求解决方案满足以下三个条件：
可以得出方程：
求得接触力 =0 或者 75N，当 =0 时， =0.00025，物理意义是，不存在接触力时，间隙为 0.00025，也就是悬臂梁在承受均布荷载时的自由端竖向位移；当 =75N 时， =0，物理意义是当接触力为 75N 时，梁端刚好与刚体接触。所得结果与直接法一致。
在以上求解过程中，附加的未知(接触力)称为拉格朗日乘子，并利用一致性条件来确定接触状态和接触力，该方法被称为拉格朗日乘子法。
罚函数法
在处理接触问题中，还有一种常见的方法，就是罚函数法。什么是罚函数法呢？
刚刚讲的拉格朗日乘子法中，不允许从面侵入主面，那换个角度想，如果侵入了，该怎么处理呢？
此时，来了一个罚函数法，允许接触过程中有少量侵入，且接触力与侵入量成正比，即 可为正也可为负，定义渗透函数 ：
接触力可由渗透函数定义：
其中， 被称为罚系数，有时也被称为罚刚度。可以想象一下，接触过程中有了少量的渗透，然后通过施加一个大的力来“惩罚”它。将接触力带入 ，得到：
通过赋予 不同的值，可得到不同的渗透量，及接触力，汇总如下表：
从上表可以看出，随着罚系数量级的增大，渗透量逐渐减小，接触力逐渐收敛于解析解 75N。
木木最近也在找有限元开发性质的工作，面试的时候发现企业类的大多喜欢有接触经验、弹塑性经验或者动力分析经验的应聘者，而且所使用的语言也多为C++，正好最近也在学习非线性相关的理论，那就顺道一起更新一下接触方面的知识~
如果你有合适的有限元开发类的工作推荐，欢迎联系木木！！！
【声明】：本次案例分享来自Kim教授的《Introduction to Nonlinear Finite Element Analysis》。
