有限元接触解析算法（之二）

罚函数法，Lagrangian乘数法，递增Lagrangian乘数法， 内点法

3  渐增Lagrangian乘数法/内点法的实装方法

前述第1节讲述的四种方法中，罚函数法，Lagrangian乘数法的解可以从所得方程直接解出。而递增Lagrangian乘数法和内点法需要使用迭代法求解。

在本节中，将以一自由度的接触问题为例，又易到难地讲解采用递增Lagrangian乘数法和内点法的接触算法。由于接触问题的解析一般用增分法求解，而在第一节中我们给出的是全量公式。在这一节中给出的算法是基于增分法的，因此第一节中给出的公式也要相应改写。

为显示算法实现，从本节开始将添加一些用python语言编写的程序附件。在之后的算法实现中。将大量使用迭代计算。为了明示迭代步数，本文将以上下标n，k表示各迭代步。

3.1  考虑接触非线性的基于递增Lagrangian乘数法的实装方法

图1 一自由度的接触问题

S0:  初始状态的求解

  1） 初期接触状态的判定:无接触，由（0.4）式得到物体m的位移量$x_0=F/K$

  2)  接触状态的判定: 如$g=g_0+x_0>=0$计算结束。否则判定进入接触，继续下一步计算。

  以下进入递增Lagrangian乘数法求解。在这里我们令罚函数$\varepsilon_N=10000$,收束判定值TOL=1.e-10

在此，对应于增分解析，我们将(3.1)改写为
  $J=1/2K(x+x_0)^2 -F(x+x_0) +\lambda (x+g)+1/2\varepsilon_N (x+g)^2=0$
在这里x表示下一步的位移。相应地，式(3.4)改变为
  $ x=(F-\varepsilon_Ng_0-\lambda-Kx_0)/(K+\varepsilon_N)  $

S1: k=0, $\lambda_0=0$, $x_0=0$
S2: if $g+x_k<=TOL$
         goto S3
      else
         $ x_{k+1}=(F-\varepsilon_Ng_0-\lambda-Kx_0)/(K+\varepsilon_N)  $
         $\lambda_{k+1}=\lambda_{k}+ \varepsilon_N *(g+x_{k+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

3.3  考虑材料非线性和接触非线性的基于递增Lagrangian乘数法的实装方法

在前述算例中我们只考虑了接触非线性，而在实际应用中，接触非线性往往和材料非线性，几何非线性等并存。这一节将考虑非线性弹簧，其弹性系数K=100-1/30x。以此演示和其他非线性现象共存时接触解析的基于递增Lagrangian乘数法的实装方法。

此时的泛函J的变分为
$J=1/2K_0^2-1/30 x^3+ (\lambda+1/2\varepsilon_Nx) x -F x$
$\delta J=(K_0-0.1x) x\delta x+ (\lambda+\varepsilon_Nx)\delta x -(F-\varepsilon_Ng_0) \delta x=0$
基于变分$\delta x$的任意性，可得
$(K-0.1x+\varepsilon_N)x+\lambda-(F+\varepsilon_Ng_0)=0$
这是在这一节要解的方程。采用Newton迭代法，上式的增分形式为
$\Delta x^{n+1}=(F-(K_0-0.1x^n)x^n-\lambda-\varepsilon_N(x^n+g_0))/(K_0-0.2x^n+\varepsilon_N)$

S0:  初始状态的求解

  1） 初期接触状态的判定:无接触，由$(K+0.1x)x+\lambda=F$求解物体m的初始位移量$x_0$。由于该式非线性，我们需要使用迭代计算。

  2)  接触状态的判定: 如$g=g_0+x_0>=0$计算结束。否则判定进入接触，继续下一步计算。

  以下进入递增Lagrangian乘数法求解。在此我们使用两重迭代(nested augmented lagrangian algorithm)，分别处理材料非线性和接触非线性
S1: k=0, $\lambda_0=0$, $x_0^0=0$
S2: 外部循环，更新计算Lagrangian乘数
      if $g+x_k<=TOL$
         goto S5
      else
         goto S3
S3   内部循环，Lagrangian乘数不变，更新计算位移增量
       n = 0
       $df =F-\varepsilon_Ng_0-\lambda-(K_0-0.1 x_k^n) x_k^n$
       if   df<TOL
          goto S4
       else
          $ \Delta x_{k+1}^{n+1}=df/(K-0.2x_k^n+\varepsilon_N)  $
           $x_{k+1}^{n+1}= x_{k+1}^n + \Delta x_{k+1}^{n+1}$
          n = n+1
          goto S3
S4   $\lambda_{k+1}=\lambda_{k}+ \varepsilon_N *(g+x_{k+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

3.4  摩擦接触问题的求解方法

     

图2 一自由度的摩擦接触问题

下面将以图2的问题为例，讨论解析摩擦接触约束的方法。与上节相同，本节也采用弹性系数K=100+1/10x的非线性弹簧。

在这里我们采用库仑摩擦定理
(4.1)$\begin{eqnarray}
f<=-\mu*mg*sign \dot x \\
x=0;\    if\   f<\mu*mg
\end{eqnarray}$
在$f<\mu*mg$是为粘结摩擦状态，否则称为滑动摩擦状态。

此时的泛函J为*1
(4.2)$\begin{eqnarray}
在粘结摩擦状态下:  J=1/2K_0 x^2 +0.1/3 x^3 -F x +\lambda x+1/2\varepsilon_T x^2 \\
在滑动摩擦状态下:  J=1/2K_0 x^2 +0.1/3 x^3 -F x +\mu*mg x
\end{eqnarray}$
在这里，我们使用递增Lagrangian乘数法处理条件$x=0;\    if\   f<\mu*mg$，该泛函与式（3.1）相同。滑动摩擦状态下的泛函则对应于式（0.1）。

由于要考虑粘结和滑动的状态变化，其解法要相应地复杂一些。这一问题的解法如下:

S1: k=0, $\lambda_0=0$, $x_0^0=0$
S2: 外部循环，更新计算Lagrangian乘数
     if $\Delta x<=TOL$
         goto S5
      else
         goto S3
S3  内部循环，Lagrangian乘数不变，更新计算位移增量
      n=0
      $fric=\lambda_k^n$
       if 粘结摩擦
          $df = F-\lambda_k-\varepsilon_T x -(K_0+0.1x_k)x_k$
          $\Delta x_{n+1} =df/(K_0+0.2x_n+\varepsilon_T)$
          $fric = fric+\varepsilon_T x$
          if $fric>\mu*mg$
              $fric = \mu*mg$进入滑动摩擦
       else
          $df = F-fric-(K_0+0.1x_k)x_k$
          $\Delta x_{n+1} =df/(K_0+0.2x_n)$
          $fric = \mu*mg$
       $ x=x+\Delta x$
        if df<TOL
           goto S4
       else
           n=n+1;  goto S3
S4   if 粘结摩擦
          $\lambda_{k+1}=\lambda_{k}+ \varepsilon_T *x_{k+1}$
      else
           $\lambda_{k+1}=\mu*mg$
      k=k+1;  goto S2
S3   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 对于摩擦状态变化的问题，没有统一的泛函表达式，需要使用虚位移原理公式。这一内容将在后续博文中论及。