细说数值计算方法：拉格朗日法，欧拉法，混合法

目前，数值计算方法在科学研究和工程技术中得到了广泛应用，与理论研究和实验研究组成了现代科学技术的三大支柱，并具有快捷、安全和低成本的优势。然而，超高速碰撞、冲击侵彻、爆炸、金属冲压成形、裂纹动态扩展、流固耦合、应变局部化等强非线性问题涉及材料的特大变形、破碎甚至熔化和汽化，包含了一系列复杂的物理过程，给数值计算方法带来了巨大的挑战。

根据所采用的运动描述方法，数值计算方法可以划分为拉格朗日法、欧拉法和混合法三大类。

1.1拉格朗日法

在拉格朗日法中，计算网格固连在物体上，随其一起变形，因此材料与网格之间不存在相对运动(即迁移运动，也称对流运动)，控制方程中不存在对流项，大大地简化了控制方程及其求解过程。单元质量在计算过程中始终不变，但其体积由于网格变形而不断发生变化。拉格朗日法具有如下优点。

(1)拉格朗日法的概念清晰。由于不存在描述质量在网格间迁移的对流项，拉格朗日法的质量、动量和能量守恒方程比较简单，且在每个时间步中所需的计算量较小。

(2)网格线（面）和物体的外表面及材料界面在求解过程中始终重合，因此拉格朗日法易于处理边界条件和跟踪材料界面。

(3)在物体的运动过程中材料与网格始终重合，易于处理与变形历史相关的材料本构模型（如弹塑性及应变硬化和应变率效应）和失效模型。

图：拉格朗日网格

上图给出了一个典型的拉格朗日法的例子。如前所述，拉格朗日法的计算网格随物质一起变形扭曲。当计算网格严重扭曲甚至相互重叠时，将引起较大的误差，甚至使计算失败。这个问题可以通过使用滑移界面(slidelines)和重分网格(rezone)技术在一定程度上

得以克服。在物体之间可能产生相互滑移的问题中，需要使用滑移界面单元，如流体与固体壁之间的相互作用问题、子弹穿透靶体过程和碰撞物体之间的相互接触问题等。滑移界面在离散可能出现非常大的剪切变形和失效区域也是很有效的。大多数滑移界面元都是将加速度和速度分解为与界面相切的切向分量和与界面垂直的法向分量。当物体之间相互接触时，法向运动是连续的，而当物体之间脱离接触后，法向运动是完全独立的。当物体之间相互分离或接触面之间无摩擦时，切向运动是相互独立的，但当物体相互接触后并存在摩擦力时，切向运动需满足一定的关系。两个相互接触的物体在满足一定条件下可以相互分离,而已分离的两个物体可以再次相互碰撞。

在拉格朗日法中，当网格扭曲非常严重时将导致很大的误差，而且六面体（或四边形）网格扭曲后，其体积（或面积）可能成为负值，使得计算毫无意义。在显式积分方法中，为保证算法的稳定性，积分步长是由最小单元的特征尺寸控制的。当单元的扭曲、变形增加时，时间步长逐步减小，并趋于零，使得计算成本急剧增加，甚至难以完成。为了解决此问题，计算中当网格严重扭曲时必须重新划分网格，并将旧网格中的物理量映射到新网格中。网格重分技术已成功地应用于一维和二维问题，但对于三维问题，网格重分的过程既费时又复杂，甚至难以实现。对于记忆材料，由于新网格可能覆盖几个旧网格的区域，

而每个旧网格都有不同的历史，因此在变换内变量时也将遇到困难。在得到一个成功的新网格之前，一般需要多次网格重分过程，而且最终结果的精度在很大程度上取决于计算者的经验。

克服网格畸变的另外一种方法是采用侵蚀(erosion)算法，当单元的变形过大，其拉力或等效塑性应变达到一定阈值时就将该单元删除。侵蚀算法在侵彻问题的模拟中取得了一定的成功，但也带来新的问题，如单元侵蚀造成的质量不守恒给系统带来扰动，忽略了被侵蚀单元对系统的作用，进而造成计算误差：另外被侵蚀的单元也无法继续参与描述材料的变形，无法模拟超高速碰撞中碎片云的形成过程及其作用效果。

网格扭曲可以用多种方法来判断。最简单的方法之一是判断时间步长是否减小到了不能接受的地步，如果是，则认为网格已发生严重的扭曲。应在此时间步存储重启动文件，然后通过手工重新划分网格或通过专用程序自动进行网格重分。对三维问题，甚至最复杂和完善的网格重分程序也难以得到令人满意的结果。TOODY和DYNA2D等程序提供了功能强大的网格重分能力。

HEMP程序是美国劳伦斯利福摩尔国家实验室（Lawrence Livermore National Laboratory,LLNL)的Wilkins在20世纪50年代到80年代研发的拉格朗日型显式有限差分程序，早期被用于计算爆炸、冲击问题。EPIC6)(Elastic Plastic Impact Computations).是由Johonson等在20世纪70年代研制的拉格朗日型显式积分有限元程序，采用三角形单元和四面体单元，并引入了重分区(rezoning)和侵蚀(erosion)算法，可以计算高速碰撞和爆炸等问题。EPIC程序的发展略早于DYNA程序，它们的区别在于EPIC主要用三角形单元，而DYNA主要用四边形单元4.EPIC在研究领域一直发挥着作用，如美国陆军研究实验室用EPIC研究材料的强度和局部化行为。

PRONTO3D是由美国桑迪亚国家实验室(Sandia National Laboratory,SNL)开发的采用拉格朗日网格的三维瞬态动力学有限元法程序，用于求解高应变率下高度非线性材料的大变形问题。PRONTO3D采用八结点六面体单元和四结点四边形壳单元，采用变步长的显式时间积分，用自适应时间步控制算法来控制时间步长。程序中包含沙漏模态控制、生死单元技术、刚体处理、程序起爆等功能，用结点约束的方式处理接触。PRONTO3D将无网格光滑质点流体动力学法（smoothed particle hydrodynamics method,.SPH)作为一种特殊单元加入到程序中，并通过接触算法和有限元耦合。

DYNA系列的程序是最为知名的显式动力学有限元分析程序，现在使用的很多程序或软件(比如MSC.DYTRAN)都是在DYNA程序的基础上发展而来的。最初的DYNA3D程序源自20世纪70年代的美国LLNL国家实验室，其发展背景和PRONTO3D程序类似，但其商业化版本LS-DYNA取得了更为广泛的应用。一些商业有限元软件（如ANSYS)对LS-DYNA求解器进行了包装，提供了部分前后处理功能。Aup公司的OASYS也为LS-DYNA提供了前后处理环境。对早期各种程序更详细的总结和介绍可参见Zukas主编的《High Velocity Impact Dynaimcs》的第9章。

1.2欧拉法

对于可能产生严重网格扭曲的问题，或者是在相互分离的材料将发生混合的问题中，需要使用欧拉法。在欧拉法中，计算网格固定在空间中，不随物体运动，而材料相对于网格运动(如下图所示)，因此不存在网格畸变问题。各时刻的速度、压力、密度和温度等物理量是在空间点上计算的，而不是像拉格朗日法那样在物质点上计算，因此质量、动量和能量等物理量将跨越单元边界在单元间输运。各单元的体积在计算过程中保持不变，但单元质量由于输运而不断发生变化。欧拉方法比较适合模拟大变形问题，因此在计算流体力学中多采用欧拉方法，早期用于冲击爆炸问题的流体动力学程序也多采用欧拉方法。

图：欧拉网格

欧拉法只计算质量、动量和能量等物理量在单元间跨越网格边界的输运量，难以精确确定材料界面和自由表面的位置，因此欧拉法施加边界条件比拉格朗日法更困难，精度也更低。除非使用特殊的材料分界面定位技术，否则材料分界面将很快在计算网格中弥散。目前已提出了多种方法，如通过引入无质量的拉格朗日示踪点来计算材料分界面的实际位置。这些示踪点也可用于提供材料历史数据，但将大大地增加计算的复杂程度和计算量。

欧拉法的网格区域需要足够大，以在所有时刻均能完全覆盖物体，因此会存在大量的空单元（没有被任何材料占据的单元），而拉格朗日法只需对物体进行离散，不存在任何空单元。欧拉法有多种实现方法，比如在计算偏导数时直接引入材料输运项，也可以在每个时间步里先进行拉格朗日计算，再单独计算对流输运项。

HELP程序是典型的欧拉方法流体动力学程序，是20世纪60年代到70年代由Valsh和Hageman开发的有限差分法程序。为了处理多物质界面问题，HELP中借助无质量的拉格朗日示踪点来显示物质界面。CTH程序14是美国桑迪亚国家实验室在CSQ程序的基础上开发的一个软件体系，采用欧拉格式下的有限体积法，主要用来处理多维、多材料、大变形、强冲击波动的物理问题。在CTH中求解分为两步，第一步是拉格朗日步，网格随物质运动变形。第二步是映射步，变形的网格被映射回欧拉网格。CTH程序近年来仍在不断更新和维护，并发挥重要的作用。Sche田er等应用CTH对斜侵彻问题进行分析,Huerta等将CTH用于聚能装药侵彻的模拟分析。针对结构网格的自适应加密网格(adaptive mesh refinement,AMR)技术是Berger等在20世纪80年代提出的),AMR在发展过程中逐渐得到广泛关注并已形成了一些比较成熟的程序包。CTH程序已经引入了AMR技术，并应用于气体爆炸后冲击波对结构的作用。以及超高速碰撞问题的模拟。美国桑迪亚国家实验室把CTH和PRONTO3D整合到一起，形成了Zapotec程序，用于侵彻和爆炸问题。Bessette等对这套程序和方法进行了评估，并与其他几种类似的程序进行了对比，为美国即将发射的火星科学实验室的安全分析做准备。

北京理工大学宁建国研究组开发了使用欧拉网格的有限差分程序MMIC2D和MMIC3D,引入了界面处理技术，并用于聚能射流等问题的模拟，同时，为了配合求解程序的使用，也开发了前处理程序MESH2D以及可视化程序VISC2D和VISC3D.北京计算数学与应用物理研究所进行了大量流体动力学程序的研究开发工作，开发了欧拉网格的有限差分程序MEPH2D和MEPH3D29,3o,并引入自适应加密网格技术和界面处理技术，用于高速碰撞、炸药爆轰、聚能射流等问题的模拟。中科院力学所张德良和北京大学刘凯欣等基于高精度时空守恒元解元算法开发了SUPER CE/SE程序用于高速侵彻和爆轰问题的模拟。杨秀敏院士课题组基于高精度、高分辨率差分格式WENO和贴体坐标编制了空气冲击波流场计算程序EF3D,可以研究不同爆炸方式、不同入射反射背景条件下的各种空气波流场参数变化规律。

1.3混合方法

拉格朗日法和欧拉法都存在严重的缺陷，但也都具有各自的优势，如果能将二者有机地结合，充分吸收各自的优势，克服各自的缺点，则可解决一大批只用拉格朗日法或欧拉法所解决不了的问题。

1.3.1任意拉格朗日-欧拉法

任意拉格朗日-欧拉法(arbitrary Lagrangian-Eulerian method,ALE)最早是由Noh(1964年)以耦合欧拉-拉格朗日法的术语提出的，并用有限差分法求解带有移动边界的二维流体动力学问题。在Noh的研究工作中，网格点可以随物质点一起运动，但也可以在空间中固定不动，甚至网格点可以在一个方向上固定，而在另一个方向上随物质点一起运动，因此ALE描述也被称为耦合欧拉-拉格朗日描述。例如在液体表面波的传播问题中，网格点在垂向随物质点一起运动，而在水平方向上固定不动，如图所示。这样很容易描述液体表面的运动，而且网格不会发生扭曲。

图：流体表面波传播问题

在ALE法中，计算网格可以在空间中以任意的形式运动，即可以独立于物质坐标系和空间坐标系运动。这样通过规定合适的网格运动形式可以准确地描述物体的移动界面，并维持单元的合理形状。类似于欧拉法，在ALE法的控制方程中也将出现对流项，因此使得系数矩阵不对称，还可能得到振荡解，需要进行相应的数值处理，如使用迎流权函数或彼得洛夫-伽辽金法(Petrov-Galerkin)法。其实质是在原微分方程中增加一扩散平衡项。纯拉格朗日法和欧拉法实际上是ALE法的两个特例，即当网格点的运动速度等于物质点的运动速度时就退化为拉格朗日法，而当网格固定于空间不动时就退化为欧拉法。

ALE法最早是为了解决流体动力学问题而引入的，并且使用有限差分法。由于核反应堆结构安全分析的需要，ALE法被引入有限元法中，用以求解流体与结构相互作用问题，目前已被应用于固体力学领域中求解大变形问题，如接触碰撞、弹性断裂力学和加工成形等。

1.3.2质点网格法(PIC)

另外一个典型的混合方法是由美国洛斯阿拉莫斯国家实验室（Los Alamos National Laboratory,LANL)的Harlow及其负责的计算流体动力学小组(T-3)于1955年提出的质点网格法(particle in cell,PIC)。PIC方法采用双重描述，即将流体离散为质点（拉格朗日描述），而计算网格仍然使用欧拉网格。每个质点具有一定的质量、速度和能量，计算过程中质点在欧拉网格之间迁移，描述了流体的运动。在每个时间步中，质点网格法的计算分为两步：第一步先忽略偏微分方程中的输运效应，进行拉格朗日计算，用差分法计算由压力分布所引起的欧拉网格上的速度（或动量）和能量的变化。若一个网格内含有多种流体，则按一定的规则将能量的变化量分配给不同流场的质点。第二步是在第一步的基础上，按一定的加权平均方法计算每个质点的速度和在本时间步结束时的新位置。第二步是质点迁移计算，其实质是对第一步计算中忽略的输运效应进行补偿。质点网格法的思想对以后许多计算方法的发展都产生了重要的影响，并成功地应用于等离子体模拟中。Julien Forest等基于Java开发了三维PIC等离子体模拟开源软件PicUp3D。

为了解决PIC法中的粒子脉动和存储量大的问题，在Harlow的指导下，Gentry、Martin和Daly于1966年进一步提出了流体网格法(fluid-in-cell method),简称FLIC法。它和PIC法一样采用欧拉网格，但在第二步中不计算质点的迁移，而是计算连续流体的迁移，即先算出通过网格边界的质量输送量，得出每个网格的新密度，再算出通过网格的质量所携带的动量和能量的输送量，最后得到每个网格的新速度和能量。FLIC法还有一套局部网格单元的计算格式，能计算一些边界形状比较复杂的问题。

为了模拟不可压缩自由表面流动问题，Harlow和Velch在PIC法的基础上进一步建立了标记网格法(marker and cell method,MAC)。MAC法仍然采用欧拉矩形网格单元，对纳维-斯托克斯方程则用差分近似，而把压力和速度分量作为基本未知量。此外，MAC法还在网格中布置适量的无质量标记点，并在整个计算中跟踪每个标记点，以判定网格里有哪种流体存在。MAC法可用于计算多种流体和带有自由面的问题，是模拟不可压缩流动问题的一个成功的数值方法。

经典的PIC法不是一种完全的拉格朗日质点方法，它只在质点上存储质量和位置信息，而其他物理量仍然存储在网格上。动量在网格和质点之间的迁移将产生较大的数值耗散，严重损害了PIC的计算精度。有两条途径可以提高PIC的精度：采用更精细的动量迁移算法，或者令质点携带所有物质信息。Brackbill等将流体的所有物质信息均赋予质点，基于有限差分法发展了一种用于求解流体流动问题的低数值耗散的PIC方法——FLIP(fluid-implicit-PIC)方法。FLIP是一种完全的拉格朗日质点格式。

1.3.2物质点法(MPM)

FLP主要用于分析本构方程与变形历史无关的材料（如流体），因此本构方程是在网格结点处计算的。Ssky等将本构方程改为在各质点上计算以便于处理与变形历史相关的材料，并通过等效积分弱形式，采用质点离散建立动量方程的离散格式，将FLP扩展到固体力学问题中，并将其命名为物质点法(material point method,MPM)。

物质点法仍然采用拉格朗日质点和欧拉网格双重描述，它将连续体离散成一组质点，如图所示。每个质点代表一块材料区域并携带该材料区域的所有物质信息，如质量、速度、应力和应变等，因此所有质点的集合代表了整个材料区域。计算网格仅用于动量方程的求解和空间导数的计算，它不携带任何物质信息。计算网格可以在空间中固定，也可以按某种方式自由布置。在每一个时间步中，质点和计算网格完全固连（拉格朗日计算步），计算网格提供了对物体的拉格朗日有限元离散，因此可以用标准的有限元法在计算网格上求解物体的运动方程。计算网格结点的运动方程可以通过将质点的运动量映射到计算网格而得到，求解后将计算网格结点的运动量映射回各质点，从而得到这些质点在下一个时刻的运动量（对流计算步）。

图：物质点法示意图(a)材料区域；(b)质点代表区域；(c)物质点法离散

物质点法是一种完全的拉格朗日质点类方法，在每步中质点和计算网格没有相对运动，避免了欧拉法因非线性对流项产生的数值困难，并且容易跟踪物质界面。质点已经携带了连续体的所有物质信息，因此物质点法在每个时间步结束时抛弃变形后的计算网格，在新的时间步中仍可以采用未变形的计算网格，从而避免了拉格朗日法因网格畸变而产生的数值困难。物质点法发挥了拉格朗日方法和欧拉方法各自的长处，克服了其弱点，在冲击、接触、流固耦合等涉及大变形和材料破坏的问题中具有明显的优势。马上和张雄等从形函数特性、影响点搜索、接触算法、稳定性和边界条件等方面详细比较了MPM与SPH法。研究结果表明，与SPH法相比物质点法具有更高的计算效率和稳定性。

图：MPM与SPH计算效率比较

上图比较了LS-DYNA的SPH模块和张雄课题组研发的冲击爆炸三维物质点法数值仿真软件MPM3D在计算方块平动时每个时间步所需的CPU时间和质点总数之间的关系，可见LS-DYNA的SPH模块的计算量随着质点总数的增加速率远高于MPM3D。在32位个人计算机上，LS-DYNA的SPH模块当质点总数超过106万时崩溃，而MPM3D单机可以求解的质点总数可达350多万。比较了由LS-DYNA的并SPH模块和MPM3D得到的Taylor杆碰撞问题的最终构形俯视图，SPH结果中出现了多处数值断裂和不稳定现象。

图：MPM与SPH求解Taylor杆冲击的比较

物质点法最初采用显式积分，适合于分析冲击爆炸等载荷作用时间短的瞬态问题。针对准静态和结构动力学问题，Cummins等和Sulsky等给出了物质点法的隐式积分算法，采用迭代方法求解动量方程，无需显式组装切线刚度矩阵。Guilkey等则采用直接法求解动量方程组，提出了需显式组装切线刚度矩阵的隐式积分算法。

在物质点法中，动量方程在计算网格结点上求解，因此需要建立计算网格单元的质量阵，既可采用集中质量阵，也可采用协调质量阵。Tan等在计算裂纹尖端的能量释放率时指出协调质量阵能够获得更好的能量守恒性质：Lov等分析了质量阵形式对能量和动量守恒性的影响，结果表明协调质量阵有较好的能量守恒性。

Bardenhagen等发现，物质点法中形函数导数在计算网格结点处不连续，使得质点穿越计算网格时产生较强的数值噪声，导致应力振荡。为此，他们采用彼得洛夫-伽辽金法对MPM进行了扩展，提出了广义插值物质点法（generalized interpolation material point method,GMP),建立了一族具有C1连续性的形函数，有效减弱了数值噪声。在GMP中，质点的影响范围由质点的特征函数表示。在计算过程中，需要实时计算该质点的形状以判断该质点的影响范围。为了避免实时计算质点形状的困难，一般假定质点的形状不变，通过其现时位置和初始形状大小判断其影响范围，该方法称为uGIMP。Zhang和Ma等通过将标准物质点法的形函数导数值与FLIP中基于背景网格结点的形函数导数值加权平均，建立了一种新的形函数导数计算方式，以降低质点穿越网格引起的噪声。该方法不需要跟踪计算质点的形状，并且质点仅影响其所在的单元，因此便于非结构化计算网格的计算。

在物质点法中，质点和计算网格之间采用单值映射函数，因此自然地满足了物体间的非穿透条件。对于无滑移粘着接触，MPM不需采取额外的处理；对于需要考虑滑移接触的问题，则需要建立相应的接触算法。York等首先提出了一种简单的接触算法：在物体发生接触分离前按正常单值速度场计算，当检测到物体发生接触分离时采用物体各自独立的速度场进行计算，从而实现了接触物体正常的分离。该算法过于简单，并且未考虑接触物体之间的摩擦作用。Bardenhagen等借鉴有限元法接触算法的思想，建立了考虑摩擦的物质点接触算法，实现了物体间的接触、滑移和分离过程的模拟，随后又对接触判断条件进行了改进。黄鹏和张雄等指出Bardenhagen等提出的接触算法中法向量计算不满足共线条件易导致动量不守恒的事实，为此提出了一种接触界面法向量计算方法，保证了接触面上法向量共线，并将其用于中低速冲击侵彻问题。Chen等则从多重计算网格的思路建立了物质点法接触算法，通过公共计算网格更新接触面处法向的速度，在物体各自的计算网格计算其他物理量，但该算法对内存需求量大。马志涛和张雄等在此基础上建立了局部多重网格接触算法，大幅降低了对内存的使用量，提高了接触算法的效率。由于物质点接触算法在计算网格结点上实施，马志涛和张雄等研究表明仅利用非穿透条件判断接触易导致接触的提前发生，为此提出了一种计算物体之间相邻最近点距离的公式，改进了物质点法的接触判断方法。

物质点法采用质点离散，各质点之间无拓扑关系，因此易开展h自适应分析。Tan等建立了质点和计算网格的层级结构化自适应分析算法，通过局部加密质点和计算网格计算了裂纹尖端的能量释放率。针对冲击爆炸中易出现的非物理现象的数值断裂问题，马上和张雄等提出了质点自适应分裂的物质点法，有效避免了数值断裂问题。

在并行算法方面，Parker等基于消息传递模型(MP)研究了MPM的并行算法，开发了Uintah程序。黄鹏和张雄等提出基于信息共享模型(OpenMP)的MPM并行算法，张衍涛和张雄等进一步建立了基于OpenMP的区域信息交错更新的MPM并行算法，精巧地避免了数据竞争问题。Ma和Lu等基于并行程序开发平台SAMRAI2建立了GIMP的并行算法。至此，已经建立了基于两种并行机制的MPM并行算法。

在物质点法中，本构方程在质点上计算，动量方程则在背景网格上计算，因此既可在动量方程更新前也可在动量方程更新后进行本构计算。Bardenhagen分析了两种次序对能量守恒性的影响，结果表明前者能量守恒性优于后者。Wallstedt和Guilkey具体分析了质点速度到计算网格映射过程中的精度问题，并提出了一种改进方案。Steffen等对MPM进行了误差分析，指出质点跨越网格时导致物理量间断引起的误差是时间步长的二阶小量。

经过不断发展和开发，物质点法已经在超高速碰撞、冲击侵彻、爆炸、裂纹扩展、材料破坏、颗粒材料流动和岩土冲击失效等一系列涉及材料特大变形的问题中体现了优势。

在超高速碰撞问题方面，马上和张雄等以空间碎片防护问题为背景采用MPM研究了弹丸碰撞薄板和厚板等一系列问题，黄鹏和张雄等则进一步采用并行的MPM方法精细模拟了超高速碰撞问题中的碎片云问题。宫伟伟和张雄等基于泡沫铝的Micro-CT扫描图片重构了具有真实微观结构的三维泡沫铝物质点模型，并模拟了两种含有泡沫铝的Whipple防护结构的超高速碰撞问题，计算结果与实验结果吻合。

在冲击侵彻问题方面，Sulsky等模拟了Taylor杆问题、金属成形问题和钢球侵彻铝靶体问题；黄鹏和张雄等进一步采用物质点接触算法模拟了中低速冲击侵彻问题，所得计算结果与实验结果吻合；马志涛和张雄等基于局部多重网格物质点接触算法模拟了同样的弹体侵彻问题；宫伟伟和张雄等将Deshpande-Fleck泡沫模型引入到物质点法中，对泡沫铝材料动态压缩性能进行了三维仿真模拟：王晓军和张雄等将JH2、HJC和RHT本构模型引入到物质点法，模拟了陶瓷和混凝土的冲击问题。

在爆炸问题方面，Hu和Chen83基于物质点法模拟了爆炸对混凝土墙的破坏作用：马上和张雄等将MPM应用于高能炸药爆轰问题；王宇新等采用MPM系统研究了爆炸焊接问题；张忠等基于MPM分析了非均质固体炸药的爆轰过程及其对金属材料的破坏作用。

在裂纹扩展问题方面，Nairn的研究组采用MPM做了大量工作，模拟了二维、三维的裂纹扩展问题。采用多速度场描述裂纹面以避免裂纹面之间的穿透，在物质点法中计算了裂纹尖端各种参数，如动态J积分、I、Ⅱ和I型裂纹的应力强度因子,模拟了脆性材料的断裂问题并计算了材料断裂过程中的能量释放率。此外，Gilabert等采用MPM模拟了陶瓷材料中裂纹萌生和扩展过程；Daphalapurkar等在GMP中实现了内聚力模型(cohesive zone model)模拟了延性和脆性材料中IⅡ型裂纹动态扩展的问题；Wang等采用非规则的计算网格计算了二维混合型裂纹扩展问题。

在模拟材料失效问题方面，Chen的研究组做了大量工作，采用MPM模拟了冲击载荷下脆性材料的动态失效问题、受压缩的薄膜脱层问题和材料在局部加热情况下的失效问题;Shen采用MPM模拟了玻璃在冲击作用下的失效破碎问题；Schreyer和Sulsky等在MPM中采用脱聚本构模型(decohesion constitutive model)研究了复合材料中的脱层失效模式；Ionescu等采用MPM模拟了各项异性的软组织失效问题；Li和Pan等基于MPM研究了脆性材料在冲击作用下的失效问题。

在岩土冲击动力学问题方面，Coetzee等采用二维物质点法模拟了岩土中错墩受力问题；Andersen等采用GIMP分析了边坡失效问题；黄鹏采用MPM模拟了半球壳体侵彻岩土和边坡失效问题。

在颗粒物流动问题方面，Wieckowski等采用MPM模拟了储料垛内的颗粒流动问题；Bardenhagen等采用MPM模拟了颗粒物受剪作用下的流动过程；Coetzee等基于MPM模拟了挖斗铲玉米颗粒堆的过程；王津龙等基于MPM模拟了柱状堆石成形过程。

此外，MPM在海冰动力学问题、切削问题、碳纳米管增强型复合材料、薄膜流固耦合问题和多孔介质的动力学响应分析等问题中也有应用。

为综合物质点法和其他方法的优势，一些学者将物质点法与其他方法耦合。Guilkey建立了多物质欧拉方法和MPM的耦合方法，分别采用多物质欧拉算法模拟流体区域、物质点法模拟固体结构。Gilmanov和Acharya将MPM与浸没边界法进行了杂交，其中MPM用于求解固体结构。Guo和Yang采用握手区方法耦合了MD和MPM,模拟了Cu-Cu和Si-Si的群簇高能碰撞问题。Lu等通过逐步细化MPM计算网格间距过渡到分子尺度，在过渡区建立质点和分子的一一对应关系，将MPM与MD进行无缝耦合，分析了硅在纳米尺度下的拉伸过程。Chen等采用耦合的MD/MPM进行了纳米尺度的切削和薄膜生产过程的多尺度分析。

为了发挥物质点法和有限元法各自的优势，张雄等提出了显式物质点有限元法。该算法初始采用有限元离散整个求解区域，同时在材料可能发生大变形的空间区域布置计算网格。被网格覆盖的结点被看成质点，其动量方程在预先设置的计算网格上求解，而在其余结点仍然是有限元结点，其动量方程在原有的有限元网格上求解。如果不设置计算网格，显式物质点有限元法退化为显式有限元法：如果计算网格在所有的时间步中都能覆盖所有的结点，则显式物质点有限元法退化为物质点法。廉艳平和张雄等将有限元法的杆单元引入到物质点法中，提出了杂交物质点有限元法，用以研究钢筋混凝土的冲击侵彻问题。杂交物质点有限元法采用质点离散混凝土，用杆单元离散钢筋，避免了在钢筋直径方向上的离散，降低了MPM的离散规模。在冲击侵彻问题中，大部分材料区域始终处于小变形阶段，而物质点法在计算材料小变形时的精度和效率均低于有限元法，因此廉艳平和张雄等提出了耦合物质点有限元法，采用有限元离散仅涉及材料小变形的物体，而用物质点法离散涉及材料特大变形的物体，通过接触算法实现两个离散体之间的相互作用。为了充分发挥有限元法计算材料小变形时效率和精度的优势以及物质点法模拟材料特大变形的优势，廉艳平和张雄等进一步提出了自适应物质点有限元法，初始采用有限元离散所有物体，在计算过程中将发生畸变或破坏的单元自动转化为质点，通过计算网格实现同一个物体内不同离散区域之间的相互作用，采用耦合物质点有限元法计算各物体之间的接触问题，实现了有限元法到物质点法求解的自适应转化。数值结果表明，在计算精度一致的情况下，自适应物质点有限元法的计算效率显著高于物质点法。

随着MPM算法上的不断完善和扩展，目前出现了一批以MPM为核心算法的计算软件，包括UCF、FLIP-MPM-MFM、nairn-mpm-fea、MPM3D和MPMsim等。美国犹他大学C-SAFE(Center for the Simulation of Accidental Fires and Explosions)研究中心l31在美国能源部“加速战略计算计划”（Accelerated Strategic Computing Initiative,ASCI)的支持下开发了UCF(Uintah Computational Framework)大规模并行可视化数值模拟工具，物质点法作为其中的核心算法之一，主要进行固体和含能易爆材料的爆炸、燃烧、相变等极端载荷环境的数值模拟。美国Los Alamos国家实验室的计算流体动力学和固体力学小组开发的CartaBlanca软件将物质点法嵌入到多相流框架(multiphase flow framework)中，开发了FLIP-MPM-MFM(FLuid Implicit Particle-Material Point Method-Multiphase Flow Method)模块,用来求解复杂力学问题，包括固体的失效和侵彻、热传导、相变、炸药化学反应和多相流等，能模拟复杂武器侵彻和弹药爆炸相互作用。Nairn研究组采用C++开发了计算力学软件nairn-mpm-fa),其主要包括有限元分析模块和物质点计算模块，物质点模块部分主要求解二维、三维复合材料的失效、裂纹扩展和面板之间的脱胶等动力学问题。清华大学张雄课题组采用C++开发了三维显式物质点法程序MPM3DsO,目前已实现了USF、USL和MUSL

求解格式，GMP算法，接触算法，自适应算法，有限元法，杂交物质点有限元法，耦合物质点有限元法和自适应物质点有限元法，具有多种材料模型、状态方程和失效模型，并分别利用OpenMP和MPI实现了SMP和MPP两种体系下的并行化。2009年，美国的ADI公司发布了ASTE-P1.0软件，将物质点法加入到了FSI、CSD、CFD模块，用来计算气动弹性问题。2010年，英国的Jveer Aerospace Ltd公司推出了商业化的物质点法计算软件MPMsim。

注：本文摘自清华大学张雄教授的专著——《物质点法》，如有侵权，请联系删除。

如果你觉得此文对你有帮助，请点赞，谢谢！

计算机技术在工程学科中得到了广泛的应用，力学方面，计算机技术成为了力学的第三次革命性技术，现在基于计算机的数据科学已经逐步成为力学等其他科学发现的第四范式。人工智能、大数据、数字孪生等概念已经逐步成为当今时代的主题。智能制造、智能算法、数据驱动力学已经在力学当中展现出巨大潜力，吸引了大量的研究人员。同时高性能显卡和多核中央处理器的出现为大规模数值模型的高性能计算提供了强大算力。然而因为计算力学的论文较多，设计内容较广，需要的知识量较大，不仅需要力学，数学，物理的知识，还需要计算机的知识。入门门槛较高，因此我建立了此微信公众号（STEM与计算机方法），希望通过自己的学习加上文献翻译和整理，帮助新手快速掌握前沿研究的热点和聚焦，轻松入门计算的相关研究（实验、理论、数值计算方法），从而吸引和聚焦更多对该技术和研究领域感兴趣的华人朋友，为推动计算力学的科学研究的发展和交流做一点儿贡献！