物质点方法的稳定性分析及不稳定性缓解措施

《Computer Methods in Applied Mechanics and Engineering》：Stability analyses and instability mitigation for the material point method

本研究通过能量分析和谱特性分析，探讨材料点法（MPM）单步更新的数值稳定性，推导出临界时间步长公式及其影响因素，并提出一种基于改进形状函数梯度的稳定化策略，有效扩展了稳定时间步范围，为实际应用提供指导。

　　

杨文嘉|黛博拉·L·苏尔斯基

国立中兴大学土木工程系，145兴达路，台中市，40227，台湾

摘要

本研究通过分析单步更新的能量行为和谱特性，探讨了材料点法（MPM）的数值稳定性。首先进行基于能量的分析，以量化每个时间步长引入的能量变化，然后通过放大矩阵进行谱分析，识别出关键的时间步长约束。推导出了临界积分参数和无量纲临界时间步长的封闭形式表达式，强调了它们对质量参数和粒子分布的依赖性。稳定性分析确定了影响MPM稳定性极限的关键因素，包括质量矩阵的选择、速度投影、部分占据的网格单元以及基于粒子的公式中的积分误差。数值实验验证了分析预测，并揭示了基于粒子的积分误差对稳定性的影响。提出了一种简单的稳定化系数，该系数修改了部分填充边缘单元中的形状函数梯度，显著扩展了稳定时间步长范围，而不会增加计算成本。所提出的框架为选择稳定时间步长和提高MPM模拟的鲁棒性提供了实用指导。

引言

材料点法（MPM）[1] [2] 是从流体力学扩展到固体力学的FLIP [3]，广泛用于模拟涉及大位移、强非线性、复杂接触和多相相互作用的问题。它用离散的材料点（粒子）表示连续体，并使用背景计算网格来更新它们的状态变量。大多数实现使用中心差分方案显式推进网格，因为该算法简单且高度可并行化。为了保持稳定性，时间步长通常是根据经验选择的，通常是具有相同显式积分器的有限元模型中使用的临界时间步长的小部分（例如，十分之一）。这些保守的选择大大增加了计算成本，并限制了大规模或实时模拟的可行性。本研究确定了真正控制稳定性的参数，推导出了临界时间步长的封闭形式表达式，通过数值实验验证了这些选择，并引入了一种实用的稳定化策略，可以在几乎不增加额外成本的情况下放宽对步长限制的要求。

大量研究已经解决了单元交叉误差问题，这是由粒子跨越单元边界时引起的积分误差问题。在标准MPM中，多线性网格形状函数产生的应变增量场只是分段连续的。单元接口处的不连续性会产生虚假的脉冲和能量尖峰。广义插值材料点（GIMP）方法[4]和传导粒子域插值（CPDI）方法[5] [6]通过将每个粒子扩展到有限域并在粒子级别减少积分误差来缓解这一问题。加权最小二乘和相关应力/传递场的重构[7] [8] [9] [10]提供了补充的解决方案。尽管这些方法显著提高了精度和收敛性，但它们并没有增加显式方案的最大稳定时间步长，这促使人们研究完全隐式的时间积分方法[11] [12] [13]。

单元交叉解决方案还有助于抑制由重复的粒子-网格映射引起的不稳定性。因为粒子可以携带比网格所能表示的更高空间频率的信息，信息可能会以两种方式被错误表示：（i）更高频率的分量可能在节点处不可见，不会贡献节点信号；（ii）其他分量可能会混叠成较低频率，并与网格上的其他低频分量无法区分。这两种类型的错误表示分别激发了虚假的零能量模式（例如沙漏和棋盘图案）并引起人工共振，通常称为“ringing”[14]。通过使用更高阶的形状函数[15] [16] [17]或显式滤波器[18] [19] [20]增加网格的表示能力，可以减少这种映射引起的不稳定性。由于这些错误表示只是扰动了而不是根本改变了基于网格的时间积分器所规定的临界时间步长，因此在解决单元交叉问题后观察到的严格步长限制主要应归因于积分器。

临界时间步长的稳定性条件可以从放大矩阵的谱分析中得出，该放大矩阵控制着网格节点值的更新。Berzins [21]将Spigler–Vianello非线性稳定性框架[22]应用于MPM/GIMP（应力最后），得到了一个包含变形、应力和粒子支持的显式步长公式。这个限制对于一维固定-固定杆是有效的，而有限元分析中常用的简单声速时间步长限制则不够充分。Ni和Zhang [23]为各种应力更新方案下的一维系统推导出了精确的临界时间步长限制。使用移动网格（完全拉格朗日）公式，他们对一个实际的两单元模型进行了谱（特征值）分析，该模型具有完全重构的节点质量。通过假设正交网格并对坐标方向取最小值，他们的限制可以扩展到更高维度。对于更新后的应力最后（USL）方案，两单元公式在移动网格、类似有限元的设置中基本上与临界步长相匹配，但在标准MPM中可能高估了稳定步长；因此建议使用约0.7的安全因子。他们还表明，当集中节点质量趋近于零时，临界时间步长趋近于零，例如当粒子即将离开节点的支持时。Bai和Schroeder [24]对具有二次和三次B样条插值的显式粒子在单元（PIC）方案进行了冯·诺伊曼稳定性分析。假设在静止状态下配置是周期性的、均匀填充的（每个单元中的粒子排列和属性相同、边界是周期性的以及布局是对称的），他们发现PIC使用二次B样条的稳定性限制约为声速限制的70%。Wallstedt和Guilkey [25]进行的时间收敛性测试中也观察到了0.7的减少，这与本研究的发现一致。

本工作旨在更深入地理解具有显式时间积分的MPM中的临界时间步长。通过关注单个计算周期，排除了粒子-网格映射引起的不稳定性。将基于粒子的公式中固有的积分误差视为对理想情况的扰动，在理想情况下，弱形式中的积分可以精确评估。分析还扩展到了GIMP和CPDI，这两种方法与标准MPM共享相同的拉格朗日网格框架和增量更新结构，但主要目的是减轻这些积分误差。首先计算单步能量增量，以了解影响稳定性的关键因素。然后对放大矩阵进行单步谱分析，得出稳定性条件，结果表明，在节点速度投影中使用集中质量可以减少相对于仅基于系统最大频率的临界时间步长。一个自由振动的、固定-自由的一维杆作为模型问题，用于验证和检查质量参数、粒子配置和时间积分参数的影响。从单单元模型中推导出了一个实用的临界时间步长公式，并且是保守的（即，它不会高估稳定步长）。对于具有集中质量的速度投影的一维问题，减少因子介于0.71（$\sqrt{0.5}$）和0.77（$\sqrt{0.6}$）之间，具体取决于粒子配置。基于单单元临界时间步长公式，提出了一种简单而有效的稳定化技术，并在理想化和真实的MPM模拟中进行了验证。

本文的其余部分介绍了标准的应力最后MPM算法，展示了单步稳定性分析，并报告了一个关于振荡杆的案例研究。结论在最后一部分进行了总结。

章节片段

材料点法

材料点法（MPM）[1] [2]使用材料点（粒子）表示连续体，并使用拉格朗日计算网格来评估运动和变形的演变。在每个分析步骤中，将粒子速度场投影到网格上，然后从柯西运动方程的弱形式中获得加速度场。然后计算粒子速度、位移和材料状态变量（例如应变和应力）的单步增量

单步稳定性分析

从全局意义上讲，数值稳定性要求初始数据中的小扰动只产生最终解中的小偏差。线性系统（或非线性系统的线性化模型）被广泛用于稳定性分析[31]。评估稳定性的一个实际方法是检查线性系统的每个数值更新步骤是否是能量守恒的或耗散的。因为材料点方法允许粒子之间的任意相对运动

案例研究和不稳定性缓解方案

为了更好地理解稳定性条件并识别影响数值稳定性的关键参数，研究了一根一维振荡杆。图1(a)展示了一根横截面积为A_r、弹性模量为E_r、质量密度为ρ_r、长度为L_r的杆的设置。杆的左端固定（$x=0$），右端自由（$x=L_r$）。初始速度定义为$v_{\text{init}}\left(x\right)=v_o\sin \frac{\pi x}{2L_r}$，对应于系统的第一个振荡模式。采用线性形状函数

结论

本研究全面检查了材料点法（MPM）的数值稳定性，确定了影响稳定性的基本因素，并提供了有效的缓解策略。通过基于能量和谱的单步分析，表明临界时间步长取决于动态质量参数（$\stackrel{?}{\gamma }$）、投影质量参数（$\hat{\gamma }$）、时间积分参数（α）以及粒子的空间分布，特别是

CRediT作者贡献声明

杨文嘉：写作 – 审稿与编辑，写作 – 原稿，可视化，验证，软件，方法论，调查，形式分析，概念化。黛博拉·L·苏尔斯基：写作 – 审稿与编辑，概念化。

利益冲突声明

作者声明他们没有已知的竞争财务利益或个人关系可能会影响本文报告的工作。