非线性弹性动力学一阶双曲型表述的隐式双时间步进：收敛区域与预处理

《Computers & Structures》：Implicit dual time-stepping for first-order hyperbolic formulation of nonlinear elastodynamics: Convergence regimes and preconditioning

　　

# 非线性弹性动力学一阶双曲{p,F}格式中DTS内部求解器的系统评估

## 摘要

非线性弹性动力学的一阶双曲{p,F}格式，其中p为线动量、F为变形梯度，使得源自计算流体力学（CFD）的有限体积方法可以应用于固体动力学。双时间步进（DTS），作为CFD中保证时间精度的标准方法，被扩展至固体求解器，为流固耦合（FSI）提供了统一的基础架构。本文所采用的近不可压缩neo-Hookean（NINH）模型使得收敛性具有挑战性，因为其特征结构随变形而变化，且其压力波速在近不可压缩极限下趋于无穷。然而，针对该系统DTS内部求解器的系统性评估尚不存在。本工作首次在一阶BDF2 DTS格式框架下表征了三种内部求解器：显式多阶Runge–Kutta（MRK）、LU对称Gauss–Seidel（LU-SGS），以及带有LU-SGS预处理的FGMRES。对时间步长、CFL数、Krylov参数、泊松比及网格拓扑的系统扫描表明，仅有带预处理的FGMRES在整个范围内收敛。增强后的{p,F,J}系统，将体积映射J作为守恒变量加入，恢复了{p,F}系统失效时的收敛性。我们早前的角动量校正方案被重新表述为BDF2形式，并能在独立于内部求解器的情况下将动量守恒至机器精度。本文结果为DTS的内部求解器选择提供了指导，并将时间谱方法（TSM）的扩展作为未来工作。

## 引言

计算固体动力学（CSD）传统上采用基于动量方程二阶位移格式的有限元方法（FEM）进行建模。这些方法已经成熟并被广泛应用，但在接近不可压缩极限时，低阶位移单元可能出现体积锁定现象。锁定的产生是因为不可压缩约束在每个积分点处施加后，留给表征变形的自由度太少，从而导致响应异常刚硬。尽管已有补救方案被提出[1][2][3][4]，但固体动力学的有限体积格式对于流固耦合（FSI）而言是更自然的选择，因为流体求解器本身已经基于有限体积方法（FVM）。将FEM结构求解器与FVM流体求解器进行耦合已是成熟的且有效的方法，数十年的任意Lagrangian–Eulerian流体格式配合专门界面处理的气弹性工程实践已充分证明了这一点[5][6]。由于大多数生产级CFD求解器都是有限体积方法，用相同的离散方法处理固体，可以使得两个场共享一套数据结构和一套求解器框架。它还允许使用共形界面，在共享节点上直接交换牵引力和速度，而不是在离散不匹配的界面之间进行传递[5]。将固体动力学引入这一框架的方法由一系列先驱性工作提出[7][8][9][10][11]，将非线性弹性动力学方程重新表述为在线动量和变形梯度上的一阶双曲守恒律系统，即{p,F}系统。该格式类似于可压缩气体动力学的Euler方程，可以直接使用为CFD开发的有限体积技术，包括近似Riemann求解器、稳定化方案和时间积分方法。由于变形梯度作为主要守恒变量携带，而非通过对位移场求导恢复，该格式达到了与几何相同的收敛阶数（针对速度和应力），并绕过了体积锁定的运动学根源[8][11][12]。近不可压缩性仍然通过压力波速影响系统。随着泊松比ν→1/2，体积模量κ无限增大，压力波速也随之增长。对显式时间积分而言，这表现为越来越严格的Courant–Friedrichs–Lewy（CFL）稳定极限，这也是本工作中开发隐式双时间步进框架的主要动机之一。

{p,F}系统的有限体积处理是近年来研究的课题。Lee等人[8]开发了基于Riemann通量计算的单元中心迎风求解器，Aguirre等人[9]引入了采用Jameson–Schmidt–Turkel（JST）格式[13]进行人工耗散的顶点中心格式。后续扩展涉及了线性重构的迎风方案[10]、稳定化Petrov–Galerkin方法[11]以及任意Lagrangian–Eulerian守恒格式[14]。另有独立的研究路线使用稳定化有限元处理相同问题[15]。这些研究已建立了守恒律方法处理大变形、波传播和几何非线性等固体动力学问题的可行性和准确性。

该格式相较于位移基有限元的优势基于开发它的工作中确立的性质。除上述等阶收敛和缺乏运动学锁定机制外，顶点中心和稳定化变体被构造为抑制在不可压缩极限附近污染等阶混合位移-压力单元的虚假压力模式[7][9][11]。这些是空间离散化的属性，而非本文所开发的隐式时间积分方法的属性。有限体积与有限元求解器之间的跨代码计时测试衡量的更多是实现方式而非方法本身，因此第6节中的运行时比较是在同一框架内比较内部求解器，而非两种方法之间的比较。

尽管与Euler方程存在结构类比，但{p,F}系统在影响隐式时间积分的方面有所不同。此处的比较是计算层面而非物理层面的，因为两个系统描述不同的连续介质。其守恒变量包含完整的变形梯度张量，因此双曲系统每节点携带nd+nd2个未知量，而非Euler方程的nd+2个，这扩大了特征系统和隐式方法必须处理的隐式算子。理想气体Euler特征结构仅依赖于局部热力学状态，而{p,F}系统特征结构通过弹性张量依赖于本构模型，并随局部变形而变化。在接近不可压缩极限时，ν→1/2，体积模量趋于无穷，压力波速无限增大，而剪切波速保持有界。这些性质使得隐式内部求解器的收敛性和鲁棒性在可压缩流中没有直接对应物，因此本文按自身标准对其进行了表征。

{p,F}族处于更为广泛的固体力学有限体积方法体系中，Cardiff和Demird?i?对此进行了综述[16]，其中基于动量方程位移形式的隐式和分裂求解程序是常见做法，双时间步进配合GMRES加速已被应用于完全耦合的FVM基FSI[17]。然而在{p,F}一阶双曲族内部，此前的工作几乎完全依赖显式时间积分，通常是二阶全变差缩减（TVD）Runge–Kutta（RK）方法。显式方案易于实现且非常适合快速瞬态问题，但它们受制于严格的CFL稳定条件，该条件基于网格尺度和弹性模型来约束时间步长。这一限制在精细分辨网格、接近不可压缩极限时体积模量很高的高刚度材料参数、需要长时间物理积分或求解稳态（显式步数总数可能因严格条件稳定而产生极大计算时间）的问题中变得越来越严重[18]。

在CFD中，双曲守恒律经过数十年深入研究，这一挑战推动了广泛隐式时间积分方案的发展。Yoon和Jameson[19][20][21]引入的LU对称Gauss–Seidel（LU-SGS）方案通过前向和后向三角扫描提供高效的近似Newton迭代，避免了完整Jacobian矩阵的构造和求逆。全矩阵分裂变体[22][23]和块LU公式[24]以额外的矩阵操作为代价提高了收敛速率。对于需要快速收敛至严格容差的问题，Krylov子空间方法，特别是柔性广义极小残差（FGMRES）算法[25]，已被开发出来，通常与利用上述近似分解方案的基于物理的预处理器结合使用[22][26][27]。对于瞬态问题，这些隐式求解器被嵌入双时间步进（DTS）框架[18][28][29]中，在该框架中，伪时间迭代驱动解趋向下一物理时间层级的隐式系统。DTS是非稳态CFD中的标准框架，因为网格在边界层和激波处高度集中，显式全局时间步长将由域中最小网格单元决定，因此是不切实际的，而伪时间中的局部时间步进通过允许每个单元以自身的伪CFL推进来恢复效率。DTS格式还将物理时间精度与伪时间收敛过程解耦，使得局部时间步进和残差平滑等收敛加速技术可以被自由使用而不影响物理时间积分。

除独立CFD之外，DTS框架也是单体式流固耦合（FSI）框架的自然基础，在该框架中流体和固体在相同的有限体积环境中处理并以共同物理时间步长推进。基于DTS的固体求解器直接适用于该场景，因为它接受共同物理时间步长，不论固体自身的稳定极限如何，其伪时间迭代在每个FSI时间步长收敛隐式固体系统。相同的伪时间结构还支持降阶时间谱方法（TSM），如用于周期性流动的谐波平衡（HB）方法（也称为傅里叶基时间谱方法，F-TSM）[30][31]和用于周期性和非周期性问题的Chebyshev基时间谱方法（C-TSM）[32]，在这些方法中，非稳态问题被重新表述为在多个时间层级耦合的数学稳态系统。因此，表征{p,F}系统的伪稳态内部求解器的收敛行为和鲁棒性是DTS时间推进、FSI耦合以及TSM向CSD扩展的共同前提。总Lagrangian（TL）{p,F}格式中计算网格固定在未变形构型上，兼容需要在固定网格上于多个时间层级进行评估的TSM，这是增量Lagrangian（UL）或任意Lagrangian–Eulerian（ALE）格式所不具备的优势。

尽管这些时间积分技术在CFD中已经成熟，但它们在CSD中{p,F}系统的应用仍然很大程度上未被探索。关于有限体积弹性动力学的现有文献尚未评估在双时间步进框架内不同内部伪稳态求解器的收敛速率、鲁棒性和工作范围。这种表征是DTS基流体求解器FSI耦合的前提，是将TSM降阶时间方法扩展到CSD的{p,F}格式的前提，也是针对DTS效率已被证明具有决定性作用的内部求解器加速技术（如多重网格）进行定向未来开发的前提。

最近在这方面的贡献是Cardiff等人[33]的工作，他们基于动量方程经典二阶位移形式开发了用于单元中心有限体积固体力学的无Jacobian Newton–Krylov（JFNK）程序，并在涵盖线性、超弹性和超弹塑性材料的一组静态和动态问题上，将其与采用GMRES配合代数预处理器的常规分离拟Newton算法进行了基准测试。

本工作表征了非线性弹性动力学一阶{p,F}双曲守恒律格式双时间步进处理中内部伪稳态求解器的收敛区域、鲁棒性边界和实际操作指导。隐式求解器以无矩阵方式实现，在采用单元-顶点对偶控制体的网格透明、基于边的有限体积代码中。中位对偶构造可扩展至一般凸多面体，本文验证的混合单元网格由三角形、四边形、四面体、六面体、棱柱和金字塔组成，对偶网格在这些网格上构建。在此环境中，显式多阶Runge–Kutta（MRK）格式、隐式LU-SGS方法以及带和不带LU-SGS预处理的柔性广义极小残差（FGMRES）算法首次针对{p,F}系统进行了报告。LU-SGS方案既作为独立求解器使用，也作为FGMRES的预处理器，与此前研究采用的代数预处理器形成对比。通过对物理时间步长、CFL数、Krylov子空间参数、泊孙比、网格拓扑和CFL启动斜坡的系统扫描来表征这些方法，映射每个求解器族在不同区域的工作范围，包括近不可压缩和高刚度材料区域。在DTS设置中，物理时间步长和CFL数是独立参数，前者设定所解析的物理瞬态，后者设定内部迭代的伪时间步长。

增广{p,F,J}系统[7]、[34]，其中体积映射J作为额外的守恒变量携带，以实现对体积模态的显式控制，在近不可压缩性研究中与标准的{p,F}系统进行了对比。该特征化为基于DTS的流体求解器的FSI耦合以及将TSM扩展到{p,F}系统奠定了基础。我们在近期研究[35]中开发的广义角动量守恒算法也被改写为与BDF2一致的形式，并在每个物理时间步内应用一次，这使得校正与内层求解器无关，并且与DTS循环中的任意时间积分方案兼容。

各节摘要

守恒定律

在总体拉格朗日框架下，控制连续体的守恒定律为[36]、[37]：
?p/?t = ???P + ρ?f_b
?F/?t = ??v = ???[v?I]
?u/?t = v

其中，p = ρ?v 表示线动量，v = ?x/?t|_X 为物质粒子X的速度，u = x ? X 为从未变形构型 X = x(0) 到当前构型 x(t) 的位移，ρ? 为参考密度。第一Piola-Kirchhoff应力张量 P = P(F) 依赖于变形梯度 F = ??x = ?x/?X = I + ??u，而 f_b 为体力。

有限体积离散化

线法将式(7)中积分守恒定律的空间和时间处理分离，在每个节点i处产生一组耦合的常微分方程（ODE）：
ddt(q_V)_i = ?r_i(q)

残差r的计算在第3.1节中讨论，用于积分该ODE系统的时间离散化方案在第3.2节中介绍。

角动量守恒

连续体不需要单独的角动量守恒定律。Cauchy应力的对称性使力矩平衡成为线动量平衡的推论，而超弹性提供了该对称性。空间离散化破坏了这一平衡。控制体积面上的数值通量，包括人工耗散在内，不具有对称性，因此离散的力矩平衡不成立，尽管每个节点的应力仍然是对称的。在涉及自由旋转或……

算例

基准问题包括低色散立方体、旋转板、弯曲柱、铍板、L形块和扭转块。以上六例均为非线性固体动力学文献中的标准案例[7]、[9]、[11]、[15]、[59]，且均在上述参考文献中建立的构型上进行计算。其中，低色散立方体具有解析解，在5.1节中用于对精度阶进行定量检验。

系统求解器性能

求解器性能评估

为了隔离求解器参数对收敛性和计算成本的影响，本评估逐项进行。在接下来的内容中，将依次考察物理时间步长、CFL数、Krylov子空间参数、近不可压缩性、网格拓扑和CFL启动斜坡，并使用第5节中最佳示例来展示每一项。

本节参数扫描中报告的运行时间在Apple M1处理器单核上测量。

结论

本文首次对非线性弹性动力学一阶双曲型{p,F}格式的隐式双时间推进（DTS）的内层伪稳态求解器进行了特征化，采用基于网格透明、边基有限体积的求解器，在单元-顶点对偶控制体积上实现。据我们所知，这是针对该系统隐式DTS内层求解器的首次系统研究。推导了{p,F}和{p,F,J}系统的控制方程、特征结构和隐式算子。内层求解器……

CRediT作者贡献声明

Burak Cona：写作—审阅与编辑，写作—原稿撰写，可视化，验证，方法学。
Kivanc Ekici：写作—审阅与编辑，写作—原稿撰写，监督指导。

利益冲突声明

作者声明不存在已知的可能影响本文工作的财务利益冲突或个人关系。

致谢

作者感谢田纳西大学诺克斯维尔分校提供的支持。

Burak Cona | Kivanc Ekici