一种基于域分解方法的解弹塑性接触问题的解耦算法及其在拱坝中的应用

《Advances in Engineering Software》：A decoupling algorithm based on domain-decomposition method for solving elastoplastic contact problems and its application in arch dam

　　•为解决弹塑性接触问题，人们开发出一种名为DD-BDNM的算法。•该算法采用了双层迭代策略。•只需重新计算存在塑性应变的子域的刚度矩阵即可。•无需反复生成、组合和调整全局刚度矩阵。•DD-BDNM可应用于拱坝的静力与动力分析。

引言
摩擦接触问题是结构分析中最常见的非线性问题之一[[1], [2], [3], [4]]。这类问题涉及两个或多个相互接触的物体之间的相互作用，需满足的物理条件如下：（1）法向接触条件：接触区域内的物体不得互相穿透；法向接触力为非负值，即只能承受压应力。（2）切向摩擦条件，即库仑摩擦定律：切向接触力的大小不超过摩擦系数与法向接触力乘积的值，且其方向与相对滑动方向相反。

弹塑性问题是与弹性问题相对应的一个概念。弹性问题描述的是材料的理想状态：当结构受到外力作用时，其应力与应变的增加始终保持线性关系，且在载荷移除后结构能够完全恢复到初始状态。而材料的弹塑性行为则表现出两种不同的应力-应变关系。具体而言，当结构应力低于弹性极限时，材料表现为弹性应力-应变行为；当应力超过弹性屈服点后，应力-应变关系的斜率会随着外力的增加而发生变化，后续行为可能以线性或非线性方式呈现[[5], [6]]。

在工程分析中，尤其是对于那些具有重大社会和经济价值的项目，如拱坝，同时考虑接触非线性和材料非线性对于提高数值结果的准确性与可靠性至关重要。然而，解决弹塑性接触问题的主要难点在于材料非线性与接触非线性之间的相互耦合。这些耦合项包括刚度矩阵和接触力：前者由弹塑性本构模型决定，决定了系统的接触状态；后者则是通过求解接触方程得到的，会影响残余力的计算。这两者相互影响，共同决定了弹塑性接触问题控制方程的解，增加了数值计算的复杂性以及算法收敛的难度[7,8]。此外，当需要求解的方程规模较大或接触节点对数量较多时，解决弹塑性接触问题不仅面临上述挑战，还会出现接触状态的变化以及塑性应变的产生，从而导致需要反复生成、组合和调整全局刚度矩阵，进一步增加计算机存储和方程求解的难度。

BDNM是一种用于解决摩擦接触问题的有效算法，它由Pang在20世纪90年代首次提出，是将牛顿法引入B可微分方程方法中而形成的[9]。该方法能为非线性规划、变分不等式以及非线性互补问题建立统一的求解框架，并且对于非光滑方程具有全局收敛性。此后，Harker等人[10]又为该方法引入了一维线性搜索技术，证明了BDNM的计算效率与稳定性。在此基础上，Qi等人[11]从理论上证明了可以使用广义雅可比矩阵来替代对控制方程关于独立变量导数的计算，从而极大推动了BDNM在数值计算领域的应用[[12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27]]。

鉴于接触条件的非光滑特性以及BDNM在处理此类问题方面的显著优势，Christian等人[[28], [29]]率先将这种方法应用于摩擦接触问题的研究。在他们的工作中，平衡方程与接触条件被同时表示为B可微分方程的形式，随后推导出相应的雅可比矩阵。之后，通过使用BDNM求解控制方程，得到了接触力与节点位移的值。基于此项工作，Christian进一步将这种方法扩展到二维弹塑性摩擦接触问题的研究中[30,31]。Str?mberg[32]则将BDNM应用于三维磨损问题的求解，并通过数值示例验证了BDNM在收敛速度与稳定性方面的优异表现。González等人[33]则采用BDNM与边界元方法共同研究二维弹性滚动接触问题，通过与参考解的对比，分析了离散化方法及积分步长对计算结果的影响。

通过回顾BDNM的发展及其在接触问题中的应用情况可以发现，现有的研究大多集中在弹性摩擦接触问题的求解上。根据已发表的文献，Christian等人[30,31]已将BDNM扩展到二维弹塑性摩擦接触问题的研究之中。不过对于二维问题而言，接触体的滑动方向被限制在一条直线上，这并不符合三维接触问题中接触体在切平面内滑动方向可为任意方向的实际情况，而这实际上也是导致三维摩擦接触问题难以求解的主要原因之一。此外，Christian等人[30,31]采用的是交替求解接触方程与弹塑性平衡方程的策略。设问题的总自由度为N，接触节点对的数量为NC。对于三维弹塑性接触问题，每一步迭代中需要求解的方程维度为N + 3 × NC。当待解决问题规模较大或接触节点对数量众多时，他们所提出的算法将会耗费大量时间。另外，Yin等人[34]也将BDNM应用于弹塑性接触问题的求解，其中通过带有双拉格朗日乘数的砂浆段间接触模型来施加接触约束，以应对接触边界处非协调网格的问题，同时还采用了变量替换策略来修改刚度矩阵。

在有限元分析中，域分解方法是处理大规模问题的常用手段，其中Farhat等人在20世纪90年代提出的有限元撕裂与互连（FETI）算法便是较为著名的方法之一。在FETI算法中，结构被分解为若干个不重叠的子域，那些具有足够多狄利克雷边界条件的子域被定义为固定子域，而那些存在刚体运动的子域则被归类为浮动子域。随后，通过引入对应于物理中界面力的拉格朗日乘数，来强制满足相邻子域人工界面上的位移连续性条件。同时，还会建立正交约束条件，以消除浮动子域的刚体运动[38]。FETI算法的实施包含四个关键步骤[39]，包括计算每个子域的系数矩阵；通过结合所有子域的系数矩阵来建立控制方程；求解控制方程以得到界面力与刚体运动；再将界面力与刚体运动代回相应的子域中，进而计算位移场。在这些步骤中，计算每个子域的系数矩阵与反向代入求解位移场是相互独立的，这一特点符合分布式内存架构的要求，因此可以采用并行计算[40]。此外，由于不再需要构建全局刚度矩阵，而是使用多个尺寸较小的刚度矩阵作为替代，FETI算法能够有效降低计算机内存的使用量以及需要求解的方程规模，从而显著提升解决大规模有限元问题的计算效率[[41], [42]]。

随着相关研究的深入，以Farhat和Dostál为首的学者们相继提出了基于FETI算法的一系列改进算法，包括双原始FETI（DP-FETI）算法[[43], [44], [45], [46]]、全FETI（T-FETI）算法[[47], [48], [49]]、混合FETI算法[[50], [51], [52], [53]]、自适应FETI算法[[54], [55]]等等。同时，通过严谨的数学证明以及大量的数值示例，这些算法的性能，尤其是在处理大规模工程问题方面的能力也得到了充分验证[[56], [57]]。与传统的单域有限元方法相比，FETI算法能够将大规模问题分解为多个局部子问题，通过并行计算获取各个子域的系数矩阵与位移场，从而实现大规模非线性问题的高效求解。例如，Cermak等人[58]基于T-FETI算法提出了一种针对动硬化率无关的弹塑性问题的并行求解算法。在FETI方法中，仅通过拉格朗日乘数来强制满足界面连续性条件，而狄利克雷边界条件则直接应用于固定子域。相比之下，在T-FETI方法中，所有的边界条件——包括狄利克雷条件——都通过拉格朗日乘数来施加，这使得所有子域都可以被统一视为浮动子域。据Cermak等人所述，首先以位移为未知量，以增量形式推导出控制方程，然后提出一种基于弹塑性应力-应变关系中切向算子线性化的外部算法，以适配T-FETI算法的并行计算特性。最后，借助FLLOP计算库，通过多个数值示例验证了所提算法的弱可扩展性。Ajith等人[44]则将DP-FETI算法与随机配置法相结合，提出了一种适用于弹塑性随机力学问题的可扩展域分解方法，有效解决了传统方法中存在的计算成本过高的问题。通过同步降低空间维度与随机维度，并充分利用并行计算的优势，大幅提升了大规模非线性问题的求解效率。在数值示例中，通过P-拉普拉斯方程与平面应变问题，展示了所提算法在计算精度与效率方面的双重优势，以及出色的可扩展性。

基于上述背景，有人在[59]中提出了一种名为FETI-BDNM的算法，用于解决大规模弹性接触问题。在作者看来，FETI-BDNM的优势包括：（1）充分利用了BDNM的全局收敛性；（2）将原有的全矩阵形式的接触柔度矩阵转换为块对角结构，显著提高了计算效率；（3）能够独立计算每个子域的系数矩阵与位移场，保持了FETI算法的并行可扩展性。本文进一步将FETI-BDNM扩展到弹塑性接触问题的求解中，提出了一种名为域分解B可微分牛顿法的算法。该方法通过将结构分解为子域，并采用双层嵌套迭代策略，实现了接触非线性与材料非线性的解耦。同时，利用并行计算的优势，高效地完成了每个子域的系数矩阵、位移场以及应力场的计算，从而快速求解弹塑性接触问题。本文的结构安排如下：第2节基于虚功原理建立了弹塑性接触问题的控制方程；第3节介绍了用于求解这些控制方程的双层嵌套迭代算法；第4节进一步推导了动态情况下的控制方程及相应的求解流程，第5节则通过数值示例展示了所提方法的性能。最后，本文在第6节给出结论。

**域分解弹塑性接触问题的控制方程建立**
在有限元方法中，单元内任意点的位移uˉ可通过节点位移插值表示为：
uˉ=Nu
其中N为形函数，u为节点位移向量。根据几何关系，应变ε与节点位移u之间的关系可表示为：
ε=LTuˉ
此处L为微分算子矩阵。将式(1)代入式(2)后可得：
ε=Bu
其中B为……

**基于域分解的B可微牛顿法求解控制方程**
式(14)体现了接触非线性和材料非线性的耦合，这要求需计算内力、界面力、接触力以及刚体运动，使得问题求解更为复杂。为解决这一问题，本文采用双层嵌套迭代策略来处理弹塑性接触问题。在该算法中，控制方程被视为弹性接触方程与弹塑性方程的结合体。

**域分解弹塑性动力接触问题的控制方程**
使用DD-BDNM求解弹塑性动力接触问题时，首先将结构分解为若干个不重叠的子域。在时间t和t+Δt时刻，子域s的平衡方程分别可表示为：
Fi,ts+Csu˙ts+Msu¨ts=Fe,ts+Fc,ts+Fj,ts
Fi,t+Δts+Csu˙t+Δts+Msu¨t+Δts=Fe,t+Δts+Fc,t+Δts+Fj,t+Δts
其中Fis、Fes、Fcs和Fjs分别代表内力、外力、接触力和界面力；Ms为质量矩阵；Cs为……

**静态接触问题**
图3展示了一个具有接触面的悬臂结构。该模型的尺寸为L × H × D = 20 × 6 × 6米。接触面将该结构分为两个相同的部分，每个部分 further被划分为4个尺寸为L? × D? × H? = 10 × 3 × 3米的子域，其中包括4个固定子域和4个自由子域。该模型通过有限元离散化处理，共包含1169个节点、720个单元以及140对接触节点对。集中载荷被施加到……

**结论**
本文提出了一种基于域分解的算法——DD-BDNM，该算法采用双层嵌套迭代策略与FETI-BDNM方法，用于求解弹塑性接触问题。所提算法具有以下四项显著优势：
(1) 弹塑性接触问题的求解被分解为多个维度较低的局部子问题，从而有效控制了需要求解的方程规模，提升了计算灵活性。
(2) 在……

**CRediT作者贡献声明**
高航铎：原文撰写、软件应用、方法设计、实验分析。何先锋：指导工作。兰艳：结果验证、资源提供。杨浩明：方法设计、形式化分析、数据整理。徐欣：资金筹集。林高：指导工作。
高航铎|何先锋|兰艳|杨浩明|徐欣|林高