物点法的理论、实现与应用系列解说之物质点法简介
本文的目录如下：
介绍
1.1 计算机科学与工程(CSE)
1.2 实验在计算机科学与工程中的地位
关于物质点法(MPM)的简介
2.1 拉格朗日粒子和欧拉网格
2.2 基本的物质点法(MPM)算法
2.3 物质点法的优点和缺点
2.4 目前的物质点法的表达式
2.5 多物理场的物质点法
2.6 接触
2.7 断裂
2.8 流体和气体
2.9 物质点法和其它方法的比较
2.10 耦合物质点法和其它的方法
全文字数1w+，非常烧脑，请各位读者耐心阅读。
1. 介绍
1.1 计算机科学与工程（CSE）
本期推送的主题属于计算科学与工程（CSE）的范畴。CSE是一门相对较新的学科，涉及计算模型的开发和应用，通常与高性能计算相结合，以解决工程分析和设计（计算工程）以及自然现象（计算科学）中出现的复杂问题。CSE被描述为仅次于理论和实验的“第三种发现模式”。在 CSE 领域中，这些是解决问题的步骤。首先，选择或发展一个最能描述该问题的数学模型。模型开发的这一步是由具有足够数学技能的人手工完成的。大多数的数学模型都是用微积分来建立的。因此，它们是不适合于数字计算机的连续模型。其次，推导出这个数学模型的计算模型。计算模型是数学模型的近似形式，是可以用计算机求解的离散形式。第三，这个离散模型是用一种编程语言(过去是 Fortran，现在是 C + + 和 Python)实现的，以便有一个计算代码或平台。对于固体力学，流行的计算平台是 Abaqus 和 LS-Dyna。最后，这些计算平台被用来进行计算机模拟或计算机实验。
计算机模拟不仅有助于解决过于复杂而无法解析解决的问题，而且越来越多地取代了昂贵和耗时的实验。此外，它们可以在难以或不可能实现实验可视化的空间和时间尺度上提供大量信息。最后，模拟在预测尚未创建的配偶和结构的行为方面也有价值；实验仅限于已经创建的材料和结构。
我们想简要地讨论一下为什么微积分在科学和工程学中占有如此重要的地位。这一切都始于伽利略和牛顿的著作，他们发现自然定律可以被数学，特别是微积分，很好地描述。如果天体的运动可以用数学来模拟，那么将其应用于人类问题就是合乎逻辑的。这正是伯努利兄弟、欧拉、拉格兰奇、柯西等天才在300年前所做的事情。这些人发展了偏微分方程，可以模拟流体、气体和固体的变形等一系列现象。正是这些偏微分方程描述的模型把人类送上了月球，给了我们手机，电脑，无线电。或者电视。或者给准妈妈做超声波，或者给迷路的旅行者做 GPS。
微积分分为微分和积分两部分。当我们从有限走向无限时，后者是有趣的。为了数值求解，我们反过来做: 从无穷回到有限。我们用一个只包含有限个自由度的网格来代替一个具有无限个自由度的实体。
1.2 实验在计算机科学与工程中的地位
可以肯定的是，计算模型需要通过实验来获得模型中使用的参数。为了强调实验在科学和工程中的重要作用，我们考虑了 Boyce 等人(2016)的有趣文章，其中提出了 Sandia 断裂对计算断裂社区的挑战。这项挑战涉及模拟复几何钢样品的断裂。只有一个最低限度的实验数据(拉力试验钢券)提供给分析师。参与挑战的不同研究小组使用了所有现有的骨折模型，没有一个与实验匹配。
因此，模拟仅仅是不够的，因此，对于任何问题都应该采用联合实验模拟方案。有趣的是，大约70年前创造了“有限元分析”这个词的 R.W. 克拉夫停止了这个方法的研究，转而进行实验(克拉夫1980)。
2. 物质点法的简介
材料点法是粒子单元法(PIC)的最新发展之一。第一个 PIC 技术是在20世纪50年代早期由哈洛(1964)和哈洛(2004)在洛斯阿拉莫斯国家实验室发展起来的，主要应用于流体力学。1986年，Brackbill 和 Ruppel 通过引入 FLIP-流体隐式粒子方法(Brackbill and Ruppel 1986; Brackbill et al. 1988)克服了第一批 PIC 受到过度能量耗散的困扰。在计算机图形学方面，PIC/FLIP 已经成为流体模拟的行业标准方法(Zhu and Bridson 2005)。FLIP 后来被新墨西哥大学的 Sulsky 及其同事(Sulsky 等，1994,1995b)修改并定制用于固体力学中的应用，此后被称为材料点方法(Sulsky 和 Schreyer 1996)。
在 FLIP 中，应变和应力存储在细胞中心。然而，在 MPM，它们是由粒子自身携带的。因此，MPM 颗粒携带材料的全部物理状态，包括位置、质量、速度、体积、应力、温度等。注意，在 PIC 中，粒子只携带位置和质量。
MPM 建立在两个已经在 PIC 中使用的主要概念之上，即使用携带物理信息的拉格朗日材料点，以及用于连续场(即位移场)离散化的背景欧拉网格。对拉格朗日和欧拉描述的简短描述，请看图1。
图1：拉格朗日描述（上）与欧拉描述（下）。在拉格朗日描述中，网格附着在实体上，因此它在实体的变形过程中变形。网格中的每个点始终只与一个材质点相关联，因此可以轻松地对依赖历史的材质进行建模。实体边界也得到了很好的定义。但是，栅格可能会扭曲。另一方面，欧拉网格在空间中是固定的，材料流经网格。网格变形从未发生
2.1 拉格朗日粒子与欧拉网格
在MPM中，连续体由在整个变形过程中跟踪的有限组np拉格朗日材料点（或粒子）离散化。粒子和物质点这两个术语将在本书中互换使用。在原始MPM中，粒子表示的子区域没有明确定义。只跟踪它们的质量和体积。在先进的MPM配方中，如GIMP或CPDI，这些子区域的形状被跟踪。每个物质点有一个关联的位置 xtp (p = 1,2，... ，np) ，质量 mp，密度 ρp，速度 vp，变形梯度 Fp，柯西应力张量 σp，温度 Tp，以及任何其他内部状态变量必要的本构模型。总的来说，这些物质点提供了连续体的拉格朗日描述。由于每个物质点在任何时候都包含一定量的质量，质量守恒就自动得到满足。
Sulsky开发的原始MPM实际上是一个更新的拉格朗日格式。对于这种MPM，模拟物体在变形过程中占据和将占据的空间由一个称为背景网格的网格离散，在该网格中求解动量平衡方程。另一方面，在总拉格朗日MPM（de Vaucorbeil et al.2020）中，背景网格仅覆盖其参考配置中物体覆盖的空间。我们参考图2，了解ULMPM和TLMPM的笛卡尔网格上叠加的材料点的图形说明。栅格是固定的，粒子在栅格上移动（图3）。
网格的使用有以下好处。首先，它消除了直接计算粒子-粒子相互作用的需要，使该方法具有相当的可扩展性。其次，通过这种背景欧拉网格可以很容易地处理碰撞（事实上，该方法固有的是防滑、不 穿透的接触）。第三，动量方程是在网格上求解的，由于网格点比粒子少得多，这是一种非常有效的替代方法。大多数情况下，由于效率的原因，在整个模拟过程中使用固定的规则笛卡尔网格。
图2：MPM离散化：空间由背景网格离散化，背景网格可以是笛卡尔网格或非结构化网格（未显示），而固体则使用粒子离散化。更新后的拉格朗日MPM网格覆盖整个变形空间，而总拉格朗日MPM栅格仅覆盖初始构型
图3：MPM离散化：网格是固定的，粒子在其上移动
2.2 MPM的基本算法
MPM最初是为了解决快速瞬态冲击固体力学问题而开发的（Sulsky等人，1994）。因此，MPM是使用显式求解器开发的，对于此类问题，它比隐式求解器更有效。然后将该方法应用于负载率较低的其他应用。对于这些问题，隐式求解器更适合。由于显式MPM算法比隐式算法更简单，因此在下文中，使用显式求解器提出了更新的拉格朗日MPM算法。隐式MPM公式也进行了讨论。根据更新后的拉格朗日MPM，只需稍作修改即可获得总拉格朗日MPM。
图4：物质点法：计算步骤由四个步骤组成：（1）P2G（粒子到网格），其中信息从粒子映射到节点；（2）网格更新，其中为节点求解动量方程，（3）G2P（栅格到粒子），其中更新的节点随后被映射回粒子以更新它们的位置和速度，以及（4）栅格重置，其中栅格被重置。ULFEM中没有虚线框中的操作
图4给出了一个典型的显式 ULMPM 计算循环。我们参考表1.1中的表示法列表(并非详尽无遗)和表1.3中的缩写。这个算法有点不成熟，因为有些术语还没有被精确地定义，但是在这个阶段提出来是为了给出一些观点。第一步是将信息从粒子映射到网格(P2G) ，因为网格在每个周期都被重置。其次，在网格节点上求解离散的动量方程(网格更新)。然后，更新颗粒的位置、速度、体积、密度、变形梯度、应力和所有相关内部变量(G2P)。最后两个步骤等价于更新的拉格朗日有限元法。因此，在 Gupta 等人(2011)中声明 MPM 使用欧拉内核是不正确的。最后，将网格重置为其原始状态。由于这种网格重置不会发生网格变形，因此 MPM 是解决大变形问题的一种很好的方法。请注意，网格不需要在每个时间步重置。例如，重置可以每 N 个时间步骤完成。N 可以是任意大小。网格也永远不能被重置(Guilkey等人2006）。此外，还可以使用全新的网格。但据我们所知，这还没有在任何代码中实现。
由于拉格朗日粒子和背景网格的结合，很难对MPM进行精确分类。在我们看来，当MPM求解弱形式的动量方程时，它可以被视为一种无网格伽辽金方法，类似于EFG、RKPM和OTM。MPM与其他Galerkin MM的区别在于构造形状函数的容易性。事实上，它们是在固定的欧拉网格上定义的简单有效的多项式。请注意，大多数无网格形状函数是在节点云上定义的计算成本高昂的有理函数。当网格不固定时，MPM与OTM非常相似（反之亦然）。不同之处在于OTM采用了最大值近似（Iaconeta等人，2017）。与OTM相反，MPM使用背景网格，如果不固定，将产生网格纠缠问题，类似于更新的拉格朗日有限元。
2.3 物质点法的优点和缺点
MPM 的优点包括:
1. 无网格纠缠问题;无误差地通过物质点的运动使物质特性平移;
2. 该方法自动采用防滑、防渗透接触算法，即不需要额外的计算代价;
3. 利用欧拉网格的背景，直接有效地处理多体的摩擦接触;
4. 适用于以图像为基础的模拟，因为将图像转换为MPM 模型;
5. 与现有的无网格方法相比，MPM 的计算机实现非常简单。通过对计算域的分解，可以方便地为并行分布式存储计算机编写 MPM 算法程序;
6. 利用现有经过充分研究的有限元算法，因为MPM 与 FEM算法的相似性。
7. 适用于难以转换成优质有限元网格的非常复杂的几何问题。在这方面，MPM 与浸入边界方法(Mittal 和 Iaccarino 2005; Schillinger 等，2012)及其最近的变体如有限细胞方法(Parvizian 等，2007)和切割 FEM (Burman 等，2015)有相似之处。所有这些方法，不管名称如何，都将任意形状的实体嵌入比实体大的立方体(3D)中。立方体由规则的结构网格组成，该网格具有光滑的 Ck 基函数，可以是 B 样条或 T 样条。
除了MPM提供的这些优点之外，与任何数值方法一样，它也有自己的一系列缺点：
1. 大的内存占用空间，因为网格必须覆盖主体占据的整个区域；
2. MPM的形式化分析（收敛性、误差和稳定性）极其困难；
3. 与有限元法相比，边界条件的实施很困难；
4. 精度低于FEM，因为材料点通常不位于数值积分的最佳位置。
第一项只应用于 ULMPM，而不是 TLMPM，因为后面的网格只覆盖初始的未变形配置。由于颗粒的不规则分布和颗粒相对于网格的相对运动，MPM 的形式分析是非常困难的。如果要进行这样的分析，需要许多假设来使分析易于管理，即，一维线性弹性材料，粒子不会从一个细胞移动到另一个细胞(York et al. 1999)。边界条件难以执行的原因是缺乏对边界的明确表示。但是，它仍然比其他一些方法(例如:SPH) (Raymond et al. 2018).然而，细胞中的广义粒子(GPIC)不存在这个问题(Nguyen et al. 2021)。最后，低精度只适用于小型和中等大型变形问题。
2.4 目前的物质点法的表达式
文献中可用的不同MPM变体之间的主要差异来自于（i）不同的网格基函数，（ii）不同的时间积分方案，（iii）不同类型的网格（笛卡尔或非结构化）和（iv）更新的拉格朗日量或总拉格朗日量公式的使用。此处的讨论仅限于ULMPM。
使用不同的网格基函数。在最初的 MPM，在这里被称为标准/传统的 MPM，网格基础函数 φI (x)是线性帽子函数。因为这些函数只有 C0，所以标准 MPM 在粒子穿过细胞边界时会遇到所谓的细胞穿越不稳定性。为了解决这个问题，Bardenhagen 和 Kober (2004)引入了广义插值材料点(GIMP)方法。在 GIMP 中，与最初的 MPM 相反，粒子不是点状的，而是有限的范围，见图1.15。GIMP 得到的网格基函数是 C1光滑函数，类似于 Steffen 等人(2008b)使用的 B 样条函数。具有 B 样条的 MPM，命名为 BSMPM，现在非常流行(Stomakhin 等，2014a; Tielen 等，2017; Gan 等，2018)。
与粒子从一个细胞移动到另一个细胞有关的另一个问题是固体的数值断裂非物理分离（图5）。当两个原本相邻的粒子分离足够高的距离时，就会发生这种情况。他们不再互相影响。此距离取决于所使用的网格函数的类型。由于这种数值断裂，很难确定预测的断裂是物理的还是数值的，见图6
图5：MPM 中由于不稳定引起的数值断裂。这是用标准 MPM (线性基准)模拟两个橡胶环的碰撞。模型中没有包括断裂模型
图6：球形弹丸穿透钢板：局部Johnson Cook损伤模型。具有线性加权函数的ULMPM
尽管在GIMP中，粒子域在一定程度上被跟踪，但它们之间的间隙仍然存在，导致在任意变形下精度较低。根据我们的现有情况，转换粒子域插值（CPDI）是GIMP的最新发展，由Sadeghirad等人提出。（2011），Sadeghirard等人（2013）解决了这个问题。在CPDI中，粒子分别在2D和3D中建模为四边形和四面体。因此，对于任意载荷，变形的固体是无间隙平铺的，见图1.15c。CPDI非常善于处理极端拉伸变形，而不会出现数值断裂，并且能够忠实地表示复杂的几何形状。然而，它也不能免除问题。首先，CPDI存在网格失真问题，这违背了传统MPM的精神，参见Wang等人（2019）。其次，并行化在CPDI中更为复杂（Homel等人，2016），因为粒子可能具有足够大的影响域，可以跨越许多不同的CPU域。
MPM的最新进展之一是Nguyen等人（2021）的细胞中的广义粒子。它将CPDI扩展到任何元素类型和背景网格，并且它不会受到CPDI的网格失真的影响，因为它采用了计算力的拉格朗日公式。我们在图7中给出了 GPIC 的一些应用。
图7：广义单元内粒子（GPIC）方法：嵌入欧拉网格的TL有限元网格。GPIC解决了基于粒子的MPM的许多问题：（i）更好地表示固体边界，（ii）狄利克雷和诺依曼边界条件的无缝实施，（iii）材料界面的无缝处理，（iii）更高的效率和（iv）更好的应力场（Nguyen等人，2021）
使用不同的时间积分方案。MPM最常见的实现方式使用显式向后欧拉积分方案。其他显式集成方案也得到了不同程度的成功使用。Sulsky等人（2007）采用了交错中心差分积分方案来模拟海冰动力学，后来将其与Wallstedt和Guilkey（2008）最常见的实施方案进行了比较。他们的结论是，这两种方案之间只有很小的区别。
隐含的MPM也存在。显式时间积分方案对于短时间动态问题具有易于实现、高效和稳定的特点。然而，它们需要遵守Courant稳定性条件dt<dx c，其中dx是网格间距，c是材料中的声速。对于低应变速率和准静态问题，例如金属轧制、镦粗和机加工等制造问题，这种情况很严重，因此使用隐式时间积分方案更有效。在隐式mpm中，可以显式地形成雅可比矩阵，也可以采用无矩阵求解器。在下文中，我们回顾了动力学和准静态问题的隐式MPM。
第一个使用隐式方案进行MPM的是Cummins和Brackbill（2002）。他们将其应用于模拟涉及摩擦接触的颗粒材料的准静态载荷。众所周知，隐式集成方案成本高昂。因此，为了减少计算时间（通过避免构造正切矩阵），他们采用了无矩阵的Newton-Krylov算法。Sulsky和Kaul（2004）报道了类似的方法，Love和Sulsky（2006b）、Love和Sulsky（2006a）将其扩展到超弹性塑料材料。
这些论文采用了三场公式来避免体积锁定（第9.6节提供了关于体积锁定的讨论）。在这些文献中详细讨论了MPM的线动量守恒和角动量守恒。Nair和Roy（2012）为GIMP实现了这种无矩阵的隐式动力学。Sanchez（2015）等人将这种无矩阵方法扩展到了准静态问题。展示了不需要复杂本构模型通常难以获得的一致材料切线的优点。
相比之下，Guilkey和Weiss（2003），Wang等人（2016）明确形成了切线刚度矩阵，并使用Newton-Raphson方法和众所周知的Newmark积分方案及时求解平衡方程。他们报告说，时间步长比显式MPM中使用的时间步长大数百倍。此外，一致切线的使用允许时间步长比Cummins和Brackbill（2002）、Sulsky和Kaul（2004）的无矩阵公式的时间步长大得多。我们参考Iaconeta等人（2017）对Guilkey和Weiss（2003）算法的详细介绍。
笛卡尔网格与非结构化网格。在MPM中，通常使用统一的笛卡尔网格。这消除了在粒子网格交互过程中对计算成本高昂的邻域搜索的需要（即，粒子将其数据映射到哪些节点，反之亦然）。此外，使用笛卡尔坐标系比使用非结构化网格更容易开发光滑的Ck（k≥1）基函数。
另一方面，在岩土工程界，通常采用非结构化网格。Wieçkowski等人（1999）和Wieçkowski（2004）是第一个将非结构化网格用于筒仓卸料应用的人。此后，德国斯图加特大学和荷兰代尔夫特理工大学的相关研究小组在后来的工作中使用了非结构化网格，例如Beuth等人（2011）、Jassim等人（2013）。与笛卡尔网格相比，在非结构化网格中搜索哪个元素包含给定粒子并非易事，而且非常耗时。然而，非结构化网格的使用有助于执行复杂的边界条件（即曲面上的边界条件）。
在与GIMP和CPDI平行的研究中，Zhang等人（2011）提出了非结构化网格的双域MPM（DDMPM）。DDMPM的动机是难以为非结构化网格开发C1函数来缓解单元交叉问题。其基本思想是将粒子应力映射到网格节点，然后对其进行插值，以在域的任何点获得连续的应力场。这种双重映射过程会产生更平滑的渐变。最近，de Koster等人（2019）使用Powell和Sabin函数在非结构化网格上开发了C1基函数（Powell and Sabin 1977）。
高阶MPM。尽管MPM及其先前提出的变体已经进行了令人印象深刻的模拟，但对简单问题的方法的严格分析表明，收敛速度很差（很少是二阶），如第9章中进一步解释的那样。为了提高收敛阶数，Wallstedt和Guilkey（2011）提出了固体的加权最小二乘MPM，Edwards和Bridson（2012）提出了流体的移动最小二乘（MLS）MPM。Sulsky和Gong（2016）提出了一种改进的MPM（iMPM）。在iMPM中，除使用MLS的从粒子到节点的速度映射外，所有映射都使用hat函数。此外，使用一点正交，即，正交点是小区中心。这消除了细胞交叉的不稳定性。
他们证明了iMPM的二阶收敛性，但仅适用于中等精细的网格。对于非常精细的网格没有观察到收敛。类似的想法可以在Wobbes等人（2019）中找到。这种高阶收敛仅在1D问题中得到证明。最近，梁等人（2019）也使用一点正交，但不是通过MLS，而是使用额外的交错网格来重建细胞中心数据
MPM 的另一个稳定性问题是所谓的零空间问题ーー非零粒子值的映射可以导致零节点值。这个问题产生于粒子数和节点网格点数之间的差异。零空间问题可能是导致 iMPM 不收敛的罪魁祸首。解决这一问题的方法在格里顿和伯津斯(2017) ，特兰和索沃斯基(2019)。
表1给出了一个表格，展示了常见 MPM 变体的特征。请注意，即使 iMPM 仅应用于标准 MPM，它也可以与其他 MPM 配合一起使用。
表一：常见MPM变体的总体特征
2.5 多物理场的物质点法
尽管MPM最初是为力学问题开发的，但它也已扩展到模拟多物理问题。这种模拟包括求解一个以上偏微分方程的耦合系统。例如，对于热机械过程，运动方程与Chen等人提出的能量平衡方程（或热扩散方程）相耦合。（2008），Nairn和Guilkey（2015），Fagan等人（2016），Tao等人（2016年），Gritton等人（2017），Tau等人（2018），Leroch等人（2018年）。他们表明，当使用显式MPM求解器时，使用交错求解器可以直接合并热扩散。在这些工作中，Fagan等人（2016）证明了MPM能够模拟后期搅拌摩擦焊（FSW）。FSW是一种最近出现的复杂的热机械工艺，用于连接不同的材料，但仍不太清楚。Grit-ton等人（2017）首次报道了硅阳极变形的化学/机械MPM耦合模拟
2.6 接触
当两个可变形实体相互接触时，就会发生接触。它是理解许多工程问题的关键因素，如桩土相互作用、磨损、金属成形过程等。接触的模拟仍然具有挑战性，并且已经开发了各种算法（例如，注意到段接触、段到段接触，惩罚和拉格朗日乘子法）（Benson 1992）。然而，由于使用了背景网格，在MPM中要容易得多。事实上，它的开发主要是为了处理接触问题
无滑动无穿透接触是MPM固有的，即固体的接触是自动处理的，而无需任何额外的数值处理。这是由于使用了单值速度场来更新材料点的位置。MPM的这种自动无滑动接触能力已被用于许多涉及复杂接触的工作中。例如，Bardenhagen等人（2005年）、Brydon等人（2005）使用GIMP模拟泡沫微观结构压缩至完全致密化。这是一个具有挑战性的问题，涉及离散复杂微观结构、模拟大变形和多重接触。这项工作表明，粒子方法非常适合由于不需要贴体网格，因此MPM25模拟具有复杂几何形状的实体的简介。奈恩（2006）应用MPM对木材的横向压缩和致密化进行了研究。刘等人（2015a）研究了高速冲击下的蜂窝夹芯板。最近，Sinaie等人（2019）模拟了涉及许多接触的细胞结构的大变形响应。
相反，不同物体之间的摩擦接触或单个固体内的自接触需要修改标准MPM算法。
多体摩擦接触。Wieçkowski等人可能提出了MPM中的第一种库仑摩擦接触算法。（1999）在孤岛放电模拟中模拟可变形的身体与刚性壁之间的接触。同时，约克等。（1999），约克等。（2000年）观察到，在MPM模拟中，有时将两个物体在应分开时毫无保证地“粘在彼此之间”。为了减轻这个问题，他们提出了MPM的第一个接触算法，该算法允许接触机构彼此释放。但是接触仍然没有滑倒。Bardenhagen等。（2000），Bardenhagen等。（2001）将约克的接触算法扩展到摩擦接触，并将其用于模拟颗粒材料中晶粒之间的相互作用。Bardenhagen的算法（称为多物质或多体触点）非常有效，因为它在谷物的数量中是线性的，并且可以分离，滑动和滚动。该算法的基本思想是修改接触节点的速度（即，从多个材料/身体接收贡献的速度）来解释碰撞。该算法非常受欢迎。例如，Nairn（2013）将其与修改一起模拟了不完美（剥离）接口。Lemiale等人给出了Bardenhagen算法的改进版本。（2010年）以建模一种称为均等信道角压力技术的金属形成过程。这项工作强调了在刚性和可变形物体之间接触的接触处理。所有先前的应用都仅限于库仑摩擦。最近，Nairn等人。（2018）概括了接触的MPM方法，以处理任何摩擦法律，并以粘附或速度依赖性摩擦系数给出摩擦的示例。
在岩土工程领域，Beuth（2012）提出了准静态MPM的首次接触算法。他们采用了有限元法中常用的零厚度界面单元来模拟土壤和刚性表面之间的接触。
该技术仅适用于在模拟开始之前已知接触表面的情况。Bardenhagen的算法在Al-Kafaji（2013）中使用显式动力学MPM代码对粘性接触进行了修改。马等人（2014）提出了另一种适用于岩土工程的Bardenhagen接触模型修改（涉及刚性结构元件和土壤之间的接触）。
由于大多数MPM变体提供了噪声应力场，de Vaucorbeil和Nguyen（2021b）提出了TLMPM的接触模型。该方法产生更平滑的应力场，这对于例如损伤/断裂问题至关重要
刘和孙（2020）提出了一种隐式MPM，它采用水平集方法来表示固体边界，并采用iMPM来具有更好的求积精度。
摩擦自接触。Homel和Herbold（2017）是第一个提出摩擦自接触算法的人。其基本思想是自动和动态地检测接触节点。这是通过定义潜在接触表面的标量场（实际上是其梯度）来实现的。同样的想法（命名为DFG）也被应用于压裂，并产生了压裂和压裂后摩擦接触的第一个模型，即裂缝面之间的接触。
2.7 断裂
断裂的定义是将一个固体分离成两块或多块。自Griffith（1920）、Irwin（1957）的开创性著作以来，人们对这种现象进行了广泛的研究。他们发展了断裂力学领域，该领域在许多工程学科，特别是航空航天工程中有着巨大的应用。
对固体中裂纹的萌生和扩展进行建模是一个具有挑战性的问题，因为需要跟踪位移场不连续的一组不断演变的内表面。基本上，裂缝建模有三种方法：不连续、连续和混合连续-不连续方法（Rabczuk 2013）。我们还参考了Boyce等人（2016）的一篇有趣的文章，该文章对模拟韧性断裂的多种不同技术的可预测性进行了分析。
在不连续方法中，裂纹被明确表示（遵循断裂力学理论）。然而，在连续方法中，裂纹没有明确建模（采用损伤力学理论（Kachanov 1958；Lemaitre和Chaboche 1994））。在后面的方法中，通过用于降低应力的坝龄变量来处理裂缝的存在。图1.19显示了使用这两种方法进行的裂缝模拟。
还有另一种方法——组合连续-不连续方法。断裂的连续-不连续组合方法将这两种方法结合在一起，顾名思义。它的开发主要是为了克服连续断裂模型的缺点，即无法捕捉到真正的材料分离。
非连续方法。在 MPM，明显的裂缝已被模拟为强烈的不连续性，通过允许多个速度场在节点的支持被裂纹切断。这类似于FEM中的重复节点。裂纹可以使用LEFM（线性弹性断裂力学）（Nairn 2003；Tan和Nairn，2002年；奈恩2007a；郭和奈恩2004；吉拉伯特等人，2011年；Wang等人2005）或内聚带模型（Daphalapurkar等人2007；Nairn 2007b）或两者的组合（Nairn 2009；Bardenhagen等人2011）。显式裂缝模拟在计算上是昂贵的（因为需要跟踪演变的表面），并且对于大规模模拟来说是禁止的。此外，实现是复杂的，特别是对于复杂的裂纹模式。这方面的最新进展是Moutsanidis等人的工作。（2019a）通过修改网格基函数引入了裂纹的单一速度场。
批注：尽管这些工作大多涉及涉及小变形的模拟，但我们认为不连续方法用于大变形没有任何限制。对于小（或中大）变形断裂力学问题，我们不认为MPM比FEM更好。这一观察结果得到了第一次桑迪亚骨折挑战的结果（Boyce等人，2014）的支持，其中一个团队使用了MPM（Yang等人，2014年），而MPM并没有给出最佳解决方案。
连续的方法。在MPM中，使用各种连续方法对裂缝进行了建模。它们都有一个共同的特点：不需要像不连续方法那样明确地表示裂纹表面。对于渗透问题，通常使用基于应变的失效准则和颗粒侵蚀来处理裂缝。也就是说，当粒子满足特定的基于应变的断裂条件时，它被设置为失效。当这种情况发生时，其偏应力设置为零，但它仍然是模拟的一部分，其质量对材料的整体惯性有贡献。因此，质量守恒被强制执行。Ionescu等人（2006）也使用了基于应变的失效与颗粒侵蚀相结合的方法来研究低速弹丸穿透软组织的失效。他们报告了精细网格的收敛结果，但没有提供定量证据。Huang等人采用了具有颗粒侵蚀的弹塑性材料。（2011）对薄钢板的穿透和厚铝板的穿孔进行了建模。研究了结果对网格大小和粒子数量的敏感性，一些量显示出收敛性。
对于低应变速率和准静态问题，通常使用标准本构模型。例如，Alonso和Zabala（2011）中使用了软化莫尔-库仑塑性模型来模拟Aznalcóllar大坝的破坏。最近，Raymond等人提出了一种与Grady-Kipp损伤模型相结合的软化Drucker Prager塑性。（2019）对骨料材料的失效进行了建模。他们报告了网格收敛的结果，尽管没有对软化进行特殊处理。Homel和Herbold（2017）采用了具有线性应变软化的Rankine损伤模型。
Schreyer等人也在静态和低应变速率动态问题上。（2002），Sulsky和Schreyer（2004）使用了弥散裂纹模型（他们称其为去粘聚性），请参见例如Rots等。（1985），rots（1991）。Chen等。（2002）通过使用基于应变的损伤扩散方程与拉伸损伤模型相结合，改善了该模型。类似的作品包括Shen和Chen（2005），Shen（2009），Yang等。（2012），Yang等。（2014）。MPM中这种剥夺模型的实施细节直到最近才在Sanchez（2011）中给出，他证明了结果是偏见的。也就是说，仅在裂纹方向与网格细胞系的方向重合的情况下，才能实现准确的溶液。桑切斯（Sanchez）进一步评论说，对MPM中弥散裂纹模型的研究将是浪费时间。
在最初的损伤力学框架中，损伤是使用局部变量（如应力和应变）来确定的。这已被证明可以创建依赖于网格的软化。为了缓解这个问题，可以使用非本地模型。所谓非局部，我们的意思是使用非局部变量来确定损伤，例如附近的平均应力和应变。这个想法是在方程中引入一个长度尺度。Burghardt等人（2012）是第一批在MPM中实现非局部塑性模型的人。
非局部损伤也可以使用相场变分坝龄/裂缝（PFF）理论进行模拟，例如参见Bourdin等人（2008）、Miehe等人（2010）、Wu（2017）、Wu和Nguyen（2018）以及Wu等人的最新综述。（2019）。Kakuris和Triantafyl-lou首次在MPM代码中实现了PFF（2017a，b）。不幸的是，当时只证明了涉及没有接触的小变形的问题。最近，Cheon和Kim（2019）也提出了类似的方法。它们使用自适应网格细化和粒子分割来捕捉相位场的梯度。在计算机图形学界，Wolper等人最近提出了一个类似的想法。（2019），在该界通常可以看到令人印象深刻的骨折模拟。另一个变分断裂模型，Schmidt等人的特征侵蚀模型。（2009）被Zhang等人使用。（2020）用于模拟小变形脆性断裂。
图8显示了TLMPM与非局部Johnson-Cook损伤模型相结合的应用，以模拟子弹穿透Weldox钢板引起的韧性断裂。
图8：初始速度为250m/s（有颗粒侵蚀）时，射弹与目标之间碰撞的时间演变。颜色显示每个粒子中存在的损坏程度。为了便于可视化，只显示了集中在撞击区域周围的模拟区域的一小部分（de Vaucorbeil等人，2022b）
备注6：随着断裂变分方法（或相场断裂模型）的发展，可以用传统的连续体有限元对固定网格上固体中复杂裂纹网络的萌生和扩展进行建模。因此，无网格方法更适合于断裂似乎不再令人信服的时候。然而，拥有一种可以处理碎片之间大型变形，接触，断裂和失败接触的方法非常有吸引力。即使必须在实现此目标之前必须做更多的工作，MPM也可能是一种方法。
2.8 流体和气体
尽管MPM是为应用于固体力学而开发的，其主要应用领域仍然是固体力学，但它也已应用于流体力学问题。这可能是因为需要仅使用MPM来模拟流体-结构相互作用（FSI），以避免将MPM固体解算器与流体解算器耦合，这不是一项容易的工作。在这里，我们只讨论流体的MPM算法。
使用MPM处理流体（和气体）基本上有两种方法。第一种是所谓的弱可压缩MPM。除了使用流体的本构模型外，这种方法几乎与固体MPM求解器相同。这是通过使用人工状态方程（EOS）描述流体压力和密度之间的关系来实现的。这种流体MPM求解器首次出现在York等人（1999）中，并在York等（2000）、Tran等人（2010）、Ma等人（2009b）中用于气体动力学问题。它也被用于流体流动问题（Li等人，2014；Mast等人2012）。赵等人（2017）、Vargas等人（2018）、Sun等人（2018年）介绍了使用该流体MPM解算器对SPH和MPM进行的流体力学对比研究。总体而言，已经表明MPM和SPH预测非常相似，但MPM计算时间较小
弱可压缩流体 MPM 存在两个主要问题。首先，如果使用一个显式的时间求解器，非常小的时间步骤是必要的，因为需要使用一个非常大体积模数。第二，发生压力振荡。要解决这些问题，第二种方法是真正不可压缩的 MPM。Stomakhin 等人(2014b)在计算机图形学社区是第一个提出这个方案的人。他们使用了柯林算子分裂法(柯林1968)。在工程界,zhang 等(2017) ，Kularathna 和 Soga (2017)提出了类似的方法。
2.9 物质点法和其他方法
我应该使用哪种无网格方法？这是每个新手每次需要解决一个大的变形问题时都要面对的问题，如果他们没有被迫使用任何特定的方法的话。可以肯定的是，没有人能够回答这个问题，我们也不能。尽管如此，对于这些初学者，我们在此讨论 FEM，SPH，离散元法(DEM) ，Galerkin MM，浸入边界方法(Mittal 和 Iaccarino 2005; Schillinger 等，2012)和基于文献的 MPM 之间的比较。
SPH 可能是最受欢迎的无网格方法，因为它已被用于工程和计算机图形学，并被纳入许多商业软件，如 AUTODYN，PAM-CRASH，LS-dYNA 和 ABAQUS。尽管不是一种基于连续统的数值方法，但是处理大变形和断裂的能力使得 DEM 成为一种非常流行的技术，特别是在岩土工程领域(Cundall and Strack 1979; Scholtès and Donzé 2012; Sinaie et al。2018a)。我们可以将这两种方法的流行归因于它们在实现上的简单性，并且它们可以在2D 和3D 中工作。
超高速撞击问题的 SPH 和 MPM 的第一次比较研究是由 Ma 等人(2009a) ，Ma 等人(2009b)进行的。仿真结果表明，MPM 比 SPH (LS-DYNA 中的 SPH)更快、更精确。这可以用以下事实来解释: MPM 中的临界时间步长取决于背景网格的单元尺寸，而不是 SPH 中的粒子空间，因此 MPM 中使用的时间步长远远大于 SPH 中使用的时间步长。此外，没有相邻粒子的搜索，这是非常耗时的。然而，请注意，这个结论是相当肤浅的，因为它只是基于一些数值试验。可以肯定的是，一个工作的 MPM 模拟需要比 SPH 更少的棘手的自组织参数。典型的 MPM 模拟只需要网格间距和时间步长等数据，其方式与有限元法相同。
有趣的是，对于众所周知的泰勒冲击问题，显式MPM比LS-DYNA显式FEM更快，参见Ma等人（2009b）。Leroch等人也报道了这种优越的效率。（2018）用于微铣削模拟。原因是一样的：随着FEM网格变形，单元尺寸变得非常小，MPM中的临界时间步长比FEM中的大得多。
在隐式动力学的背景下，Iaconeta 等人(2017)对标准 MPM 和所谓的伽辽金无网格方法进行了比较研究（GMM）在Espluga（2014）中提出。该GMM与Li等人的OTM非常相似。（2010），只是使用了前向欧拉进行时间积分（图9）。作者发现，使用最大值和MLS的GMM对非常大的变形问题缺乏鲁棒性。这种鲁棒性的缺乏是由于这些形状函数的无网格性质，即必须确保在给定的材料点周围总是有足够的节点。Weißenfels和Wriggers（2018）提出了一种稳定的OTM。Wobbes等人（2020）对MPM和OTM进行了比较研究，他们发现尽管这两种方法非常不同，但它们有相似之处。
图9：最优无网格传输（OTM）方法：除了计算固体力学界新的最优传输理论外，OTM方法还采用了经典的离散化概念：节点和材料点。物质点起着积分点的作用，在节点的速度场中运动。由于没有网格，节点形状函数更为复杂。在MPM中没有从材质点到栅格的映射
批注7：最近的工作结合了MPM和SPH。例如，Raymond等人（2018）将SPH耦合到MPM。他们在主体中使用SPH，在表面周围使用MPM，以便于执行边界条件。He等人（2019）开发了所谓的SMPM，其中使用SPH函数对粒子速度和应力进行平滑。对于冲击模拟，SMPM的性能优于传统的MPM和SPH，并且比后者更有效。另见何和陈（2019），了解SMPM在应变定位中的应用。
对于颗粒流，Coetzee（2003），Coetze等人（2007），Ceccato等人（2018），Gracia等人（2019）报道了一些关于DEM和MPM性能的研究。基本上，如果采用适当的组成模型，MPM可以捕捉DEM的精度。有人可能会争辩说，MPM应该用于使用DEM模拟得出的本构模型的大规模模拟。在Dunatunga和Kamrin（2015）中，提出了一种颗粒介质的本构框架，该框架能够在一个环境中模拟跨越多个阶段的广泛颗粒行为：类固体静态行为、塑性流动（高达非常大的应变）以及材料的分离和再固结。蒋等人最近的工作(2019年)提出了一个组合 MPM-DEM 方法。
尽管MPM中最初使用的背景欧拉网格是为了处理碰撞，但最近使用有限元求积来评估弱形式积分的趋势（见刘和孙（2020））（这一趋势始于Sulsky和Gong（2016）与前面讨论的iMPM的工作）使MPM类似于浸入式边界方法，反之亦然。在浸入边界法中，所考虑的固体浸入比其大且不符合固体边界的网格中，见图10。一些流行的浸入式边界方法是有限元法（Parvizian等人，2007）和剪切有限元（Burman等人，2015）。这些方法中棘手的部分是由实体边界切割的单元的求积。最近，刘和孙（2020）使用物质点来表示实体域，水平集来捕捉边界，物质点的信息被映射到标准正交点进行数值积分。
图10：边界协调方法(左)与浸入式边界方法(右)
2.10 将MPM与其他方法耦合
由于没有一种方法是没有缺点的，因此很自然地希望将两种（甚至两种以上）方法结合在一起，以利用它们的最佳特性。在此，我们回顾了MPM与其他数值方法的耦合工作，无论是FEM还是分子动力学。
在变形适中的区域采用高效准确的有限元法，在变形较大的区域采用SPH或MPM等粒子法是合理的。张等人（2006）开发了一种显式材料点有限元方法（MPFEM），以利用有限元和MPM的优点。其基本思想是通过有限元网格离散实体，并在潜在的大变形区额外预定义计算网格。网格覆盖的节点被视为MPM粒子，其余节点被视为由FE节点。Lian等人后来对MPFEM进行了改进。（2011a），Lian等人（2012），Chen等人（2015）。这里没有介绍这些算法，我们详见张等（2016b）的教材。这些论文中报道的结果表明，耦合的MPFEM确实比纯MPM更快。MPM已与分子动力学相结合，以实现多尺度原子到连续体方法（Ayton等人，2001；Lu等人2006；Chen等人2012；刘等（2013）。正是MPM的粒子性质使耦合更加直接。由于我们缺乏这方面的经验，我们的讨论只是为了完整。我们请读者参考张等人的教科书。（2016b）了解详细信息。
