极端变形问题的物质点法研究进展

极端变形问题的物质点法研究进展
张雄;刘岩;张帆;陈镇鹏
【摘要】By employing both the Lagrangian and Eulerian descriptions,the material point method (MPM) discretizes a material domain into a set of particles (Lagrangian description) moving in a grid fixed in the space (Eulerian description).Thus,the MPM combines the advantages of both the Lagran-gian methods and Eulerian methods,which makes MPM especially effective in modelling problems with extreme deformation such as impact and explosion problems.The research achievements of our group on MPM algorithms for impact,explosion,fluid-structure interaction and other problems with extreme deformation are systematically summarized in this paper.The three-dimensional parallel explicit simulation software MPM3D developed by our group is briefly introduced.The applications of MPM3D in hyper-velocity impact,perforation,explosion,slope failure,metal cutting,and fluid-structure interaction problems are described as well.%物质点法采用了拉格朗日和欧拉双重描述,将物体离散为一组在固定于空间的网格(欧拉描述)中运动的质点(拉格朗日描述),有效地综合了拉格朗日法和欧拉法的优点,是分析冲击爆炸等极端变形问题的一种有效方法. 本文系统总结了本课题组近年来针对冲击、爆炸和流固耦合等极端变形问题在物质点法算法方面的研究进展,并简要介绍了开发的三维显式并行物质点法数值仿真软件MPM3D及其在超高速碰撞、侵彻、爆炸、边坡失效、金属切削和流固耦合等问题中的应用.
【期刊名称】《计算力学学报》
【年(卷),期】2017(034)001
【总页数】16页(P1-16)
【关键词】物质点法;冲击;爆炸;超高速碰撞;流固耦合
【作者】张雄;刘岩;张帆;陈镇鹏
【作者单位】清华大学航天航空学院,北京 100084;清华大学航天航空学院,北京100084;清华大学航天航空学院,北京 100084;清华大学航天航空学院,北京100084
【正文语种】中文
【中图分类】O351
数值计算在科学研究和工程技术中得到了广泛应用，与理论和实验一起成为现代科学技术的三大支柱[1]，并具有快捷、安全和低成本的优势。

然而，极端变形问题(如超高速碰撞、冲击侵彻、爆炸、金属加工成型、应变局部化、边坡失效和液体
晃动等)是几何、材料和边界条件均为非线性的多物理场强耦合问题，涉及高速、
高压、高温、相变和化学反应，气体、液体和固体等多种物质间相互耦合甚至混合，材料发生严重扭曲、破碎、熔化甚至气化，给数值计算带来了巨大的挑战。

拉格朗日法的网格与材料固连，因此格式简单，效率高，可以准确跟踪材料界面，易于引入与变形历史相关的材料模型。

但是，在模拟极端变形问题时，材料的极端变形将会导致拉格朗日法的网格畸变，严重降低结果精度，甚至会使单元雅可比行列式在高斯点处为0或负值，使计算异常中止。

网格畸变还会使显式时间积分的
临界时间步长降低1～2个数量级，显著增加计算量。

另外，拉格朗日网格类方法也难以有效地模拟材料的破碎、熔化和气化等复杂行为。

LS-DYNA等商用软件引
入侵蚀算法，将发生畸变或者破坏的单元人为地删除，以克服网格畸变带来的困难。

但是侵蚀算法是一种非物理方法，会造成质量、动量和能量损失，且结果在很大程度上取决于侵蚀阈值的选取。

欧拉法的网格固定在空间中，因此不存在网格畸变困难，适合于分析极端变形和流体流动问题。

但由于网格与材料独立，不易准确跟踪材料界面，难以引入与变形历史相关的材料模型，且非线性对流项也会导致数值求解困难。

由质点网格法PIC(Particle in Cell)[2,3]发展起来的物质点法MPM(Material Point Method)[4-7]采用了拉格朗日和欧拉双重描述，将物体离散为一组在固定
于空间的网格(欧拉描述)中运动的质点(拉格朗日描述)。

质点携带了所有物质信息，因此便于跟踪材料界面和引入与变形历史相关的材料模型。

固定于空间的网格只用于求解动量方程，不携带任何物质信息。

物质点法是一种基于有限元法的拉格朗日粒子类方法，有效地综合了拉格朗日法和欧拉法的优点，是分析冲击爆炸等极端变形问题的一种有效方法。

本课题组自2004年开始，致力于冲击爆炸等极端变形问题的数值模拟方法和软件研究，先后建立了物质点法的高效实现方案、改进的物质点接触算法、自适应物质点法、并行物质点算法、不可压物质点法、物质点有限元法、杂交物质点有限元法、耦合物质点有限元法、自适应物质点有限元法以及物质点有限差分法等，并在算法研究的基础上，基于C++，Qt，VTK，CMake，OpenMP和MPI开发了可运行于Windows，Linux和Mac OS等多种平台上的三维显式并行物质点法数值仿真软件MPM3D，可应用于超高速碰撞、侵彻、爆炸、边坡失效、金属切削和流固
耦合等问题中。

本文系统总结了本课题组近年来针对极端变形问题在算法、软件和应用方面的研究工作和成果。

基于更新拉格朗日格式的控制方程的弱形式为
式中下标i和j指示空间坐标分量，并遵守Einstein求和约定，为当前时刻的密度，
bi为体力，ui为位移，δui为虚位移，σi j为柯西应力张量，Γt为现时构形中的指定面力边界，i为指定的面力。

2.1 物质点法的基本格式
物质点法将材料区域离散为一组相对背景网格运动的质点，如图1所示。

每个质
点表示一个材料团簇，并携带其所有物质信息，如质量、速度、应力和应变等。

背景网格仅用于求解动量方程和物理量的空间导数，一般取为固定于空间的规则网格。

将材料区域离散为质点后，其密度可近似为
式中np为物质点总数，mp为物质点所代表区域的质量，δ为Dirac函数，xp为质点p的坐标。

在求解动量方程时，背景网格和质点完全固连，可以将背景网格看成是对材料区域的有限元离散，因此域内任一点x的位移可以近似为
式中ng为背景网格结点总数。

将式(2，3)代入式(1)，并采用集中质量阵，可得
式中vi I为结点I的速度，
是结点I的质量，其中，NI p=NI(xp)，
是结点力，其中
是结点内力，
是结点外力，σi j p和bi p分别为质点p的应力和体积力。

为了公式简洁，这里
没有考虑自然边界条件。

背景网格结点动量方程式(4)可以采用显式积分或者隐式积分求解，参考文献[6，7]。

由以上格式可以看出，物质点法在每个时间步中可以看成是采用了质点积分(即将
质点取为积分点)的拉格朗日有限元法。

当材料区域中的某个单元内不存在任何质
点(即积分点)时，该单元为空单元，即材料在该单元处发生了数值断裂。

在超高速
碰撞等问题中，可以利用这个特性对材料破坏进行粗略模拟，但对某些问题可能会产生很大的误差。

例如在模拟聚能射流问题时，在射流形成过程中材料会产生很大的拉应变，质点之间的距离不断增加，最终产生空网格，造成数值断裂。

为此，马上等[8]提出了质点自适应分裂技术，当质点在某个方向的累积应变达到给定阈值时，将其在该方向上分裂为两个质点，有效地避免了数值断裂问题。

2.2 接触算法
在MPM中，质点和背景网格之间采用单值映射函数，自动满足无穿透条件，可
以处理无滑移接触问题。

但若要使界面间可以相互滑移和分离，则需要引入接触算法。

York等[9]提出，在物体相互接近时，采用标准物质点格式，当物体有相互分离趋势时，各物体采用独立的速度场进行计算，从而实现了接触物体的正常分离。

Bardenhagen等[10]基于拉格朗日乘子法建立了物质点接触算法，并考虑了摩擦。

Hu等[11]从多重背景网格的思路建立了物质点法接触算法，接触点的法向运动在
公共背景网格上计算，切向运动则在物体各自的背景网格上计算，保证了接触点处法向运动连续，又允许切向运动独立。

物质点法的接触算法可分为三个步骤，(1) 进行接触检测，确定可能发生接触的位置；(2) 计算接触点处接触面法线向量，以确定接触力的方向；(3) 计算接触力大小，并通过背景网格施加到相互接触的物体上。

物质点法中的接触关系是通过背景网格定义的。

接触检测时，需要找出所有发生接触的背景网格结点。

若不同物体的质点对同一个背景网格结点有贡献，则认为两个物体可能发生接触，如图2所示。

物质点法的接触算法一般通过背景网格结点的质量梯度来计算物体r的表面法向量，即物体r在背景网格结点I处的表面法向量为
由式(9)得到的接触面法向量不共线，因此不满足动量守恒条件。

黄鹏等[12，13]
提出了接触面法向量计算的修正公式，保证接触面上法向量共线，解决了原接触算法中动量不守恒的问题。

在施加接触力时，首先假设没有发生接触，不同物体的动量方程在各自独立的背景网格上独立求解，得到各自的试探速度。

若试探速度
表明两个物体将在该结点处发生穿透，需要施加法向接触力
fc,nor,nI=-
对试探解进行修正。

马志涛等[14]指出式(10)会导致物体提前发生接触。

当物体之间的距离小于两个背景网格单元尺寸时，虽然两个物体的材料区域并没有发生重叠，但按式(10)则需要施加接触力。

为此，提出了一种辅助接触判断条件
式中为两个物体表面与接触结点I最近的距离在该结点法线方向上的投影之和，如图3所示，dc为背景网格的结点间距，λ是人为选取的系数，可取为0.5。

式(12)表明，当两个物体的间距时，才需要施加接触力。

上述接触算法需要建立多重背景网格用于独立求解各个物体的动量方程，即每个物体需要一套背景网格，大大地增加了内存量。

马志涛等[14]提出了局部多重背景网格接触算法，仅在发生接触的区域为每个物体建立单独的背景网格，其余区域仍然共用一套背景网格，如图4所示。

该方法大大地降低了内存使用量，提高了接触算法的效率。

2.3 网格自适应技术
背景网格的尺寸是影响物质点法计算结果精度的主要因素。

加密网格可以提高结果的精度，但同时会大幅度地增加计算量。

为了在提高精度的同时尽可能减少物质点法的计算量，本课题组先后建立了动态网格、移动网格、多重网格、多级网格和自适应网格等技术。

2.3.1 动态网格
与欧拉法类似，物质点法的背景网格需要覆盖物体可能运动到的所有区域。

在大多数实际问题中，质点一般不会充满整个背景网格，即存在大量的不含有任何质点的空单元，采用静态网格(图5(a))时将会浪费大量的计算机内存和计算时间。

如在模拟某斜侵彻问题时，在每个时间步中大约有87%的空单元，即只有13%左右的单元中存在质点。

马志涛等[14]基于C++提出了动态背景网格算法，在初始时刻只
创建结点和背景网格单元的指针数组，不创建其对象，因此不占内存，也不参与计算。

在计算过程中，只有当某个单元中存在质点时，才创建该单元及其相关结点的实例，同时删除空单元及其相关结点的实例(图5(b))。

动态背景网格技术可大幅度降低内存耗费量并提升计算效率。

2.3.2 移动网格
在计算过程中，某些质点可能会飞出预设的背景网格区域。

为了避免这一现象，需要在初始时刻设置足够大的背景网格，但会大幅增加内存使用量和计算量。

马志涛等[14]提出了移动网格技术，在计算过程中根据质点的当前位置自动调整背景网格在各方向上的单元尺寸，以在不改变网格拓扑的情况下覆盖整个材料域(图6(b))。

保持网格拓扑结构不变而改变单元尺寸可能会导致单元在某方向的尺寸过大，如图6(b)所示，影响计算精度。

杨鹏飞[15]对移动网格技术进行了改进，通过调整网格的拓扑结构使得单元尺寸保持不变，同时删除没有包含质点的背景网格单元(图
6(c))。

2.3.3 多重移动网格
在用物质点法模拟多个物体时，即使物体之间只在局部发生作用，也需要将所有物体的质点信息映射到同一套背景网格上求解。

杨鹏飞[15]提出了多重移动网格技术，允许每个物体具有各自的背景网格，根据物体的运动范围进行实时调整，并且可以关闭或者激活个别网格。

如在模拟多层侵彻问题时，在弹体与靶体发生接触前，动
量方程在各自的背景网格中求解(图7(a))，接触后则在同一套网格中求解(图7(b))，而在弹体穿透靶体并与靶体分离后,又可以在各自的网格中求解(图7(c))。

在弹体与靶体发生接触前，靶体没有外载荷和初速度，因此可以关闭靶体的背景网格，以进一步提高计算效率；在弹体和靶体分离后，当靶体的动能不再发生变化时也可以关闭其背景网格(图7(d))。

2.3.4 多级网格及其自适应
为了提高局部化问题的求解精度和效率，杨鹏飞[15]提出了多级网格法，允许在同一套求解体系下存在不同级别(密度)的背景网格，即在局部大变形区域采用较密的物质点和背景网格，而在其余区域采用较疏的物质点和背景网格，如图8所示。

多级背景网格需要用户事先指定要加密的区域和加密的级数，并给定相应的多级背景网格。

杨鹏飞[15]进一步提出了自适应背景网格技术，根据各时间步的计算结果(如应变梯度或能量范数)，在梯度较大的区域自动加密背景网格并分裂物质点，提高了求解精度和效率。

如图9所示,当单元2满足分裂条件时，自动将其分裂为4
个单元，并且单元可以进行多级分裂。

2.4 不可压物质点法
自由液面的流动问题在工程实际中有着大量的应用，如航天器在变轨、交会以及对接过程中液体推进剂在外界刺激下产生的剧烈晃动和核反应堆堆芯在地震及海啸等突发激励下的晃动等。

自由液面流动问题涉及液体表面的卷曲、破碎和融合等复杂的几何非线性，用拉格朗日和欧拉类方法模拟均较为困难。

课题组提出了弱可压物质点法[16]和不可压物质点法[17]两种方案来模拟自由液面的流动问题。

在弱可压物质点法中，质点的压力通过人工状态方程[18]来计算。

弱可压物质点法能较好地模拟溃坝、楔形体入水和晃动等问题[16，20]，但遇到了如下困难,(1) 为了将流体密度变化控制在1%之内，人工声速至少要大于流场中最大流速的10倍，
导致显式时间积分的临界值过小，计算成本偏高；(2) 没有跟踪物质界面，难以准确施加压力边界条件；(3) 质点在穿越网格时会导致数值振荡，弱可压状态方程放大了网格穿越误差，甚至会导致非物理飞溅。

为了解决弱可压物质点法的上述缺点，基于算子分裂技术[21，22]提出了完全不可压物质点法[17]。

利用算子分裂技术可将更新拉格朗日格式的控制方程
分裂为
式(14)只和剪切变形有关，因此显式积分的临界时间步长取决于流体的剪切波速。

流体的剪切波速远小于压力波速，因此式(14)可以利用显式积分采用大时间步长高效求解。

将求解式(14)得到的试探速度场v*i作为式(15)的初始条件，并要求由式(15)得到的速度场应满足散度为0的条件，可得到压力泊松方程
将背景网格作为差分网格，采用中心差分法对式(16)进行离散，得到离散压力泊松方程组，可以利用预处理共轭梯度法高效求解。

另外，采用level set函数隐式跟踪自由表面，并采用虚拟网格法施加压力边界条件(含表面张力)，同时引入沙漏阻尼来抑制压力-速度半交错格式产生的沙漏模态[17]。

3.1 与有限元法的耦合
物质点法在每个时间步开始都会抛弃上个时间步已经变形的背景网格，重新采用新的规则背景网格，因而避免了有限元法中的网格畸变问题，可以很方便地处理极端变形问题。

物质点法可以看成是采用质点积分且在每个时间步中都进行网格重划分的有限元法。

与高斯积分不同，质点积分无法高效准确地积分多项式，因而对于小变形问题，物质点法的精度和效率均低于采用单点高斯积分的有限元法。

在很多工程问题中，极端变形可能只发生在局部区域，其他部分仍然是小变形。

将物质点法和有限元法按一定方式进行耦合，如仅在大变形区域采用物质点法，而小变形区域仍采用有限元法，可以充分发挥二者的优势。

3.1.1 物质点有限元法
有限元法在求解超高速碰撞问题时，为了模拟材料破碎和解决网格畸变，通常采用非物理的单元侵蚀算法，不但无法模拟超高速碰撞薄板形成碎片云的过程，且会造成质量、动量和能量不守恒。

针对有限元法在大变形区域网格畸变的问题，张雄等[23]借鉴物质点法的思想，在大变形区域引入背景网格，提出了物质点有限元法(MPFEM)。

计算过程中，所有物体仍用有限元离散，受背景网格覆盖的所有有限
元结点的动量方程均映射到背景网格上求解，其余结点的动量方程仍在原有限元网格上求解。

如果背景网格没有覆盖任何节点，则该方法退化为有限元法，如果背景网格覆盖所有结点，则该方法退化为物质点法。

为了保证计算精度和计算效率，一般只需在局部大变形区域布置背景网格。

以二维Taylor杆碰撞问题为例，如图10所示。

杆的底部在冲击载荷作用下将产
生很大的变形，导致严重的网格畸变。

为了解决单元畸变问题，可在杆的底部区域引入背景网格。

图11给出了MPFEM的基本原理，图中结点a～f是有限元结点，结点g～i是质点(进入背景网格区域的有限元节点)，A～D是背景网格上的结点，称为网格点。

当有限元结点在背景网格区域以外时，动量方程仍在有限元网格上求解，如a～f。

对于背景网格覆盖的材料子域，如结点g～i，动量方程在背景网格上求解。

求得
背景网格结点上的变量后，将其映射回质点以更新质点信息，然后丢弃变形后的背景网格，在下一时间步重新使用规则的背景网格，从而避免了网格畸变。

在大变形区域，时间步长取决于背景网格尺寸，不再是畸变的网格，因而提高了计算效率。

该算法适合于超高速碰撞等问题的模拟。

由于物质点法可以自动处理无滑移接触，所以在背景网格内发生碰撞的两物体不会相互嵌入，无需引入接触算法。

另外，在物质点法中，当两个质点的间距大于背景网格间距时，它们之间没有相互作用，可以作为对破坏现象的粗略描述。

3.1.2 耦合物质点有限元法
在中低速冲击侵彻和流固耦合类问题中，弹体和固体的变形一般较小，而靶体和流体的变形可能很大。

针对此类问题，廉艳平等[24]建立了耦合物质点有限元法CFEMP(Coupled Finite Element Material Point Method)。

该方法采用物质点
法离散大变形物体，用有限元法离散小变形物体，并通过基于背景网格的接触算法实现两者的耦合。

考察图12所示的两个物体r和s的接触问题。

Ωr和Ωs分别表示物体r和s所占据的区域，Γr
和Γs表示物体r和s的边界，Γc=Γr∩Γs表示接触边界，接触界面条件为非嵌透条件和接触面力条件。

假设物体r采用物质点离散，物体s采用有限元离散，如图13所示。

在每个时间
步中，首先假设接触没有发生，求得两个物体动量方程的试探解；然后通过背景网格进行接触探测，若两物体发生接触，则计算接触力并修正试探解以消除嵌透，若未发生接触，则试探解即为真实解。

在CFEMP中，分别采用有限元法和物质点法模拟求解不同的物体，以吸收二者各自的优点。

此外，耦合接触算法还可以模拟界面的接触/滑移/分离。

在最初的耦合物质点有限元法中，物质点和有限单元之间的耦合是基于背景网格的接触算法实现的，因此要求有限元区域和物质点区域之间的网格离散尺寸相互匹配。

该限制可能导致有限元区域网格过密，从而显著减少时间步长，增加计算成本和数据存储量。

为了解决网格匹配限制问题，陈镇鹏等[20]提出了耦合物质点有限元的改进算法。

该方法采用点面接触算法实现物质点粒子和有限单元之间的耦合，而非采用基于背景网格的接触算法，因此有限元网格的离散尺寸不再需要和背景网格匹配。

数值研究表明，改进的耦合物质点有限元法(ICFEMP)稳定性更好，效率和精
度均比纯物质点法和耦合物质点有限元法更高。

如在模拟楔形体入水问题时，ICFEMP使用不同有限元尺寸的模型均获得了理想的结果，但CFEMP在使用较大
尺寸的有限元模型时发生了严重的粒子穿透现象，如图14所示。

3.1.3 自适应物质点有限元法
在许多工程实际问题(如超高速碰撞和侵彻)中，远离载荷作用点的材料区域并没有发生大变形，材料大变形和破碎等现象仅发生在物体局部区域内；另外，初始时全部材料区域均为小变形状态，而在加载过程中某些局部区域逐渐发生大变形直至破碎。

3.1.1节的MPFEM需要事先指定可能发生大变形的区域，且不区分材料的变形状态，即使进入背景网格区域的单元仍处于小变形状态也会强制转化为物质点。

3.1.2节的CFEMP分别用物质点法和有限元法离散不同的物体，而不是根据各区域的当前变形状态采用相应的离散方法。

为了充分发挥物质点法和有限元法各自的优势，廉艳平等[23,24]充分吸收MPFEM和CFEMP两种方法的思想，提出了自适应物质点有限元法AFEMP(Adaptive Finite Element Material Point Method)[25,26]。

初始时所有物体均用有限元离散(图15(a))，当某单元的变形满足给定条件时自动将其转化为物质点(图15(b))。

同一个物体内物质点区域和有限元区域通过背景网格实现内力的传递(与MPFEM类似)，而不同物体之间的相互作用通过接触算法实现(与CFEMP类似)。

AFEMP充分发挥了MPFEM和CFEMP 的优势，非常适合模拟特大变形和破碎等问题。

与CFEMP类似，AFEMP也要求有限元网格尺寸和背景网格匹配。

陈镇鹏等[20]利用点面接触算法对自适应有限元物质点法进行了改进，消除了背景网格和有限元网格的尺寸匹配限制，同时也具有更高的精度和效率。

3.1.4 杂交物质点有限元法
钢筋混凝土在工程中有着广泛的应用，钢筋在混凝土的受拉承载中起着非常重要的作用。

一般钢筋直径的尺寸较小而混凝土结构尺寸较大，若钢筋和混凝土均采用物质点离散将导致计算规模过于庞大。

为了高效模拟钢筋混凝土问题，廉艳平等[27]。