数值算法：解决接触问题的优化方案

接触是一类极度不连续的约束，因此有限元法求解接触问题的过程是搜寻接触状态的反复迭代过程。为此，先假设一个可能的接触状态，然后代入约束边界条件，得到接触点的接触内力和位移，判断是否满足接触条件。当不满足接触条件时修改接触点的接触状态重新求解，直到所有接触点都满足接触条件为止。其过程如图2.15所示。

图2.15　接触问题的有限元法计算流程

有限元法求解接触问题的关键为接触算法，主要可以分为直接迭代法、接触约束算法和数学规划法等。其中，接触约束算法是通过对边界约束条件的适当处理，将约束优化问题转化为无约束优化问题求解，主要可以分为罚函数法、拉格朗日乘子法和增广拉格朗日乘子法。

1.罚函数法

罚函数法（Penalty Method）是用罚函数引入接触约束定解条件的附加泛函

Πp=αg2　（2.54）

式中，α称为罚参数或罚数；g为间隙（包括切向间隙）。

与原问题不包含接触约束条件的总位能Π（U）一起构成修正的泛函。将约束条件极值问题转化为无条件极值化问题

minΠ*(U)=Π(U)+Πp(U)　（2.55）

罚函数法的最大优点在于引入接触条件时无须增加系统的自由度，不增加计算机原存储量和计算量，而且使求解方程的系数矩阵保持正定，因此得到较广泛的应用；缺点是约束条件只能近似地被满足。理论上，增大罚参数可以使计算精度提高，但是罚参数过大容易导致方程病态。

2.拉格朗日乘子法

对于包含接触界面的接触问题，泛函可表示为

Π=Πu+Πc　（2.56）

式中，Πu是原问题中不包括接触约束条件的总位能；Πc是用拉格朗日乘子法引入接触约束条件的附加泛函。

Πc=λg　（2.57）

式中，g的意义同前。(www.xing528.com)

与原问题不包含接触约束条件的总位能Π（U）构成修正的泛函，将约束条件极值问题转化为无条件极值问题

minΠ*(U)=Π+Πc　（2.58）

拉格朗日乘子法可以精确地满足接触约束条件。但是由于拉格朗日乘子的引入，即对每个接触对引入了一个新的未知量，增加了方程的自由度数，扩大了系统的求解规模。而且求解方程的系数矩阵不再保持正定，对每一个乘子项，都在对角线上为零值，必须采取适当的方法以保证方程求解的收敛性和稳定性。

为了克服拉格朗日乘子法中系数矩阵存在零主元的弱点，在修正泛函式（2.58）中增加一惩罚项

式中，Ep为罚因子，Ep＞0。将式（2.58）改写为

minΠ*(U,λ)=Π+Πc+Πε(Ep)　（2.61）

与式（2.58）相比，由于多了λ的二次项，经变分后得到的有限元计算格式中对应于λ的主对角元素不再为零，总体刚度矩阵不再为奇异矩阵。当罚因子Ep趋于无穷时，式（2.61）的解收敛于式（2.58）的解。这就是摄动拉格朗日乘子法（Perturbed Lagrangian Method）的思想。

3.增广拉格朗日乘子法

由于拉格朗日乘子法和罚函数法各有优缺点，将它们联合起来，便形成了各种增广拉格朗日乘子法（Augmented Lagrangian Method）。如在式（2.61）的基础上通过引入附加泛函

Πp=αg2/2　（2.62）

构造修正的势能泛函

Π*=Π+Πc+Πp　（2.63）

通过应用

λk+1=λk+αg　（2.64）

计算拉格朗日乘子的更新，而不增加系统自由度。由于引入了高次罚因子项Πp，增加了方程矩阵的对角优势，可以改善求解收敛性。因此，增广拉格朗日乘子法成为目前较为流行的方法。但与罚函数方法一样，其求解精度同样依赖于罚参数值的大小。