一种用于等几何配置方法的非重叠对偶-原始域分解求解器

《Computer Methods in Applied Mechanics and Engineering》：A non-overlapping dual-primal domain decomposition solver for isogeometric collocation methods

　　克里斯托斯·格克里扎利斯|马诺利斯·帕帕德拉卡基斯
希腊雅典15780，雅典国立技术大学佐格拉富校区，结构分析与抗震研究研究所

摘要
等几何配置法（IGA-C）作为一种替代等几何伽辽金法（IGA-G）的模拟方法被提出，因为它能在保持相同精度的同时降低计算成本。然而，为大规模IGA-C问题设计高效数值方案在算法和计算层面都比IGA-G模拟更为复杂。对于复杂几何区域，等几何分析需要将区域表示为多个子区域的并集，这为应用非重叠域分解方法提供了天然基础。尽管这类域分解方案已用于IGA-G模拟，但将其应用于IGA-C问题并非易事，因为不仅需要通过代数修改来处理非对称矩阵，还需对解决方案进行深入的数值评估，而目前针对此类问题的数学框架尚不完善。本研究提出了基于双-原始域分解方法（P-FETI-DP）的非重叠域分解方案，该方案经过适当修改以应对非对称矩阵，并作为预处理器加速广义最小残差法（GMRES）对相应代数方程组的求解。通过在典型测试问题上对所提方案的计算效率与数值可扩展性进行评估，验证了其有效性。

引言
等几何分析最初在[1]中被提出，作为基于多项式插值的有限元方法的替代方案。其核心思想是使用相同的基函数既表示精确的几何模型，又用于近似未知场。这样一来，等几何分析提供了一个统一框架，涵盖了物理域的几何设计及其离散化模型，从而省去了耗时的网格生成步骤，同时也使得在工程设计分析过程中应用自动自适应网格细化技术成为可能。虽然等几何分析最初基于广泛用于CAD领域的NURBS技术，但它也可以与其他样条技术结合作为基函数。由于NURBS的张量矩形拓扑结构不支持局部操作[2]，因此人们在[3]中提出了NURBS的泛化版本——T样条，它能够实现局部细化、面片的防水合并以及处理修剪曲面问题。而在[4]中，则研究了结合局部细化策略的B样条在自适应等几何分析中的应用。在多种工程应用中，等几何分析在精度方面优于传统有限元分析[5]，但其更高的精度也带来了相应代价：由于NURBS基函数具有更强的元素间连续性，导致在相同自由度数量下，基于伽辽金的等几何分析（IGA-G）的刚度矩阵比基于伽辽金的有限元分析（FEA-G）更为密集。这一特点不仅增加了刚度矩阵构建的成本，也提升了代数方程组求解的难度。因此，许多研究致力于降低等几何分析在刚度矩阵构建及方程组求解方面的计算成本。在[6,7]中，提出了适用于NURBS基函数的替代积分规则，以解决因元素边界处连续性过高而导致的计算效率问题，因为标准高斯积分规则在NURBS基等的几何分析中效率较低。此外，[8]还提出了利用并行计算和GPU编程环境，通过基于控制点交互而非元素交互的并行构建算法来加速系统矩阵的构建过程，进而提升性能。

等几何配置法（IGA-C）作为一种更快且能弥补IGA-G缺点的替代方案被开发出来。IGA-G基于边值问题的弱形式，而IGA-C则基于其强形式。在[9]中，有研究比较了IGA-G与IGA-C在求解标量及向量椭圆问题时的计算效率。德·洛伦齐斯等人[10]探讨了具有诺伊曼边界条件的IGA-C变体的精度与收敛性问题，以及其在接触力学问题中的应用。在[11]中，有研究分析了配置点的合理选择及其对模拟精度与收敛特性的影响。阿尼泰斯库等人[11]还提出了一种称为超配置法（IGA-SC）的IGA-C变体，该方法依据超收敛理论通过选择特定配置点来实现与IGA-G相近的收敛速度。类似地，戈麦斯和德·洛伦齐斯[12]以及莫塔尔迪尼等人[13]也研究了通过调整配置点位置来提升IGA-C计算效率的途径。虽然通过将单个面片域分解为相互重叠的子域，重叠施瓦茨预处理器已被用于IGA-G和IGA-C对应的线性方程组求解[14, [15], [16]，但非重叠域分解方案尚未在IGA-C框架中得到应用。在[17]中，有人提出了基于基本与增强型配置方案的两级混合DD型预处理器，其中GMRES求解法与施瓦茨型预处理器被调整为对应矩阵的舒尔补形式。

在几何上较为复杂的计算域通常由多个不同的NURBS面片组成，每个面片都有独立的几何参数映射。这一多面片参数化特性不仅为应用非重叠域分解方法提供了可能，还能生成质量更好的网格[1]。在[18,19]和[20]中，等几何非重叠撕裂与互连法（IETI）以及基于约束的域分解平衡法（BDDC）已在IGA-G模拟框架中被实现。不过，IGA-C的应用主要限于单面片域的情况[21]，而将非重叠域分解方法扩展到多面片IGA-C领域仍是一个未解决的问题。IGA-C公式的一个特征是会导致非对称刚度矩阵，且在混合边界条件下这些矩阵可能变为不定矩阵[17]，这一特点阻碍了非重叠域分解方法的直接应用。此外，与对称系统相比，针对非对称及不定系统的非重叠域分解方法的数学框架尚未得到充分研究[22]。在IGA-C框架中应用有限元撕裂与互连法（FETI）以及优化后的施瓦茨方法[23]，由于这些方法相比舒尔或施瓦茨等传统域分解算法需要更少的处理器间通信，同时仍能保持相同的并行水平[24]，有望进一步提升IGA-C处理实际大型工程问题的能力。

本研究旨在针对多面片IGA-C离散化问题，基于双-原始域分解方法（P-FETI-DP）开发一种非重叠域分解迭代求解方案，该方法可用于有限元模拟的求解[25,26]。P-FETI-DP是原始子结构化方法（PSM）与有限元撕裂与互连法版本FETI-DP的结合体[27,28]。通过对P-FETI-DP进行修改，以适应IGA-C对应的代数方程的特殊特性，研究发现该方法能够有效利用粗粒度和细粒度并行计算资源，从而显著提升IGA-C的计算效率。此外，本文还提出了一种针对多面片IGA-C问题的非重叠域分解方法，其思路与[18,29]中用于IGA-G公式的类似方法相似，即利用多面片模型中的各个面片来定义IETI及双-原始等几何撕裂与互连法（IETI-DP）子域。

本文结构如下：第2节简要概述了应用于单面片及多面片域线性弹性边值问题的等几何配置方案。第3节重点讨论IGA-C刚度矩阵的性质，以及用直接法与迭代法求解相应方程组的方法。第4节介绍域分解公式的整体框架，第5节则探讨重叠加性施瓦茨法（OAS）及其向多面片IGA-C的扩展可能性。第6节介绍了IGA-C离散化所需的舒尔补公式与P-FETI-DP方法。最后，在第7节中，通过典型基准问题进行数值测试，比较所提求解方案与现有方法的计算性能。

章节摘录
单面片域
本节简要介绍了应用于线性弹性边值问题的等几何配置方法。对于具有混合边界条件的情况，位移公式可表示为以下强向量形式：
{L(u)=?FinΩ?B(u)=tonΓ?Nu=gonΓ}
其中，L和B算子定义为：
L(u):=μ?2u+(λ+μ)?(?·u)?B(u):=2μEn+λ(?·u)n
这里，Ω表示受体力F作用的弹性体所在的开集，t表示给定的牵引力。

IGA矩阵方程
IGA-C刚度矩阵通常具有非对称性，且条件数往往很大，因此需要采用高效的数值方法来解决大规模问题。此外，相应产生的代数方程也是非对称的，且维度远高于相同精度下的IGA伽辽金法方程。混合边界条件的存在往往会导致矩阵变为不定矩阵[17]，进而影响迭代解法的收敛性。

域分解公式
设物理域Ω可由N个单独的面片Ωi组成，即：
Ω=∪i=1NΩi，且当i≠j时，Ωi∩Ωj=?
对于每个面片，可通过应用NURBS形状函数的几何映射F(i)，实现从参数空间到物理空间及反向的转换：
Ωi=F(i)(Ω^i)，?Ω^i=(F(i))?1(Ωi)
其中，Ω^i表示第i个面片的参数域。对于每个几何映射F(i)，可在相应子域上定义一个样条空间V^i。

重叠加性施瓦茨法的局限性
本节将简要介绍用于多面片IGA-C的加性施瓦茨预处理器。当域Ω由N个独立面片Ωi组成时，方程(14)的解u在Ωi范围内的限制可表示为ui=Riu，而u在整个域Ω范围内的扩展则表示为u=RiTui，其中Ri和RiT分别表示从全局自由度到局部自由度及反向的限制与扩展算子。由此，配置矩阵K的可限制形式为……

舒尔补公式的推导
根据关系Ωi′=Ωi∪Γi，对图3中与扩展子域矩阵相关的局部配置矩阵Ki、位移向量ui及力向量fi进行重新排列，可得：
Ki=(KIiIiKIiΓiKΓiIiKΓiΓi)，?ui=(uIiuΓi)，?fi=(fIifΓi)
其中，Ii和Γi分别表示子域i的内部自由度与边界自由度。同样，限制算子Ri的形式为：
Ri=(RIi00RΓi)
因此，对于矩阵K而言，有：
K=∑i=1NRiTKiRi
此处，边界自由度集合为……

数值测试
本节通过数值测试，评估所提出的迭代求解方案在解决用IGA-C方法离散化的多面片NURBS线性弹性问题时的收敛性能与计算效率。为此，将该迭代方案与Eigen C++模板库中提供的稀疏直接求解器SUPERLU进行了比较。改进后的P-FETI-DP方法相较于直接求解器所具有的真正计算优势将在后续测试中显现。

结论
对复杂结构的等几何分析通常需要对其几何形态进行多面片描述。目前，应用于等几何配置法（IGA-C）公式中的域分解方法仅限于单面片问题的重叠分解形式。此外，IGA-C离散化所得的刚度矩阵通常为非对称且往往不定，其元素缺乏传统有限元分析中那种清晰的物理解释。

数据提供方式
如需数据，可随时提出申请。