有限元中约束方程的施加方法：拉格朗日乘子法和罚方法

2021-02-26 2258

最近经常有人问到有限元编程中关于约束方程应该如何处理的问题，于是决定借助本文从简单的数值实现上来探讨通用的含有线性约束方程时的计算方法，具体介绍最常用的两种方法：拉格朗日乘子法和罚方法。

关于方法的原理可以在任何一本数值计算或者有限元相关的书上找到，这里仅仅介绍如何数值实现。

首先让我们考虑一个线性方程组：

$\bf{KU=F} \$

${\bf{K}}=\begin{bmatrix} 4.5 & 1.2 &-3.3\ 1.2 & 6.0 & 1.9 \ -3.3 & 1.9 & 4.7 \end{bmatrix} ,\quad {\bf{F}}=\begin{Bmatrix} 2.1\0.0\-0.5 \end{Bmatrix} ,\quad {\bf{U}}= \begin{Bmatrix} u_1\ u_2\ u_3 \end{Bmatrix} \$

这样一个方程显然有解：

${\bf{U}}={\bf{K}}^{-1}{\bf{F}}= \begin{Bmatrix} 1.72769\ -0.798162\ 1.42934 \end{Bmatrix} \$

该方程可以类比有限元中组装完成后的总刚和荷载列向量并通过消去法施加约束后的方程。如果，因为某些原因，导致我们希望在模型上施加一些额外的自由度之间的约束，例如这里我们想要对 $u_1$ 和 $u_2$ 两个自由度施加一组约束方程(1)为：

$-1.2u_1+0.5u_2=-2.0 \$

这时候就可以采用拉格朗日乘子法，或者罚方法来施加约束。此时求出的自由度是满足约束方程条件下的势能泛函的极值点，也被称为拉格朗日条件极值。

拉格朗日乘子法

先将约束方程中缺少的自由度补充完整，得：

$-1.2u_1+0.5u_2+0.0u_3=-2.0 \$

写出约束方程的系数矩阵：

${\bf{Q}}=\begin{Bmatrix} -1.2\ 0.5\ 0.0 \end{Bmatrix} \ f=-2.0$

然后可以得到考虑约束后的刚阵和右端项：

${\bf K}_\lambda=\begin{bmatrix} {\bf K}&{\bf Q} \ {\bf Q}^{\rm T}&{\bf 0} \end{bmatrix} =\begin{bmatrix} 4.5&1.2&-3.3&-1.2\ 1.2&6.0&1.9&0.5\ -3.3&1.9&4.7&0.0\ -1.2&0.5&0.0&0.0 \end{bmatrix} \ {\bf F}_\lambda=\begin{Bmatrix} {\bf F}\ f \end{Bmatrix} = \begin{Bmatrix} 2.1\ 0.0\ -0.5\ -2.0 \end{Bmatrix}$

此时可以解出：

$\begin{Bmatrix} u_1\ u_2\ u_3\ \lambda \end{Bmatrix} = {\bf K}_\lambda^{-1}{\bf F}_\lambda = \begin{Bmatrix} 1.40723\ -0.622654\ 1.13338\ -0.21235 \end{Bmatrix} \$

这个解是满足约束方程(1)的。其中拉格朗日乘子自由度的值表示约束内力。

罚方法

采用罚方法施加约束方程(1)，首先仍然写出系数矩阵：

${\bf{Q}}=\begin{Bmatrix} -1.2\ 0.5\ 0.0 \end{Bmatrix} \ f=-2.0$

将罚系数取为一个较大的数，一般可以取刚度最大元素的1e6倍，例如本例取为：

$p=10^7\$

则罚方法的刚度矩阵和右端项为：

${\bf K}_p={\bf K}+p{\bf Q}{\bf Q}^{\rm T}= \begin{bmatrix} 4.5 & 1.2 &-3.3\ 1.2 & 6.0 & 1.9 \ -3.3 & 1.9 & 4.7 \end{bmatrix}+10^7 \begin{bmatrix} 1.44 & -0.6 &0.0\ -0.6 & 0.25 & 0.0 \ 0.0 & 0.0 &0.0 \end{bmatrix} \ {\bf F}_p= {\bf F}+pf{\bf Q}= \begin{Bmatrix} 2.1\0.0\-0.5 \end{Bmatrix} + 10^7\times(-2.0) \begin{Bmatrix} -1.2\0.5\-0.0 \end{Bmatrix} \$

同样可以求得：

$\begin{Bmatrix} u_1\ u_2\ u_3 \end{Bmatrix} = {\bf K}_p^{-1}{\bf F}_p = \begin{Bmatrix} 1.40723\ -0.622654\ 1.13338 \end{Bmatrix} \$

罚方法得到的解只是近似的满足约束方程，这里贴出软件算出的罚方法的结果：

$u_1=1.4072275746568566 \ u_2=-0.6226538632934947$

其与约束方程的误差为：

$g=-1.2u_1+0.5u_2-(-2.0)=-2.1235\times10^{-8} \$

该误差乘以罚系数则得到拉格朗日乘子的值：

$\lambda=gp=-0.21235 \$

这里贴出采用mathematica软件计算的代码，首先是原始线性方程组的解：

Clear["`*"]
(*Linear Equations*)
K = {
   {4.5, 1.2, -3.3},
   {1.2, 6.0, 1.9},
   {-3.3, 1.9, 4.7}
   };
Det[K]
F = {
   {2.1},
   {0},
   {-0.5}
   };
U = LinearSolve[K, F];
Print["u1 = ", U[[1, 1]], "\nu2 = ", U[[2, 1]], "\nu1 = ", U[[3, 1]]];

然后是拉格朗日乘子法：

(*
constraint:
-1.2*u1+0.5*u2+0.0*u3=-2.0
*)
(*Lagrangian multiplier Method*)
Q = {
   {-1.2},
   {0.5},
   {0.0}
   };

K\[Lambda] = ArrayFlatten[{{K, Q}, {Transpose[Q], 0}}];

Print["K\[Lambda] = ", MatrixForm[K\[Lambda]]];

F\[Lambda] = ArrayFlatten[{{F}, {-2.0}}];

U = LinearSolve[K\[Lambda], F\[Lambda]];

Print["u1 = ", U[[1, 1]],
  "\nu2 = ", U[[2, 1]],
  "\nu1 = ", U[[3, 1]],
  "\n\[Lambda] = ", U[[4, 1]]];

Print["-1.2*u1+0.5*u2 = ",
  -1.2*U[[1, 1]] + 0.5*U[[2, 1]]];

最后是罚方法：

(*
constraint:
-1.2*u1+0.5*u2+0.0*u3=-2.0
*)
(*Penalty Method*)
penalty = 10^7;
Q = {
   {-1.2},
   {0.5},
   {0.0}
   };

Kp = K + Q.Transpose[Q]*penalty;

Print["K\[Lambda] = ", MatrixForm[Kp]];

Fp = F + Q*(-2.0)*penalty;

U = LinearSolve[Kp, Fp];

Print["u1 = ", U[[1, 1]],
  "\nu2 = ", U[[2, 1]],
  "\nu1 = ", U[[3, 1]]];

Print["-1.2*u1+0.5*u2 = ", -1.2*U[[1, 1]] + 0.5*U[[2, 1]]];

g = -1.2*U[[1, 1]] + 0.5*U[[2, 1]] - (-2);

Print["\[Lambda] = ", g*penalty];

添加评论

登录后可评论