有限元接触解析算法（之二）
3 渐增Lagrangian乘数法/内点法的实装方法
前述第1节讲述的四种方法中，罚函数法，Lagrangian乘数法的解可以从所得方程直接解出。而递增Lagrangian乘数法和内点法需要使用迭代法求解。
在本节中，将以一自由度的接触问题为例，又易到难地讲解采用递增Lagrangian乘数法和内点法的接触算法。由于接触问题的解析一般用增分法求解，而在第一节中我们给出的是全量公式。在这一节中给出的算法是基于增分法的，因此第一节中给出的公式也要相应改写。
为显示算法实现，从本节开始将添加一些用python语言编写的程序附件。在之后的算法实现中。将大量使用迭代计算。为了明示迭代步数，本文将以上下标n，k表示各迭代步。
3.1 考虑接触非线性的基于递增Lagrangian乘数法的实装方法
图1 一自由度的接触问题S0: 初始状态的求解
1） 初期接触状态的判定:无接触，由（0.4）式得到物体m的位移量x0=F/K
2) 接触状态的判定: 如g=g0 x0>=0计算结束。否则判定进入接触，继续下一步计算。
以下进入递增Lagrangian乘数法求解。在这里我们令罚函数 ,收束判定值TOL=1.e-10
在此，对应于增分解析，我们将(3.1)改写为
J=1/2K(x x0)2−F(x x0) λ(x g) 1/2εN(x g)2=0
在这里x表示下一步的位移。相应地，式(3.4)改变为
x=(F−εNg0−λ−Kx0)/(K εN)
S1: k=0,λ0=0,x0=0
S2: if :g xk<=TOL
goto S3
else
xk 1=(F−εNg0−λ−Kx0)/(K εN)
λk 1=λk εN∗(g xk 1)
k=k 1; goto S2
S3: end
计算源码在此（好像没有百度账号不能上传）。
计算过程如下。
iter, disp incr , gap , Lagrangian
0, 0.000000000, -1.0000000e-01, 0.000000000
1, 0.090909091, -9.0909091e-03, -9.090909091
2, 0.099173554, -8.2644628e-04, -9.917355372
3, 0.099924869, -7.5131480e-05, -9.992486852
4, 0.099993170, -6.8301346e-06, -9.999316987
5, 0.099999379, -6.2092132e-07, -9.999937908
6, 0.099999944, -5.6447393e-08, -9.999994355
7, 0.099999995, -5.1315812e-09, -9.999999487
8, 0.100000000, -4.6650737e-10, -9.999999953
9, 0.100000000, -4.2409770e-11, -9.999999996
Displacement= -0.100000000042
Final gap= -4.24097701401e-11
Contact force= -9.99999999576
3.2 考虑接触非线性的的基于内点法的实装方法
这里的算法在（1.4）节给出的primal-dual算法。计算源码在此（好像没有百度账号不能上传）。
计算过程如下。
iter, displacement, Lagrangian, barrier parameter
1, -0.038196601, 16.180339887, 1.0000000e 00
2, -0.075838015, 12.416198487, 3.0000000e-01
3, -0.091690481, 10.830951895, 9.0000000e-02
4, -0.097369211, 10.263078947, 2.7000000e-02
5, -0.099196457, 10.080354318, 8.1000000e-03
6, -0.099757588, 10.024241236, 2.4300000e-03
7, -0.099927153, 10.007284693, 7.2900000e-04
8, -0.099978135, 10.002186522, 2.1870000e-04
9, -0.099993439, 10.000656057, 6.5610000e-05
10, -0.099998032, 10.000196826, 1.9683000e-05
11, -0.099999410, 10.000059049, 5.9049000e-06
12, -0.099999823, 10.000017715, 1.7714700e-06
13, -0.099999947, 10.000005314, 5.3144100e-07
14, -0.099999984, 10.000001594, 1.5943230e-07
15, -0.099999995, 10.000000478, 4.7829690e-08
16, -0.099999999, 10.000000143, 1.4348907e-08
17, -0.100000000, 10.000000043, 4.3046721e-09
18, -0.100000000, 10.000000013, 1.2914016e-09
19, -0.100000000, 10.000000004, 3.8742049e-10
20, -0.100000000, 10.000000001, 1.1622615e-10
Displacement= -0.0999999999884
Final gap= 1.16226195335e-11
Contact force= 10.0000000012
4. 考虑材料非线性和接触非线性的基于递增Lagrangian乘数法的实装方法
在前述算例中我们只考虑了接触非线性，而在实际应用中，接触非线性往往和材料非线性，几何非线性等并存。这一节将考虑非线性弹簧，其弹性系数K=100-1/30x。以此演示和其他非线性现象共存时接触解析的基于递增Lagrangian乘数法的实装方法。
此时的泛函J的变分为
此时的泛函J的变分为
J=1/2K02−1/30x3 (λ 1/2εNx)x−Fx
δJ=(K0−0.1x)xδx (λ εNx)δx−(F−εNg0)δx=0
基于变分δx的任意性，可得
(K−0.1x εN)x λ−(F εNg0)=0
这是在这一节要解的方程。采用Newton迭代法，上式的增分形式为
Δxn 1=(F−(K0−0.1xn)xn−λ−εN(xn g0))/(K0−0.2xn εN)
S0: 初始状态的求解
1） 初期接触状态的判定:无接触，由(K 0.1x)x λ=F求解物体m的初始位移量x0。由于该式非线性，我们需要使用迭代计算。
2) 接触状态的判定: 如 计算结束。否则判定进入接触，继续下一步计算。
以下进入递增Lagrangian乘数法求解。在此我们使用两重迭代(nested augmented lagrangian algorithm)，分别处理材料非线性和接触非线性
S1: k=0,λ0=0,x00=0
S2: 外部循环，更新计算Lagrangian乘数
if g xk<=TOL
goto S5
else
goto S3
S3 内部循环，Lagrangian乘数不变，更新计算位移增量
n = 0
df=F−εNg0−λ−(K0−0.1xnk)xnk
if df<TOL
goto S4
else
Δxn 1k 1=df/(K−0.2xnk εN)
xn 1k 1=xnk 1 Δxn 1k 1
n = n 1
goto S3
S4 λk 1=λk εN∗(g xk 1)
k=k 1; goto S2
S5 end
计算源码在此（好像没有百度账号不能上传）。
计算过程如下。
Initial displacement and contact gap -0.196152422706 0.1
iter, disp incr , gap , Lagrangian
0, -0.196152423, -9.6152423e-02, 0.000000000
internal loop, disp= 0.0870619064302
internal loop, disp= 6.89056837473e-07
1, -0.109089827, -9.0898272e-03, -9.089827219
internal loop, disp= 0.0082633153912
internal loop, disp= 6.20737540768e-09
2, -0.100826506, -8.2650562e-04, -9.916332840
internal loop, disp= 0.000751354971978
internal loop, disp= 5.13203608441e-11
3, -0.100075151, -7.5150597e-05, -9.991483437
internal loop, disp= 6.83174816365e-05
internal loop, disp= 4.24294993083e-13
4, -0.100006833, -6.8331151e-06, -9.998316552
internal loop, disp= 6.21180988284e-06
5, -0.100000621, -6.2130523e-07, -9.998937857
internal loop, disp= 5.64812673181e-07
6, -0.100000056, -5.6492560e-08, -9.998994350
internal loop, disp= 5.13559388072e-08
7, -0.100000005, -5.1366210e-09, -9.998999486
internal loop, disp= 4.66957053199e-09
8, -0.100000000, -4.6705044e-10, -9.998999953
internal loop, disp= 4.24583585173e-10
9, -0.100000000, -4.2466849e-11, -9.998999996
Displacement= -0.100000000042
Final gap= -4.24668494814e-11
Contact force= -9.99899999575
5. 摩擦接触问题的求解方法
图2 一自由度的摩擦接触问题下面将以图2的问题为例，讨论解析摩擦接触约束的方法。与上节相同，本节也采用弹性系数K=100 1/10x的非线性弹簧。
在这里我们采用库仑摩擦定理
(5.1) x=0;iff<μ∗mg
:f<=−μ∗mg∗signx˙
上面两式分别对应于是为滑动摩擦（f<μ∗mg， signx˙ 表示运动方向）和，和粘结摩擦状态(x=0)。
此时的泛函J为*1
(5.2) 粘结摩擦状态 J=1/2K0x2 0.1x3−Fx λx 1/2εTx2
滑动摩擦状态下: J=1/2K0x2 0.1x3−Fx μ∗mgx
上式中的 εT 是一个大数值的惩罚因子。用于迫使粘结摩擦状态的位移近似于零。但是如前所述，取大数对解析的收敛性不利。实际上各大商业软件很鸡贼，大都允许取较小的数。另外，也可以把εT当成augmented Lagrangian处理。
在这里，我们使用递增Lagrangian乘数法处理条件x=0;if f<μ∗mg，该泛函与式（3.1）相同。滑动摩擦状态下的泛函则对应于式（0.1）。
由于要考虑粘结和滑动的状态变化，其解法要相应地复杂一些。这一问题的解法如下:
S1: k=0,λ0=0,x00=0
S2: 外部循环，更新计算Lagrangian乘数
if Δx<=TOL
goto S5
else
goto S3
S3 内部循环，Lagrangian乘数不变，更新计算位移增量
n=0
fric=λkn
if 粘结摩擦
: df=F−λk−εTx−(K0 0.1xk)xk
: Δxn 1=df/(K0 0.2xn εT)
: fric=fric εTx
if fric>μ∗mg
fric=μ∗mg 进入滑动摩擦
else
: df=F−fric−(K0 0.1xk)xk
: Δxn 1=df/(K0 0.2xn)
fric=μ∗mg
x=x Δx
if df<TOL
goto S4
else
n=n 1; goto S3
S4 if 粘结摩擦
: λk 1=λk εT∗xk 1
else
: λk 1=μ∗mg
k=k 1; goto S2
S5 end
在这里我们显示的是一维问题。在2维，3维计算中，我们需要同时考虑3.3节所示不可贯入条件和本节所示摩擦条件，此时我们需要结合3.3和3.4节的内容。其详细将在后续博文中论及。
为考虑粘结摩擦到滑动摩擦变化的处理方法，所付计算源码在此（好像没有百度账号不能上传）考虑了一个3步加载的问题，其计算结果如下:
-----Step: 1
External Force= 10.0 0.0
Iter, Displacement , Lagrangian
internal loop, disp= 0.00909090909091 0.00909090909091
internal loop, disp= 0.00909090157777 0.00909090157777
internal loop, disp= 0.00909090157777 0.00909090157777
1, 0.009090902, 9.090901578
internal loop, disp= 0.000826459258203 0.000826459258203
internal loop, disp= 0.000826453049022 0.000826453049022
internal loop, disp= 0.000826453049022 0.000826453049022
2, 0.000826453, 9.917354627
internal loop, disp= 7.51322082622e-05 7.51322082622e-05
internal loop, disp= 7.51321569456e-05 7.51321569456e-05
internal loop, disp= 7.51321569456e-05 7.51321569456e-05
3, 0.000075132, 9.992486784
internal loop, disp= 6.830197019e-06 6.830197019e-06
internal loop, disp= 6.83019659489e-06 6.83019659489e-06
internal loop, disp= 6.83019659489e-06 6.83019659489e-06
4, 0.000006830, 9.999316980
internal loop, disp= 6.20926970885e-07 6.20926970885e-07
internal loop, disp= 6.20926967378e-07 6.20926967378e-07
5, 0.000000621, 9.999937907
internal loop, disp= 5.64479061917e-08 5.64479061917e-08
internal loop, disp= 5.64479061623e-08 5.64479061623e-08
6, 0.000000056, 9.999994355
internal loop, disp= 5.13162783284e-09 5.13162783284e-09
internal loop, disp= 5.13162783284e-09 5.13162783284e-09
7, 0.000000005, 9.999999487
internal loop, disp= 4.66511621933e-10 4.66511621933e-10
internal loop, disp= 4.66511620318e-10 4.66511620318e-10
8, 0.000000000, 9.999999953
internal loop, disp= 4.24101488013e-11 4.24101488013e-11
internal loop, disp= 4.24101488013e-11 4.24101488013e-11
9, 0.000000000, 9.999999996
internal loop, disp= 3.85546854484e-12 3.85546854484e-12
internal loop, disp= 3.85546692997e-12 3.85546692997e-12
10, 0.000000000, 10.000000000
Total Displacement & Displacement of current step= 3.85546692997e-12 3.85546692997e-12
frictional force= 9.99999999961
Spring force= 3.85546692997e-10
Reaction force= 10.0
-----Step: 2
External Force= 20.0 9.99999999961
Iter, Displacement , Lagrangian
From stick to slip state: 19.0909090905 > 10.0
internal loop, disp= 0.00909090909476 0.00909090909091
internal loop, disp= 0.0999982644944 0.0999982644905
internal loop, disp= 0.0999900019996 0.0999900019957
internal loop, disp= 0.0999900019995 0.0999900019956
11, 0.099990002, 10.000000000
internal loop, disp= 0.0999900019995 0.0999900019956
12, 0.099990002, 10.000000000
Total Displacement & Displacement of current step= 0.0999900019995 0.0999900019956
frictional force= 10.0
Spring force= 10.0
Reaction force= 20.0
-----Step: 3
External Force= 10.0 10.0
Iter, Displacement , Lagrangian
From slip to stick state: -89.980005998 < 10.0
internal loop, disp= 9.99600149944e-06 -0.099980005998
internal loop, disp= 0.0909000016525 -0.00909000034697
internal loop, disp= 0.0908992506656 -0.00909075133393
internal loop, disp= 0.0908992506656 -0.00909075133393
13, 0.090899251, -9.090751334
internal loop, disp= 0.0991634334768 -0.000826568522726
internal loop, disp= 0.0991634272681 -0.000826574731406
internal loop, disp= 0.0991634272681 -0.000826574731406
14, 0.099163427, -9.917326065
internal loop, disp= 0.0999148452942 -7.51567052527e-05
internal loop, disp= 0.0999148452429 -7.51567565817e-05
internal loop, disp= 0.0999148452429 -7.51567565817e-05
15, 0.099914845, -9.992482822
internal loop, disp= 0.0999831683259 -6.83367359656e-06
internal loop, disp= 0.0999831683255 -6.83367402093e-06
internal loop, disp= 0.0999831683255 -6.83367402093e-06
16, 0.099983168, -9.999316496
internal loop, disp= 0.0999893806435 -6.21356025034e-07
internal loop, disp= 0.0999893806435 -6.21356028544e-07
17, 0.099989381, -9.999937852
internal loop, disp= 0.0999899455023 -5.64971807554e-08
internal loop, disp= 0.0999899455023 -5.64971807833e-08
18, 0.099989946, -9.999994349
internal loop, disp= 0.0999899968625 -5.13704107356e-09
internal loop, disp= 0.0999899968625 -5.1370410737e-09
19, 0.099989997, -9.999999486
internal loop, disp= 0.0999900015324 -4.67088634118e-10
internal loop, disp= 0.0999900015324 -4.67088632216e-10
20, 0.099990002, -9.999999953
internal loop, disp= 0.099990001957 -4.24703217975e-11
internal loop, disp= 0.099990001957 -4.24703219933e-11
21, 0.099990002, -9.999999996
internal loop, disp= 0.0999900019956 -3.86164012211e-12
internal loop, disp= 0.0999900019956 -3.8616413703e-12
22, 0.099990002, -10.000000000
Total Displacement & Displacement of current step= 0.0999900019956 -3.8616413703e-12
frictional force= 3.86238596661e-10
Spring force= 9.99999999961
Reaction force= 10.0
第一加载步处于粘结摩擦状态，可以看出此时的位移为零，但出现摩擦力。第二加载步则进入滑动摩擦，可以看出明显的位移，同时摩擦力满足库仑摩擦定理。在第3加载步，物体滑动方向逆转，再次进入粘结摩擦状态。在这一步，位移增分为零，摩擦力回零。
附注
*1 对于摩擦状态变化的问题，没有统一的泛函表达式，需要使用虚位移原理公式。这一内容将在后续博文中论及。
