时间离散化及其MPM算法

上式是物质点法的半离散化形式，全离散方程是通过对其中描述的方程进行时间积分而获得的。

本期推送目录如下：

1. 集中质量矩阵
2. 节点速度（动量）的计算
3. 标准形式（USL）
4. 改进的更新应力（MUSL）
5. 优先更新应力（USF）

我们在第1节中提出了质量集中，以避免反转质量矩阵，这是一个耗时的步骤。与存储在节点处的速度场的 FEM 相反，在 MPM，网格速度在网格重置后的每一步都是无效的。因此，在时间步长 t 的开始，需要将粒子速度投影到网格节点，作为时间推进的起点。这个关键的步骤在2节中讨论。第3节介绍了第一个完整的 MPM 算法 USL (更新压力最后) ，它可能是最流行的MPM 算法。在第4节中给出了对 USL 和被称为 MUSL (修改的 USL)的轻微修改，以增强 MPM 模拟的稳定性。

1. 集中质量矩阵

从方程中可以明显看出。上式为了获得加速度，需要在每个时间步长求解线性方程组。这个系统的大小在3D中是    ，其中    ——网格的节点数量——可以是一个非常大的数字。为了避免这种情况，采用了集中质量矩阵。这个集总质量矩阵是一个对角矩阵，其对角项由下式给出

其中上式用于第二等式，FE形状函数的单位分割（PU）性质，    ，    用于第三等式。

集中质量矩阵不仅仅是计算上的简化；对于冲击载荷，它也给出了更好的结果。还应注意，集中质量矩阵的使用会导致能量耗散

由于动量方程是在网格节点处求解而不是在物质点处求解，所以必须在节点处满足质量守恒，所以可以得到：

证明了质量在网格节点处也是守恒的。

有了这个集总质量矩阵，就可以得到下面的时间常微分方程组（ODE）

对于所有节点I。

让我们用    表示模拟时间。然后，对于时间离散化，时域    被划分为    个时间步长，时间增量    。

因此，上述方程式必须满足每次实例    ，    。在时间上推进解的最直接的方法，即求解半离散方程，是一个显式的公式，其中的解是推进时间从t(即    )到    (或    ) ，而不解一个线性代数方程组。正向欧拉法就是这样一种方案。注意显式方案需要使用相当小的时间步长    来保证稳定性。

我们用上标t表示时间瞬间t的已知量，下标    表示下一个时间瞬间的未知量。

2. 节点速度（动量）的计算

在每个时间步的开始，粒子速度需要映射到节点。这一步骤不存在于有限元中，是必要的，因为网格是在一个时间步骤结束时重置和该步骤的节点速度丢失。更准确地说，节点动量是用形状函数从粒子映射出来的。

或者粒子动量被投影到网格节点

在接下来的讨论中，我们证明了动量投影，因为粒子比节点多。为了便于表示，让我们考虑一个由一个单元和两个节点组成的一维网格，其速度用    表示。此外，在这个细胞内，有两个粒子在    和    。

这个想法是为了最小化以下函数

这是一个加权最小二乘，其中的权重是粒子质量。在上述方程中，    是粒子的速度(而不是体积) ，    是粒子的质量。微分    与    和    的关系，使它们成为0,1得到

如果使用简短的符号    ，则上述方程可以改写如下

这是一个解 v1和 v2的线性代数方程组。可以看出，这个线性系统的系数矩阵正是一致的质量矩阵。为了避免这种系统的解决方案，使用了行和集总技术，其结果如下：

其可以被推广以获得下述等式，该等式结束了证明。

我们正在用上述公式检验线性动量是否守恒。用上述公式可以写成：

证明了线性动量在网格节点处也是守恒的(只要    做单位分区)。

尽管这种粒子速度投影几乎在所有 MPM 模拟中都得到了应用，但我们将在后面看到，它不能为任意粒子位置提供线性速度场的精确投影。

3. 标准形式（USL）

上式是用于获取节点加速度    的。

然后使用显式欧拉正演方法更新节点速度场，如下所示：

其中    表示时间    处的节点速度，这是已知的；    是时间增量。理论上，节点被移动到下式给出的新位置

请注意，这种节点位置更新很少实现，因为无论如何，网格将在下一个时间步骤的开始重置，就像大多数 MPM 实现所做的那样。然而，应该强调的是，网格重置不是一个必要条件，可以使用上述等式更新网格节点。直到网格失真或通过任何其他合适的栅格替换网格。这是可能的，因为所有信息都已存储在粒子中。

网格更新后，网格速度用于更新粒子状态，包括位置、速度、体积、变形梯度、应力等。我们首先讨论粒子位置和速度的更新，因为有不同的选项可能会混淆新加入的粒子。选项包括PIC、FLIP和混合PIC-FLIP。

在 PIC 方法中，粒子速度是通过网格总速度获得的。另一方面，根据 FLIP 方法，粒子速度是通过网格速度增量得到的。这些概括在以下公式中

由于PIC用网格速度代替粒子速度，信息是丢失的（粒子比网格节点多得多），因此PIC具有数值耗散（或者换句话说，对于没有耗散的问题，能量是不守恒的）。FLIP通过仅向粒子添加栅格速度增量来克服这一点。

PIC和FLIP的结合最早由Zhu和Bridson（2005）提出；Stomakhin等人（2013）在计算机图形学社区。它后来在工程界被使用，例如Leroch等人（2018）。在这种混合的PIC-FLIP中，PIC和FLIP的混合物用于粒子速度

其中    对应于FLIP，    对应于PIC。我们参考下图来说明    的影响。

图：杆的振动: PIC 与 FLIP

请注意，Leroch 等(2018)将粒子位置更新为    ，我们发现给出了与标准方法相似的结果。

粒子速度更新实际计算如下

最后，对粒子应力进行了更新。这在 MPM 文献中被称为更新应力最后公式(USL)。根据所使用的本构模型，可能需要计算梯度变形    、速度梯度    和变形速率    等。例如，我们可能需要计算粒子速度梯度，然后使用关系式    计算梯度变形，最后使用连续性方程    来确定更新的体积。这是典型的超弹性固体。它们是由

其中    是一个    的矩阵，其中的分量是    (在三维中)。在上述方程中，    是恒等矩阵，    是    矩阵，初始矩阵为    即    。我们增加了更新粒子密度的方程，例如需要计算状态方程。请注意，还有其他方法来计算更高的精度变形梯度，虽然更复杂。

对于三维问题，速度梯度    实际上计算为：

一维和平面应变/应力二维问题的简化是直接的。

对于次弹性本构模型，需要计算应变增量    并用它来计算应力增量    更新后的颗粒应力由下式给出：

对于复杂的本构模型，可能需要计算其他量来更新颗粒应力。由于材料点是拉格朗日函数，现有的应力更新算法主要由有限元社区开发，可以很容易地在 MPM 中重用。我们参考 Simo 和 Hughes (1998)的教科书; de Souza Neto 等人(2011)的细节。这就引出了算法1中给出的第一个完整的显式 ULMPM 算法。可以看出，这是一个非常简单的算法，可以直接进行编码。然而，它已经被用来解决许多具有挑战性的固体力学问题。请注意，    没有时间标签，因为它从来没有改变，以确保质量守恒。

有一些关于算法1的注释。首先，我们忽略了由于非零牵引力对外力的贡献。这是因为它比之后讨论的FEM更难处理。其次，我们假设使用了恒定的时间步长    。在实践中，经常采用不同的时间步长。最后，我们提出了所谓的动量公式。注意，这种动量公式非常常见，但它并不能提高MPM的稳定性。

小问题。在这个公式中有一个数值问题: 在    时，除法算符会产生无穷大的加速度，如果    的质量很小。反过来速度梯度也是无穷大的，这会破坏粒子的应力(Sulsky et al. 1995b)。当一个粒子非常接近一个只有一个粒子支撑的节点时，就会发生这种情况。这通常发生在靠近材料界面的节点上（见下图）。对于多体模拟来说，这些接口节点经常处于接触状态，因此接触算法必须仔细解决小节点质量问题。

图：质量几乎为零的故障节点导致无限加速度

为了更清楚地识别故障，我们转向上图中给出的示例：我们假设只有一个元素和一个粒子，质量由    表示，位置非常接近节点1。此外，我们假设旧的节点速度为零，也就是    ，节点质量由

因此，节点加速度：

其中    和    是节点力。更新的节点速度由下式给出

更新后的粒子位置和速度由下式给出（    ）

请注意，更新的粒子位置和速度是正常的，因为他们是平滑的形状函数。

接下来，计算应力更新所需的速度梯度

总之，问题在于使用    的计算速度梯度，但不会更新粒子速度和位置

截止技术。在这种技术中，引入了一个小的正阈值来解决小质量问题。因此，节点速度计算为

这个算法需要一个额外的参数(截止值)。然而，如何选择还不清楚。即使可以选择一个好的截止值，它也会产生一个不应该出现在系统中的不需要的约束。下面将介绍更先进的技术。

Wallstedt 和 Guilkey (2008)研究了许多用于 GIMP 的时间集成方案系列，包括 Runge Kutta，Runge-KuttaNystrom，Adams-Bashforth-Moulton (ABM)和 Prector-纠正 Newmark 方法。他们报告说，这些方法很少能够达到其正式的准确程度。MPM 方法不仅用于处理高度不连续和非线性的问题，而且该方法的空间误差往往压倒了暂时高阶方法可能提供的任何改进。他们还表明，中心差分格式，通常用于非线性有限元程序(Belytschko et al. 2000) ，与 USL 完全相同，但有一个关键的差异:初始化粒子速度为负半时步。

4. 上次修改的更新应力（MUSL）

在 Sulsky 等人(1995b)的 MUSL 中，更新后的粒子速度被映射回节点，以获得节点速度

因此：

在第二个等式中，形函数在分子和分母中的出现抵消了形函数的作用，从而解决了速度梯度很大的数值问题。将此应用到上图多体问题给出的例子中，现在给出了在问题节点处的更新速度。

上式的的第一个用来简化    。显然    不是无穷大，可以安全地用来计算速度梯度。算法2给出了最终的算法，称为 MUSL。这个算法的另一个名称是双映射 USL，因为粒子动量被两次外推到节点上ーー在步骤的开始和更新节点动量之后。我们用一个波浪线    来表示临时的网格速度(第14行)。MUSL 与 USL 的区别在第16-21行。

正是 Bardenhagen (2002)引入了术语更新压力最后(USL) ，并提出了更新压力首先(USF)算法，在 P2G 步骤中更新压力。他发现虽然 USL 是耗散的(例如，遭受数值耗散) ，USF 保存能量很好。奈恩(2003)分析了 USF 和 MUSL 的二维弹性振动问题。他发现 USL 是非常不稳定的，MUSL 方法缓慢地消耗能量，而 USF 方法缓慢地增加能量。Wallstedt 和 Guilkey (2008)使用人工解的方法测试了 GIMP 与 USF 和 USL 的时间和空间收敛性。结果表明，USL 在稳定性和收敛性方面具有优势。

到目前为止，我们只提出了显式的动力学MPM。开发这些算法是非常简单的，因为 MPM 与更新的拉格朗日有限元非常相似。这就是为 ULFEM 开发的几何和材料切线矩阵的表达式可以很容易地重用，见 Belytschko 等人(2000)。Iaconeta 等人(2017)给出了一个隐式动力学 MPM 的清晰表述。

显式 MPM 代码也可用于准静态问题。尽管使用显式代码进行简单的静态模拟效率不高(因为长时间的模拟需要很多步骤) ，但显式代码是挑战静态模拟的唯一选择，因为隐式求解器可能会崩溃(例如，求解器不会收敛)。可以添加全局和局部阻尼以减轻静态模拟中应力波的影响(Al-Kafaji，2013)。

5. 更新应力优先（USF）

USF 是最后介绍的方法(Bardenhagen 2002)。在 USF，压力更新在时间步骤的开始，而不是在结束(见算法3)。根据奈恩(2003) ，USF 是另一种缓解4节讨论的小质量问题的方法。

Bardenhagen (2002)发现 USF 比 USL 更节省能量。