0. 本文的用意和写法

该系列博文面对对有限元解析有所了解，但未接触过接触解析的初学者。本文以实用性为第一目标，试图在概念准确的前提下追求直观描述，特别强调给出某些实用算法的细节，并附有相应程序。由于这些算法细节都经过本人实装测试，其有效性是有一定的保证的。希望读完此系列博文可以自己动手组装自己的解析程序，并希望能够指出其中问题或错误，共同提高。

本文的解析对象是泛函的极值问题，这是本文的起点。对泛函的极值问题不了解的读者可以跳过泛函的变分运算步骤，着重于其后的物理概念和算法的讨论。

1. 考虑约束条件的最优化算法

泛函的极值问题的求解从数值计算方法的角度来说是最优化算法，而接触问题为其约束条件。罚函数法，Lagrangian乘数法，递增Lagrangian乘数法*1，内点法(Interior point method, barrier method或内部罚函数法)都是标准的考虑约束条件的最优化算法。这一节将以图1的问题为例解释这几种解法。

图1所示问题如没有接触约束，其泛函J为

(0.1)

这里我们将m的位移记为x，外力F<0。该泛函的极值条件为其变分为零，即

(0.2)

基于变分δx的任意性，要求

(0.3)

即得

(0.4)

在此我们假设足够大使得，但是接触条件要求物体m不能进入接触面, 即。如 , 可以用下述方法导入接触约束。

1.1 罚函数法

采用罚函数法考虑接触约束，其泛函J为

(1.1)

这里为罚函数。该泛函的极值条件为其变分为零，即

(1.2)

基于变分δx的任意性，要求

(1.3)

即得

(1.4)

当罚函数足够大时，从该公式可以得到逼近于精确解的近似解。

1.2 Lagrangian乘数法

采用Lagrangian乘数法考虑接触约束，其泛函J为
(2.1)
这里λ为Lagrangian乘数,这是一个待定自变量。泛函J的极值条件为其变分为零，即
(2.2)
基于变分δx和δλ的任意性，要求
(2.3) &
即得
(2.4)
这是该问题的精确解。x为m的位置，λ则为约束反力。

1.3 递增Lagrangian乘数法

采用递增Lagrangian乘数法考虑接触约束，其泛函J为*3
(3.1)
注意到该泛函既包括Lagrangian乘数λ也包括罚函数  。但与前述Lagrangian乘数法和罚函数法不同的是，在该算法中
* Lagrangian乘数λ不是自变量
* 罚函数  可以不是很大的数值

-----------------------------------------------------------------
注：递增Lagrangian乘数法可以看成是罚函数法或Lagrangian乘数法的变种。如从Lagrangian乘数法出发，我们将（2.1）改写为

这里，我们记  为Lagrangian乘数的准确值。考虑  为Lagrangian乘数的增分并将此式与（3.1）比较，可以看出

如此Lagrangian乘数的准确值可以如下逼近

对应于下面的（3.5）式。
-----------------------------------------------------------------

泛函J的极值条件为其变分为零，即
(3.2)
在此由于Lagrangian乘数λ不是自变量，我们不需要对其变分。基于变分δx的任意性，要求
(3.3)
由此可得
(3.4)
但是这里λ也是待定值，我们需要用迭代方法解此方程。
* 假设Lagrangian乘数的初值  ，n=0
* 带入上式求得
* 判断接触条件  是否满足。如不满足，则更新λ，
(3.5)
令n->n 1; 并返回上一步重求x。

在上述计算过程中，Lagrangian乘数λ从零出发逐渐递增，这是我们称之为递增Lagrangian乘数法的原因。同时，x的值从  出发逐渐递减，当  接近于0时计算结束。

可以看出，这里给出的x的初值与上述罚函数法给出的公式相同。因此，如果  的值足够大，只需一步计算就可以使  趋于0，如此递增Lagrangian乘数法将退化为罚函数法。但是在递增Lagrangian乘数法中罚函数  一般不是一个很大的数值，其原因将在下节讨论。

1.4 内点法

内点法是从可行域（即未接触区域g>0）内某一初始内点出发，在可行域内求极值的方法。其泛函J为

(4.1)

这里称为围墙（barrier）函数，因为该函数在未接触区域为一正值，但当其靠近接触面时急剧增大，趋近接触面时，其值趋向于正无穷大。它像围墙一样阻止迭代点越出约束边界。围墙函数可以有不同形式，如 . 另外系数r>0的作用类似于前述罚函数。但在此它是个待定变数。在下面（4.3）式中可以看到，当r趋于零时，泛函J趋于极值。

泛函J的极值条件为其变分为零，即

(4.2)

基于变分δx的任意性，要求

(4.3) &

在这里我们导入了一个新变量 . 可以看出该式与方程（2.3）非常类似，当r=0时两者等价。该方程组有三个变量x，r和λ, 需要用迭代方法求解。

* 取初值r〉0，递减系数0<γ<1，

* 带入(4.3)式求得x和λ. 注意到这是一个非线性方程组，需要将其线性化

(4.4) &

采用Newton迭代法求解。从此联立方程中可以直接求得x和λ的增量（primal-dual法），也可仅以x为自变量，靠迭代求解λ的增量（primal算法）。

* 令 , 回到上一步更新x和λ，直到r趋于0。

2. 算法评价

2.1 算法评价基准

评价一个数值算法的基本判据有

* 精度
* 收敛性
* 计算量
* 内存使用量
* 并行算法的相容性

对于有限元计算来说，其大部分计算时间消耗于联立方程的求解计算。求解联立方程大致可分为两类

* 直接法

* 以CG法为代表的迭代法

其比较如下

从上表可以看出，对于大型计算而言，迭代法是较好的选择。在本文中将考虑各种算法对迭代法的收敛性的影响。如上所述，迭代法的收敛性依赖于问题类型，其具体的指标则是联立方程A的条件数。简单地，条件数可以理解为矩阵对角项的最大值与最小值的比

条件数增大则迭代法收敛性变差，这意味着迭代次数和计算时间的增加甚至计算发散。

另外，我们希望方程是正定的，主对角集中的，这些都可以提高迭代法求解器的收敛性。

2.2 算法评价

罚函数法:

缺点：近似解，要求极大罚函数值，使方程组条件数变差。

优点：计算量不大，实现最简单。

Lagrangian乘数法:

缺点：增加了计算变量（计算量增加），方程性能变差（参见式（2.3）,可以看出其导入了零对角项，该方程变为非正定方程）。

优点：精确解。

递增Lagrangian乘数法:

缺点：需要迭代求解Lagrangian乘数（计算量增加）。

优点：回避了计算变量的增加；选用合适的罚函数值时，可以回避或减缓方程组条件数的恶化；对于接触问题解析来说，可以利用此方法把非对称的接触刚性矩阵对称化，大幅度节省内存和计算时间。

内点法:

缺点：不能处理初始点已经位于接触面内的问题；采用primal-dual算法增加了计算变量；需要迭代求解。

优点：内点法应用于接触问题解析主要是其可以回避解析过程中，接触点—〉非接触-〉接触，即所谓接触active set的变化问题，因为该方法对接触面附近的所有点都加上了约束。因此可以期待其提高计算的收敛性。

3 递增Lagrangian乘数法的实装方法（待续）

附注

*1 Augmented Lagrangian的译法有多种。因为augment的原意为 To make(something already developed or well underway) greater,as in size,extent,or quantity （http://www.thefreedictionary.com/augment），Make (something) greater by adding to it; increase (http://www.oxforddictionaries.com/definition/english/augment). 而Augmented Lagrangian法中，Lagrangian乘数的值正是逐渐增大的，因此提出该译法。

*2 接触条件一般在接触点为零点，接触面法线方向为正向的局部坐标系下表述。由于本例题只有一个自由度，本文简单地取局部坐标系为全局坐标。

*3 由于习惯的不同，(3.1)式也可以写为。这时得到的λ与(3.1)式相比差一个负号。