Mortar算法介绍

简介

结构仿真中的绑定 ( Tie ) 希望主从面对应位置的位移相同，数学表达形式是：

\[\boldsymbol u_s(\boldsymbol x) - \boldsymbol u_m(\chi(\boldsymbol x)) =0 \]

其中 \(\chi(\boldsymbol x)\) 表示从面点在主面上的对应位置。

当主从节点重合时，Tie的约束方程可以直接写成 \(\boldsymbol d_i^s = \boldsymbol d_i^m\) 的形式。但是非匹配网格下，主从节点通常不重合，无法直接写出数学表达式。Mortar Method 是一种在主从面网格不匹配时，用弱积分形式建立位移约束的方法。

它在总势能中增加了界面约束项：

\[\Pi_c = \int_{\Gamma_s} \boldsymbol\lambda^T \left( \boldsymbol u_s-\boldsymbol u_m \right)d\Gamma \]

其中 (\boldsymbol\lambda) 的物理意义是界面牵引力。

对拉格朗日乘子取变分：

\[\delta_{\lambda}\Pi_c = \int_{\Gamma_s} \delta\boldsymbol\lambda^T \left( \boldsymbol u_s-\boldsymbol u_m \right)d\Gamma =0 \]

也就是：

\[ \int_{\Gamma_s} \delta\boldsymbol\lambda^T \boldsymbol g_u\,d\Gamma =0 \qquad \forall\,\delta\boldsymbol\lambda \]

它不再要求 (\boldsymbol g_u) 在每一点都严格为零，而是要求它对乘子空间中的所有测试函数都正交。

有限元推导

为了方便推导，我们暂时把三维向量看成一个标量分量。设从面有 (n_s) 个节点，主面有 (n_m) 个节点：

\[ u_s^h(\xi) = \sum_{j=1}^{n_s} N_j^s(\xi)d_j^s \]

主面的形函数需要在从面积分点对应的主面坐标 (\eta(\xi)) 上计算：

\[u_m^h(\eta(\xi)) = \sum_{k=1}^{n_m} N_k^m(\eta(\xi))d_k^m \]

拉格朗日乘子离散为：

\[\lambda^h(\xi) = \sum_{i=1}^{n_\lambda} \Phi_i(\xi)\Lambda_i \]

这里用 (\Lambda_i) 表示乘子节点值，避免把“nodal multiplier”与法向乘子 (\lambda_n) 混淆。

乘子变分为：

\[ \delta\lambda^h(\xi) = \sum_{i=1}^{n_\lambda} \Phi_i(\xi)\delta\Lambda_i \]

代入弱约束：

\[\int_{\Gamma_s} \left( \sum_i\Phi_i\delta\Lambda_i \right) \left[ \sum_jN_j^sd_j^s - \sum_kN_k^md_k^m \right]d\Gamma =0 \]

展开求和：

\[\sum_i\delta\Lambda_i \left[ \sum_j \left( \int_{\Gamma_s}\Phi_iN_j^s\,d\Gamma \right)d_j^s - \sum_k \left( \int_{\Gamma_s}\Phi_iN_k^m\,d\Gamma \right)d_k^m \right] =0 \]

由于每个 (\delta\Lambda_i) 都是任意且相互独立的，所以方括号必须逐行等于零：

\[\sum_jD_{ij}d_j^s - \sum_kC_{ik}d_k^m =0 \]

其中：

\[D_{ij} = \int_{\Gamma_s} \Phi_iN_j^s\,d\Gamma \\ C_{ik} = \int_{\Gamma_s} \Phi_iN_k^m(\eta(\xi))\,d\Gamma \]

写成矩阵形式就是：

\[D\boldsymbol d_s-C\boldsymbol d_m=0 \]

Dual Mortar --- 更高效的Mortar

如果直接选择：

\[ \Phi_i=N_i^s \]

那么 \(D_{ij} = \int_{\Gamma_s}N_i^sN_j^s\,d\Gamma\)​ 通常是非对角矩阵，求解时需面对复杂的鞍点问题或者求解大型稠密矩阵的逆，计算成本极高。为了解决这个问题，对偶 Mortar 法出现了。

对偶 Mortar 法的关键在于，它在从边界（Slave / Non-mortar side）上为 Lagrange 乘子（通常代表接触压力或交界面通量）构造了一组特殊的基函数 (\Phi {j})。这组基函数与该边界上有限元位移的标准形函数 (N) 满足双正交条件（Bi-orthogonal condition）： [1, 2, 3]

\[\int _{\Gamma }\Phi _{j}\cdot N_{i}\,d\Gamma =d_{j}\delta _{ji} \]

其中 \(\delta _{ji}\) 是 Kronecker delta 符号（当 i=j 时为1，否则为0），\(d_{j}\) 为缩放常数。 [1]

于是：

\[D= \begin{bmatrix} d_1&&\\ &d_2&\\ &&\ddots \end{bmatrix} \]

第 (i) 行约束变成：

\[d_i\,\boldsymbol d_i^s - \sum_kC_{ik}\boldsymbol d_k^m =0 \]

定义\(w_{ik}=\frac{C_{ik}}{d_i}\)，最终就是代码需要的多点约束：

\[\boxed{ \boldsymbol d_i^s = \sum_k w_{ik}\boldsymbol d_k^m } \]

附录

三维位移怎么表示

对于三维位移：

\[\boldsymbol d_j^s = \begin{bmatrix} d_{jx}^s\\ d_{jy}^s\\ d_{jz}^s \end{bmatrix} \]

位移插值矩阵实际是：

\[\boldsymbol N_s = \begin{bmatrix} N_1^sI_3 & N_2^sI_3 & \cdots & N_{n_s}^sI_3 \end{bmatrix} \]

如果三个方向使用相同标量形函数，那么完整矩阵为：

\[\left(D\otimes I_3\right)\boldsymbol d_s - \left(C\otimes I_3\right)\boldsymbol d_m =0 \]

代码通常只存标量 \(D_{ij}\)、\(C_{ik}\)，装配约束时再作用到 \(x,y,z\) 三个位移分量。