有限元接触解析方法: 罚函数法，Lagrangian，递增Lagrangian乘数法*1， 内点法
0. 本文的用意和写法
该系列博文面对对有限元解析有所了解，但未接触过接触解析的初学者。本文以实用性为第一目标，试图在概念准确的前提下追求直观描述，特别强调给出某些实用算法的细节，并附有相应程序。由于这些算法细节都经过本人实装测试，其有效性是有一定的保证的。希望读完此系列博文可以自己动手组装自己的解析程序，并希望能够指出其中问题或错误，共同提高。
本文的解析对象是泛函的极值问题，这是本文的起点。对泛函的极值问题不了解的读者可以跳过泛函的变分运算步骤，着重于其后的物理概念和算法的讨论。
1. 考虑约束条件的最优化算法
泛函的极值问题的求解从数值计算方法的角度来说是最优化算法，而接触问题为其约束条件。罚函数法，Lagrangian乘数法，递增Lagrangian乘数法*1， 内点法(Interior point method, barrier method或内部罚函数法)都是标准的考虑约束条件的最优化算法。这一节将以图1的问题为例解释这几种解法。
图1 一自由度的接触问题*2图1所示问题如没有接触约束，其泛函J为
(0.1) J=1/2Kx2−Fx
这里我们将m的位移记为x，外力F<0。该泛函的极值条件为其变分为零，即
(0.2) δJ=Kxδx−Fδx=0
基于变分δx的任意性，要求
(0.3) Kx−F=0
即得
(0.4) x=F/K
在此我们假设 F/K 足够大使得 |x|>g0 ，但是接触条件要求物体m不能进入接触面, 即 g=g0 x≥0 。如 g<0 , 可以用下述方法导入接触约束。
1.1 罚函数法
采用罚函数法考虑接触约束，其泛函J为
(1.1) J=1/2Kx2−Fx 1/2εN(x g0)2
这里 εN>0 为罚函数。该泛函的极值条件为其变分为零，即
(1.2) δJ=(K εN)xδx−(F−εNg0)δx=0
基于变分δx的任意性，要求
(1.3) (K εN)x−(F−εNg0)=0
即得
(1.4) x=(F−εNg0)/(K εN)
当罚函数 εN 足够大时，从该公式可以得到逼近于精确解 x=−g0 的近似解。
1.2 Lagrangian乘数法
采用Lagrangian乘数法考虑接触约束，其泛函J为
(2.1) J=1/2Kx2−Fx λ(x g0)
这里λ为Lagrangian乘数,这是一个待定自变量。泛函J的极值条件为其变分为零，即
(2.2) δJ=(Kx λ−F)δx δλ(x g0)=0
基于变分δx和δλ的任意性，要求
(2.3) Kx λ=F & x=−g0
即得
(2.4) x=−g0;λ=F Kg0
这是该问题的精确解。x为m的位置，λ则为约束反力。
1.3 递增Lagrangian乘数法
采用递增Lagrangian乘数法考虑接触约束，其泛函J为*3
(3.1) J=1/2Kx2−Fx λ(x g0) 1/2εN(x g0)2
注意到该泛函既包括Lagrangian乘数λ也包括罚函数 εN>0 。但与前述Lagrangian乘数法和罚函数法不同的是，在该算法中
* Lagrangian乘数λ不是自变量
* 罚函数 εN 可以不是很大的数值
-----------------------------------------------------------------
注： 递增Lagrangian乘数法可以看成是罚函数法或Lagrangian乘数法的变种。如从Lagrangian乘数法出发，我们将（2.1）改写为
J=1/2Kx2−Fx λ(x g0) (λ∗−λ)(x g0)
这里，我们记 λ∗ 为Lagrangian乘数的准确值。考虑 λ∗−λ 为Lagrangian乘数的增分并将此式与（3.1）比较，可以看出
εN(x g0)≈λ∗−λ
如此Lagrangian乘数的准确值可以如下逼近
λ∗≈λ εN(x g0)
对应于下面的（3.5）式。
-----------------------------------------------------------------
泛函J的极值条件为其变分为零，即
(3.2) δJ=Kxδx (λ εNx)δx−(F−εNg0)δx=0
在此由于Lagrangian乘数λ不是自变量，我们不需要对其变分。基于变分δx的任意性，要求
(3.3) (K εN)x λ−(F εNg0)=0
由此可得
(3.4) x=(F−εNg0−λ)/(K εN)
但是这里λ也是待定值，我们需要用迭代方法解此方程。
* 假设Lagrangian乘数的初值 λ0=0 ，n=0
* 带入上式求得 xn 1
* 判断接触条件 g0 x>=0 是否满足。如不满足，则更新λ，
(3.5) λn 1≈λn εN(g0 x)
令n->n 1; 并返回上一步重求x。
在上述计算过程中，Lagrangian乘数λ从零出发逐渐递增，这是我们称之为递增Lagrangian乘数法的原因。同时，x的值从 x=(F−εNg0)/(K εN) 出发逐渐递减，当 g0 x 接近于0时计算结束。
可以看出，这里给出的x的初值与上述罚函数法给出的公式相同。因此，如果 εN 的值足够大，只需一步计算就可以使 g0 x 趋于0，如此递增Lagrangian乘数法将退化为罚函数法。但是在递增Lagrangian乘数法中罚函数 εN 一般不是一个很大的数值，其原因将在下节讨论。
1.4 内点法
内点法是从可行域（即未接触区域g>0）内某一初始内点出发，在可行域内求极值的方法。其泛函J为
(4.1) J=1/2Kx2−Fx−rlog(g0 x)=0
这里 G(g)=log(x g0) 称为围墙（barrier）函数，因为该函数在未接触区域为一正值，但当其靠近接触面时急剧增大，趋近接触面时，其值趋向于正无穷大。它像围墙一样阻止迭代点越出约束边界。围墙函数可以有不同形式，如 G(g)=1/(g0 x) . 另外系数r>0的作用类似于前述罚函数。但在此它是个待定变数。在下面（4.3）式中可以看到，当r趋于零时，泛函J趋于极值。
泛函J的极值条件为其变分为零，即
(4.2) δJ=Kxδx−(F r/(g0 x))δx
基于变分δx的任意性，要求
(4.3) Kx−λ−F=0 & λ(g0 x)=r
在这里我们导入了一个新变量 λ=r/(g0 x) . 可以看出该式与方程（2.3）非常类似，当r=0时两者等价。该方程组有三个变量x，r和λ, 需要用迭代方法求解。
* 取初值r〉0，递减系数0<γ<1， λ0=0,x0=0
* 带入(4.3)式求得x和λ. 注意到这是一个非线性方程组，需要将其线性化
(4.4) KΔx−Δλ=F−Kx λ & λΔx (g0 x)Δλ=r−λ(g0 x)
采用Newton迭代法求解。从此联立方程中可以直接求得x和λ的增量（primal-dual法），也可仅以x为自变量，靠迭代求解λ的增量（primal算法）。
* 令 r=rγ , 回到上一步更新x和λ，直到r趋于0。
2. 算法评价
2.1 算法评价基准
评价一个数值算法的基本判据有
* 精度
* 收敛性
* 计算量
* 内存使用量
* 并行算法的相容性
对于有限元计算来说，其大部分计算时间消耗于联立方程的求解计算。求解联立方程大致可分为两类
* 直接法
* 以CG法为代表的迭代法
其比较如下
从上表可以看出，对于大型计算而言，迭代法是较好的选择。在本文中将考虑各种算法对迭代法的收敛性的影响。如上所述，迭代法的收敛性依赖于问题类型，其具体的指标则是联立方程A的条件数。简单地，条件数可以理解为矩阵对角项的最大值与最小值的比
cond(A)=max|Adiag|/min|Adiag|
条件数增大则迭代法收敛性变差，这意味着迭代次数和计算时间的增加甚至计算发散。
另外，我们希望方程是正定的，主对角集中的，这些都可以提高迭代法求解器的收敛性。
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)式也可以写为 J=1/2Kx2−Fx−λ(x g0) 1/2εN(x g0)2 。这时得到的λ与(3.1)式相比差一个负号。
