多尺度接触力学的迭代接触求解混合方法

《Computer Methods in Applied Mechanics and Engineering》：Iterative contact-resolving hybrid methods for multiscale contact mechanics

　　Eric T. Chung|Hyea Hyun Kim|Xiang Zhong
中国香港特别行政区沙田，香港中文大学数学系

摘要
对具有高对比系数的接触力学问题进行建模面临着重大的数学和计算挑战，尤其是在为混合公式实现强对称应力近似方面。由于接触问题本身具有非线性特征，那些将整个区域视为整体系统的传统方法往往会导致极高的全局复杂度。为解决这一问题，我们开发了一种迭代型接触求解混合方法，该方法将非线性接触约束限制在较小的子域内，而较大的子域则由线性系统描述。该系统运用了变分不等式理论、最小化原理以及罚函数方法。更重要的是，我们在双子域框架下提出了四种离散化类型，从在整个区域应用标准/混合有限元法，到在较大子域采用标准/混合多尺度方法，同时在较小子域应用标准/混合有限元法。通过采用多尺度简化技术，该方法避免了传统方法在较大区域中存在的过多自由度问题；而混合公式则能够直接计算应力、确保局部动量守恒，并有效防止几乎不可压缩材料的锁定现象。我们还为所有情况提供了收敛性分析及相应算法，并通过大量数值实验验证了这些方法的有效性。

引言
在高对比系数下的接触力学问题出现在许多工程和地球物理应用中，这类问题涉及机械性能差异显著的材料之间的相互作用。典型例子包括橡胶密封件与金属表面的接触、轮胎与路面的接触，以及不同岩层之间的地质断层。这类问题的特点是接触界面处材料参数（如杨氏模量、拉梅常数）存在巨大差异，从而导致复杂的局部变形模式，给数值模拟带来很大挑战。高对比系数加上非线性接触条件带来了诸多难题，包括解的不稳定性和边界层问题，这需要强大的离散化和求解技术。Kikuchi和Oden[1]建立了基于变分不等式的弹性接触问题基础框架。Wriggers[2]、[3]、[4]，Han[5]、[6]，Haslinger[7]以及Laursen[8]等人的现代研究探讨了求解这些非线性且可能不稳定的系统所需的计算方法。Chouly和Hild等人[9]、[10]总结了接触问题领域的最新进展。更多相关研究成果可见于[11]、[12]、[13]、[14]、[15]、[16]、[17]等近期文献。接触约束可以通过多种方法实现，包括拉格朗日乘子法[17]、[18]，罚函数法[15]、[19]、[20]、[21]、[22]，增广拉格朗日法[13]、[23]，以及砂浆公式法[24]。针对弹性接触问题，人们已经发展出了多种数值分析方法，例如优化后的Schwarz方法[25]，混合有限元法[FEM][19]、[26]、[27]、[28]、[29]，增广拉格朗日法[13]、[23]、[30]、[31]，Nitsche有限元法[4]、[32]、[33]、[34]，最小二乘法[35]，间断Petrov-Galerkin方法（DPG）[36]等等。

在接触力学中，所得到的方程组通常是一个大规模、非线性且不稳定的问题，尤其是在材料对比度较高的情况下。为应对这些挑战，我们提出了一种高效的迭代型接触求解混合框架，用于处理带罚函数的接触问题。该方法采用Robin边界条件作为传递条件，并结合一种无需导数的技术来更新界面上的传递数据（参见[37]、[38]）。在我们的模型中，接触边界完全位于较小的子域内，并与罚函数法结合使用，而较大的子域则允许灵活的离散化选择。例如，可以运用约束能量最小化的广义多尺度有限元法（CEM-GMsFEM）[39]、[40]，显著降低较大区域内高对比特征带来的计算成本，同时仍保持粗网格下的收敛性。更具体地说，我们在迭代型接触求解混合框架中引入了四种离散化类型：两种子域都采用标准/混合有限元法，或者较大子域采用标准/混合CEM-GMsFEM，而较小子域采用标准/混合有限元法。每种情况都配有相应的算法，为设计可扩展的求解器和预处理器提供了有效策略。

我们的主要贡献有三点。首先，我们开发了一种新颖的迭代型接触求解混合方法，用于多尺度接触力学问题，该方法将非线性接触约束限制在较小的子域内，而较大的子域则保持线性系统形式（见算法1、图2、图3、算法4）。传统方法通常会将整个区域的非线性接触问题一起处理，从而导致全局系统极为复杂。而在我们的方法中，非线性子问题被限制在一个宽度仅为一个粗网格单元的极小区域内，可通过半光滑牛顿法[41]高效求解。这种局部处理方式大大缩短了迭代时间。与此同时，较大的子域则支持灵活的离散化选择，无需考虑非线性问题。其次，我们提出了混合公式（见算法2和算法4），用于求解涉及高度异质性和高对比系数的线性弹性接触问题。通过采用应力-位移混合公式，我们的方法能够直接计算应力场，避免了后续处理步骤，同时确保了局部动量守恒，从而提升了物理一致性。此外，这种公式还能有效防止泊松比锁定现象，使得几乎不可压缩材料的模拟更加可靠。最后，我们将CEM-GMsFEM框架引入较大子域中（见算法3–4）。众所周知，在混合方法中实现强对称应力近似并非易事。虽然现有的混合方法[42]、[43]、[44]、[45]、[46]需要大量的自由度才能实现这种对称性，但我们的多尺度模型简化技术能够在保持精度的同时大幅降低计算成本。

我们需要强调的是，我们工作的核心创新之处在于提出了一种用于处理多尺度接触力学的混合方法——据我们所知，现有文献如[9]、[11]、[12]、[13]、[14]、[16]、[17]、[23]、[27]、[32]、[33]、[34]、[35]中尚未充分探讨这种方法。在这一混合框架下，非线性接触行为得到了有效处理。我们的方法与[37]、[38]、[47]中描述的方法有明显区别。需要注意的是，[38]研究的是应用于简单二阶线性椭圆问题的非重叠域分解方法，由于其线性特性，该方法的传递条件易于更新。而我们的工作重点则是高效处理非线性接触约束，沿界面更新传递条件是一项极具挑战性的任务。这一困难源于需要通过非线性迭代子程序来计算与接触区域相关的解，因此，我们算法中的核心数据更新步骤与接触约束的固有非线性紧密相关，使得该步骤变得极为复杂。

本文结构如下：第2节介绍模型问题及部分符号说明。弹性接触问题通过最小化定理和罚函数方法被重新表述。第3节介绍了迭代型接触求解混合框架以及各种离散化方案。具体而言，第3.1节和第3.2节分别详细描述了在整个区域应用标准有限元法和混合有限元法的方法，而第3.3节和第3.4节则讨论了在较大子域应用标准/混合CEM-GMsFEM、在较小子域应用标准/混合有限元法的方法。四种情况都配有相应的算法（即算法1、图2、图3、算法4）。第4节对第3节提出的迭代型接触求解混合方法的收敛性进行了分析。第5节通过两个测试模型进行了数值实验，以展示所提方法的性能。最后，第6节总结了本文的结论。

模型问题
我们考虑定义在Ω?Rn中的弹性接触问题：
Aσ￣=?￣(u) 在 Ω 内，
??·(σ￣(u))=f 在 Ω 内，
u=0 在 ΓD 上，
u·nc≤0，
nc·(σ￣·nc)≤0，
(nc·(σ￣·nc))u·nc=0 在 ΓC 上，
其中σ￣:Ω→Rn×n是对称应力张量，u表示位移，?￣(u)=12(?u+(?u)T)是线性化应变张量。nc是ΓC指向外的单位法向量。Ω?Rn（n=2,3）是一个有界且连通的Lipschitz多面体区域，其中填充着各向同性且线性弹性的固体。A是弹性算子的逆矩阵。

与各种离散化方式相关的迭代型接触求解混合方法
在本节中，我们介绍了迭代型接触求解混合框架的各种离散化方案。这些方案包括在整个区域应用标准/混合有限元法（见第3.1节和第3.2节），以及将较大区域中的标准/混合CEM-GMsFEM与较小区域中的标准/混合有限元法相结合的混合方案（见第3.3节和第3.4节）。

迭代型接触求解混合方法的分析
在本节中，我们对第3节提出的算法进行了收敛性分析。首先，我们在第4.1节证明了分解后的离散子问题与全局问题之间的等价性（定理1）。随后，第4.2节至第4.5节分别介绍了算法1、图2、图3、算法4的收敛性分析。

数值实验
在本节中，我们通过数值实验来展示所提出的迭代型接触求解混合方法的性能。实验涵盖了在两个子域都应用混合有限元法的情况，以及在较大子域应用混合多尺度方法、在较小子域应用混合有限元法的情况。所有计算均使用MATLAB 2021a在配备Intel Core i9-13900T处理器和32 GB内存的Lenovo ThinkCentre M80q Gen 4台式计算机上完成。

结论
在这项工作中，我们提出了一种高效的迭代型接触求解混合方法，专门用于处理具有高对比材料系数的多尺度接触力学问题。该方法将非线性接触约束限制在较小的子域内，而较大的子域则由线性系统描述。在该框架下，我们提出了四种不同的离散化策略，这些策略要么是在整个区域应用标准（混合）有限元法，要么是……