基于速度的空间-时间有限元方法在波动传播问题中的数值离散分析

《Journal of Computational Science》：Numerical dispersion analysis of velocity-based space–time finite element method for wave propagation problems

　　清水翔|维卡斯·夏尔马|藤泽和纪
京都大学农业研究生院，日本京都北区北白川追分町6068502

**摘要**
本研究构建了一个连贯的框架，将基于速度的时空有限元方法（v-ST/FEM）在波动传播问题中的理论色散分析与数值行为解释相结合。该分析分别评估了空间离散化和时间离散化对色散误差的贡献，并通过相对相位速度误差（RPSE）量化了它们的综合效应。分析结果通过RPSE等高线图进行了可视化，展示了不同空间和时间插值阶数组合对误差分布的影响：正的RPSE表明空间误差占主导，负的RPSE表明时间误差占主导。一个关键发现是v-ST/FEM主要受空间误差的影响，因此提高空间元素的阶数对于提高整体精度非常有效。值得注意的是，与传统的时间积分方法相比，结合高阶空间元素的v-ST/FEM在广泛的离散化参数范围内（包括元素大小Δx和时间步长Δt）具有更强的鲁棒性和精度。为了验证色散分析得出的结论，研究解决了一维和二维的波动传播问题。比较数值结果表明，RPSE等高线图与数值解中的振荡模式之间存在强烈的相关性：正的RPSE导致头部振荡，而负的RPSE导致尾部振荡。这些发现证实，使用高阶空间元素的v-ST/FEM提供了最有效的数值方案，能够实现几乎无振荡的波动传播问题解。

**引言**
波动传播现象对于理解结构（包括建筑物和地质结构）的瞬态响应至关重要。这些现象在材料的无损检测、声学和电磁散射以及土壤与结构在动态载荷条件下的相互作用（如地震事件中遇到的情况）中也起着重要作用。标量场u的波动方程的经典形式通过以下方程表示，并附有适当的边界和初始条件[1]、[2]：
?2u/?t2 ? c2?2u/?xi2 = 0
其中c表示波速。波动传播问题的解析解可以表示为平面谐波的叠加，即
u(x,t) = p1e^(ik1·x ? ω1t) + p2e^(ik2·x ? ω2t) + ……
其中i = ?1, kj, ωj, pj分别表示波矢、角频率和第j模式的系数。在非色散介质中，c仅由介质的固有属性决定，保持不变且与波数和角频率无关。这一特性表明所有模式都以相同的波速传播。因此，在非色散介质中，kj和ωj可以通过kj = ωj/c精确确定。

数值方法主要用于求解波动传播问题[3]。在这种情况下，基于网格的技术，特别是有限元方法（FEM）[4]、[5]被广泛使用。这些方法使用标准元素对运动方程进行空间离散化，得到一个常微分方程系统
MU?(t) + KU(t) = F(t)，
其中M和K分别表示全局质量和刚度矩阵[5]。同样，U、U?、U?和F分别表示对应于位移、速度、加速度和力的节点值向量。
方程(3)通常使用直接时间积分方法求解，其中有限差分方法（FDM）用于离散化时间域。这些时间积分方法可以分为显式和隐式方法[6]、[7]。通常，由于计算效率和实现简便性，显式时间积分方案更受青睐。然而，显式方法具有条件稳定性，需要将时间步长（Δt）保持在临界阈值以下，该阈值通常受网格中最小元素的影响[8]。相反，隐式方法可以实现无条件稳定性，确保算法的稳定性不受Δt选择的影响。因此，这些方法常用于计算中到长时间段的瞬态响应。

当应用有限元方法（FEM）解决波动传播问题时，方程(1)中提到的连续系统在时空域中被表示为离散系统（参见方程(5)）。因此，对于给定的波模式kj, ωj, pj，数值解uh(x,t)在波模式的表示上可能会不准确。这些差异可以通过数值波矢kjh、数值角频率ωjh和数值振幅pjh来量化，具体内容将在后续章节中详细说明。
uh(x,t) = p1e^(ik1h·x ? ω1ht) + p2e^(ik2h·x ? ω2ht) + ……
来自离散系统的数值解uh(x,t)在波的相位和振幅上表现出数值色散[3]。特别是，kjh和ωjh通常与kj和ωj不一致。这些差异的幅度显著受到kj, ωj的影响：对于高波数和频率模式，误差变得更加明显（此后称为高频区域）。这种现象是由于时空网格无法准确解析这些高模式所致。因此，每种模式的数值波速?cjh?ωjh/kjh与实际值c不同。在这种情况下，数值解表现出与色散介质相关的特性（见图1）。数值色散误差的存在在计算广阔时空域的解时带来了相当大的挑战。这是由于这些误差的累积性质，会损害精度。当波在介质中传播时，它们会污染数值解[9]。

许多研究调查了一维（1D）线元素[9]、[10]、[11]、二维（2D）矩形和三角形元素[12]、[13]、[14]、[15]以及三维（3D）域中使用六面体元素[16]、[17]的空间离散化相关的色散特性。这些发现表明，使用一致质量矩阵的标准Galerkin有限元方法（FEM）会导致数值波速ch≥c增加和波数kh≤k减小。相比之下，使用集中质量矩阵会产生不利影响。

Mullen和Belytschko[12]证明，在二维问题中，相位色散误差表现出各向异性，表明误差受波传播角度的影响，并随元素纵横比增加而增加。此外，观察到四边形网格比三角形网格具有更好的色散特性。Marfurt[13]利用一致质量矩阵和集中质量矩阵的广义平均值来减轻色散误差。随后，Christon[14]在Marfurt[13]的工作基础上，对双线性四边形进行了扩展，研究了各种质量矩阵对相位误差和群误差的影响。

众所周知，高阶空间插值显著减少了空间色散误差[9]、[11]、[18]、[19]、[20]、[21]。然而，这种方法引入了一个新的挑战，即在色散误差谱中出现通过带和停止带。在这方面，对于二阶元素，Belytschko和Mullen[11]在误差谱中识别出两个通过带和一个停止带。已经证明，停止带的存在会导致数值解中的虚假振荡。Thompson和Pinsky[18]使用Gauss–Legendre–Lobatto点对分层Legendre、分层Fourier和Lagrange多项式进行了详细研究，发现p阶空间插值表现出p个通过带和停止带，其中包括一个截止区域。他们推测p阶元素在小波数区域提供2p阶的精确近似，这一点后来被Ainsworth[22]证实。

在二维和三维中，Seriani和Priolo[23]引入了谱元素方法（SEM），该方法使用在高斯–Chebyshev–Lobatto点处评估的张量积Lagrange多项式。这种方法有效缓解了声波传播问题中常见的网格色散问题。此外，Mulder[24]证明使用Gauss–Legendre–Lobatto节点的SEM也可以获得更高的精度。

除了之前讨论的平均质量矩阵技术[14]和高阶有限元方法（FEM）[3]外，还出现了多种新的策略来减轻FEM中的色散误差[25]、[26]、[27]、[28]、[29]。例如，Harari和Hughes[30]、Thompson和Pinsky[31]、Oberai和Pinsky[32]、Oberai和Pinsky[33]以及Magoulès和Zhang[34]为Helmholtz方程、声波传播和时谐弹性波引入了一种稳定的FEM方法，显著减少了色散误差。此外，还开发了传统FE基函数的增强方法来改善FEM的色散行为[35]、[36]。对于时变的一维和二维波动传播问题，Kohno等人[37]、Ham和Bathe[37]、[38]提出了全局增强有限元方法，通过使用谐波函数（波包）来控制数值色散和耗散，并表明该方案可以通过增加增强函数的数量来控制色散和耗散。最近，Kumar和Kapuria[39]提出了一种新的增强FEM方法，该方法结合了单位分割[40]和波包增强的概念来定义局部元素域中的增强函数。这种方法已应用于各种类型的一维和二维波动传播问题[41]、[42]、[43]、[44]，证明了其在高频激励问题上的准确性和鲁棒性。

尽管空间离散化技术取得了显著进展，但完全消除数值色散误差仍然是一个持续的挑战[45]。这一挑战源于瞬态波动传播问题中空间和时间离散化的固有相互依赖性[5]。尽管如此，这种相互依赖性激发了通过研究各种时空离散化方法组合来研究色散误差的兴趣。Krieg和Key[10]对将有限元方法（FEM）与各种时间积分方法结合的效果进行了开创性研究。他们的发现表明，在一维波动传播问题中，使用FEM与集中质量结合中心差分时间积分是最优的。这是因为前者的相反色散特性（kh≥k）和后者的色散特性（ωh≥ω）可以有效地相互抵消。此外，他们提出一致质量（kh≤k）与梯形时间积分（ωh≤ω）的组合是另一种有利的替代方案。后续研究证实了这些发现[11]、[46]、[47]。类似类型的色散分析自然扩展到使用各种FEM和时间积分方案组合的多维问题[13]、[25]、[27]、[36]、[48]、[49]、[50]、[51]、[52]、[53]、[54]、[55]、[56]。

最小化总色散误差的一个关键因素是空间和时间精度之间的对齐。Donea等人[57]强调，高阶空间近似不应受到低阶时间近似的损害。然而，广泛使用的时间积分方法通常被归类为单步方法。因此，这些方法通常只能达到二阶精度[58]。这种差异往往导致次优结果，特别是在长期模拟的背景下[59]、[60]、[61]。

时间不连续的Galerkin方法作为一种有前途的高阶精确时间积分方法出现，可以无缝与高级高阶空间插值技术结合使用。这种方法的基础可以在[62]、[63]中找到，最初关注的是一阶双曲方程。随后，这种方法被扩展到包括二阶双曲方程[8]、[64]、[65]、[66]、[67]、[68]、[69]，从而形成了时空有限元方法（ST/FEM）。

最近，Sharma等人[70]引入了一种专门为弹性动力学设计的时空有限元方法（ST/FEM）的变体，称为基于速度的ST/FEM（v-ST/FEM）。v-ST/FEM的关键特征是它仅将速度场作为主要未知数。与原始ST/FEM相比，这种策略导致的线性系统相对较小[8]。类似策略在ST/FEM中也有应用[71]、[72]。对于线性时间元素，v-ST/FEM在时间上是无条件稳定的并且三阶精确的。然而，它不能像原始ST/FEM那样实现高频的渐近消除[73]。为了克服这一限制，引入了耗散控制的v-ST/FEM[74]，该方法集成了一种人工衰减机制来抑制高频模式。然而，这种策略将时间的精度从三阶降低到了二阶。

出于控制高频成分耗散的愿望，Sharma等人[75]引入了位移函数的概念，以开发广义基于速度的ST/FEM。他们通过使用适当的位移函数，展示了如何以统一的方式研究和实现原始的两场ST/FEM [8] 和经典v-ST/FEM [70]，其中唯一的未知数是速度场。此外，通过这两种方法的线性组合，可以衍生出另一种方法，该方法具有用户定义的参数α来控制高频的耗散。到目前为止，ST/FEM的色散行为通常仅从时间离散化的角度进行考察。然而，如前所述，评估时空离散化的综合效应对于最小化总色散误差至关重要。由于缺乏这种评估，ST/FEM的数值色散行为尚未被完全理解，因此用于解决波动传播问题的ST/FEM的发展指标也很少。本研究的主要目标是通过全面调查基于速度的时空有限元方法（v-ST/FEM）的数值色散特性，为使用ST/FEM解决波动传播问题提供指导方针。此外，该研究还旨在通过观察数值结果中的误差模式来识别色散误差的来源，从而将理论色散分析与实际解释联系起来。

本文的其余结构安排如下：第2节简要概述了传统的半离散FEM和广义v-ST/FEM。然后，在第3节中描述了色散分析的方法，并详细介绍了误差结果以及与广泛使用的时间积分方法的比较。第4节提供了验证理论误差分析的数值模拟。最后，第5节总结了研究结果并提出了结论性意见。

**时间积分方案和广义v-ST/FEM**
在深入探讨色散分析之前，首先提供了方程（1）的完全离散化版本的简要总结，这些总结利用了半离散有限元方法（FEM）和时空有限元方法（ST/FEM）。半离散形式方程（3）是这两种方案的基础。值得注意的是，用于推导半离散形式的空间离散化阶数不会影响后续的实施过程。

**广义v-ST/FEM中的连续时间时空误差**
本节研究了由于ST/FEM中涉及的空间-时间离散化而产生的相位和幅度色散误差。色散分析的过程遵循Kwon等人[54]和Chai与Zhang[55]为Bathe方法提出的方法。这里考虑了无界域中的一维和二维标量波动方程。由于两种维度中的过程非常相似，因此仅展示二维情况。

**数值结果**
本节的目的是验证第3节中关于所提方法相关的时空色散误差的结论。为此，使用基于第2节所述原理开发的开源软件研究了两个数值问题。该软件使用EASIFEM [96]框架开发，这是一个用现代Fortran实现的开源有限元方法（FEM）框架，由作者创建。

**结论性意见**
为了建立使用时空有限元方法解决波动传播问题的指导方针，本研究调查了基于速度的时空有限元方法（v-ST/FEM）的数值色散特性，并提供了从理论色散分析和数值结果解释中得出的实际建议。时空离散化的标量波动方程的色散分析是通过分离空间和时间部分来进行的。

**CRediT作者贡献声明**
Shimizu Shion：撰写——原始草稿、可视化、验证、软件、方法论、资金获取、概念化。
Vikas Sharma：撰写——审阅与编辑、监督、软件、方法论、资金获取。
Fujisawa Kazunori：撰写——审阅与编辑、监督、项目管理、资金获取。

**利益冲突声明**
作者声明他们没有已知的可能会影响本文所述工作的竞争性财务利益或个人关系。

**致谢**
第一作者获得了JSPS KAKENHI研究员资助（DC2）（资助编号JP25KJ1575），第二作者获得了JSPS KAKENHI早期职业科学家资助（资助编号JP24K17983）。该研究还得到了JSPS KAKENHI科学研究资助（A）（资助编号JP25H00951）的支持。