冲击爆炸问题的三维物质点法数值仿真

张雄，廉艳平，杨鹏飞，李金光，张衍涛，王汉奎，刘岩

 

摘 要：材料与结构在冲击爆炸载荷作用下的动态响应问题是一个非线性的多物理场强耦合问题。物质点法采用拉格朗日和欧拉双重描述，吸收了拉格朗日法和欧拉法的优点，适合于分析特大变形及流动问题，是模拟超高速碰撞和爆炸问题的一种有效方法。本文简要地介绍了物质点法的基本理论和本课题组针对冲击爆炸问题在物质点法算法方面的研究工作，着重介绍了本课题组开发的冲击爆炸问题的三维显式并行物质点法数值仿真软件MPM3D^®及其应用。MPM3D^®采用C++语言研发，并基于VTK和Qt开发了图形用户界面PeneBlast，可运行于Windows、Linux和Mac OS等多种平台下。

 

关键词：物质点法；无网格法；超高速碰撞；侵彻；爆炸；

 

3D Simulation Based on Material Point Method for Impact and Explosion Problems

ZHANG Xiong，LIAN Yanping，YAN Penfei，LI Jinguang，ZHANG Yantao, WANG Hankui, LIU Yan

 

Abstract: Dynamic response of material and structure under impact and blast loading involves multi-physics coupling and nonlinear relations. Material Point Method (MPM) consists of both Lagrangian and Eulerian description of material and is suitable for problems with extreme large deformation by taking advantage of the merits of both Lagrangian method and Eulerian method. The basic formulation of MPM and our extension on MPM for impact and explosion problems are briefly presented. Our 3D explicit parallel MPM code, MPM3D, and its applications in impact and explosion problems are described. MPM3D is developed by using C++, and its GUI PeneBlast is developed based on Qt and VTK, and can be run on different platforms including Windows, Linux and Mac OS.

Key words: material point method; meshfree method; impact; explosion

0引言

材料与结构在冲击爆炸载荷作用下的动态响应问题是一个几何、材料和边界条件均为非线性的多物理场强耦合问题，涉及高应变率、高压、高温、相变乃至化学反应，气体、液体和固体等多种物质间相互耦合甚至混合，材料不但会发生严重扭曲和破碎，还会融化甚至汽化。在这类问题的研究中，数值模拟和其他研究手段相比，更加经济，便于观察，也能突破实验和理论研究的局限性，正在发挥着越来越重要的作用。

数值模拟方法可分为拉格朗日方法和欧拉方法两大类。拉格朗日法中计算网格随物质一起变形，可方便地跟踪材料界面和引入与变形历史相关的材料模型，但对于涉及特大变形的问题会因网格严重畸变而产生数值求解困难，且难以有效地模拟材料的破碎、融化和汽化等行为。此类方法代表性程序为DYAN。欧拉法中计算网格固定在空间中，不存网格畸变问题，但不易跟踪材料界面，且非线性对流项也会导致数值求解困难。国内学者在此方面做出了大量工作，宁建国等开发了欧拉型三维爆炸与冲击问题数值模拟软件EXPLOSION-3D^[1]，北京计算数学与应用物理研究所开发了流体动力学程序MEPH2Y和MEPH3D ^[2,3]，刘凯欣等基于高精度时空守恒元解算法开发了SUPER CE/SE程序^[4,5]。

近年来，无网格算法在国内引起了广泛的关注并将其应用于冲击爆炸领域。刘桂荣及刘谋斌运用光滑质点流体动力学方法（SPH）对爆炸问题做出了大量的工作^[6]。杨秀敏等运用SPH并结合有限元开发了爆炸冲击现象数值模拟软件EF3D和EP3D^[7]。王宇新等运用物质点法（MPM）模拟了冲击爆轰等问题^[8,9]。

本文简要介绍了物质点法和本课题组针对超高速碰撞、侵彻和爆炸问题在物质点法算法方面的研究工作，重点介绍了本课题组研发的三维显式并行物质点法数值仿真软件MPM3D^®及其在超高速碰撞、侵彻和爆炸等方面的应用。MPM3D^®采用C++语言编制，并基于 Qt、VTK和CMake开发了相应的图形用户界面系统PeneBlast，可运行于Windows、Linux和Mac OS平台下。

1   物质点法

基于FLIP(Fluid Implicit PIC)方法^[1]发展起来的物质点法(material point method, MPM) ^[11]采用拉格朗日和欧拉双重描述，将物体离散为一组在空间网格中运动的质点。这些质点携带所有物质信息，如质量、速度、应变和应力等，其运动代表了物体的变形和运动。空间网格（即背景网格）用于动量方程的求解和空间导数的计算，不携带任何物质信息。背景网格可以在空间中固定，也可以根据问题的特点按照某种方式自由布置。在每一个时间步中，首先将物质点和背景网格完全固连，通过将物质点的物理量映射到背景网格结点上建立背景网格结点的运动方程，求解运动方程后再将背景网格结点的物理量映射回物质点，得到下一时刻物质点所携带的物质信息。这一步物质点和网格结点一同变形，没有相对运动，是拉格朗日求解，避免了欧拉法中处理对流项的困难，且很容易跟踪物体的界面。物质点携带了所有物质信息，因此在下一个时间步中可以丢弃变形后的背景网格，仍采用规则的背景网格，从而避免了拉格朗日法中因网格畸变而产生的数值困难。可见，MPM结合拉格朗日方法和欧拉方法各自的优势，适合于分析特大变形及流动问题。

MPM采用线性形函数，其导数在格子(cell)之间不连续，因此当质点跨越格子时会对系统产生扰动，称之为跨格子噪音（cell-crossing noise），导致精度降低，甚至使部分问题结果完全失真。Bardenhagen等^[12]针对此问题，发展了广义插值物质点法（Generalized Interpolation Material Point method，GIMP），有效地抑制了由于质点跨越格子而引起的数值噪音。

由于速度场的单值性，MPM自动满足无滑移接触条件，便于处理碰撞接触问题，特别是在超高速碰撞问题的模拟中，弹靶间的滑移可以忽略，采用MPM不需采用特别的接触算法即可达到较好的模拟效果。为了处理滑移接触， Hu等人^[13] 引入了多重背景网格的概念，各物体在接触结点处的法向运动量在公共背景网格上求解，而切向运动量则是在各物体的背景网格上求解的，既保证了各物体在接触点处的法向速度和法向加速度连续，又允许它们的切向运动相互独立。Bardenhagen等人^[14]基于多重速度场的思想给出了物质点方法的接触/摩擦/分离算法，用于颗粒材料的数值模拟。Nair等人^[15]将接触算法扩展到物体内部表面之间的自接触中，用于分析裂纹扩展问题。

物质点法目前已应用于冲击、爆炸、裂纹扩展、材料破坏(脆性破坏、分层复合材料脱层破坏、层裂破坏)、多孔材料与颗粒材料、生物力学等的数值模拟中。

2   本课题组对物质点法的扩展

2.1 超高速碰撞

在空间碎片防护设计中，研究碎片云对航天器的作用是非常重要的。通过侵蚀算法，有限元法可以模拟弹丸对靶板的侵彻过程，但无法得到超高速碰撞所形成的碎片云的形貌。通过引入Johnson-Cook材料模型和Mie- Gruneisen状态方程，我们将MPM扩展应用于弹丸超高速碰撞薄板问题，较好地模拟了薄板背后形成碎片云的形貌，展示了物质点法在这一领域的应用前景^[16]。

 

2.2 物质点有限元法

物质点法虽然避免了拉格朗日有限元法的网格畸变困难，但对于小变形问题其效率和精度均低于有限元法，因此有必要将两种方法结合。我们通过在大变形区域中引入背景网格，提出了物质点有限元法(Material Point Finite Element Method, MPFEM) ^[17]。物体初始用有限元离散，但在可能发生大变形的区域中预先布置一背景网格。将进入背景网格区域的有限元结点视为物质点，动量方程在背景网格上求解，而在其他区域中动量方程仍在有限元网格上求解。物质点有限元法充分发挥了物质点法处理特大变形的能力和有限元法的高效性。

 

2.3 杂交物质点有限元法

在将MPM扩展应用于钢筋混凝土的侵彻问题时，钢筋的尺寸相对于混凝土的尺寸较小给离散带来了困难。为此，我们将有限单元（如杆单元）引入到MPM中，提出了杂交物质点有限元法^[18]。单元的节点如质点一样，其动量方程在背景网格上积分，但其本构计算则在单元的高斯点上进行计算而非质点上。通过引入杆单元，将其成功应用于钢筋混凝土的侵彻问题的数值模拟中。

 

2.4 质点自适应分裂

在MPM中，当两个物质点间的距离大于背景网格的结点间距时，它们之间将没有相互作用，从而发生数值断裂。针对这一问题，我们提出了MPM的自适应分裂质点方案^[19]，使质点可以根据累计应变，在变形最大的方向上自适应分裂，从而避免了数值断裂，并将其成功应用于激波管、聚能射流等爆炸问题的数值模拟。

 

2.5 接触算法

对于中低速侵彻问题，标准物质点法中的无滑移接触约束产生较大的侵彻阻力，因此在分析侵彻问题时需要引入接触条件^[20]。我们建立了基于局部多重背景网格的接触算法^[21]，即仅在接触区域建立局部多重背景网格并引入多重速度场的思想在背景网格上调整速度场并求解接触力。同时，为了避免MPM中提前接触的现象，我们给出了计算即将接触的物体间物理距离的算法。该方法减少了对存储空间的需求并提高了接触算法的效率改善了接触判断的精度。

 

2.6 并行物质点法

使用并行计算技术可以实现更大规模问题的模拟。MPI为分布式内存并行技术，适合于在Cluster机器上运行。我们基于MPI给出了针对Cluster机器的MPM并行算法，该算法涉及到局部背景网格分解，每一个时间步均需要进行子区域重复节点之间的通信，而且还要建立复杂的并行I/O。OpenMP为共享内存并行技术，适合于在多核计算机或者SMP机器上运行。我们基于OpenMP提出了针对共享内存计算机的MPM并行算法^[20]。由于OpenMP支持增量并行，而且不会涉及到复杂的区域分解，因此已经逐步在显式有限元计算、分子动力学计算和流体力学计算并行化中得到应用。

3   物质点法数值仿真软件MPM3D

在算法研究的基础上，本课题组采用C++语言研发了冲击爆炸问题的三维显式并行物质点法数值仿真软件MPM3D^®，并为其开发了图形用户界面系统PeneBlast，如图1所示。两者均采用了跨平台的开发工具和程序库，可运行于Windows、Linux和Mac OS等多种平台下。

 

图 1 PeneBlast 主界面图

Figure 1 The main interface of PeneBlast software

 

3.1 图形用户界面(GUI)系统—PeneBlast

MPM3D^®的图形用户界面(GUI)系统PeneBlast是基于QT和VTK开发的。Qt是一个跨平台的C++图形用户界面应用程序框架，具有优良的跨平台特性，模块化程度高，可重用性较好。VTK是一个开放源码的跨平台可视化工具包，可用来进行3D计算机图形、图像处理和可视化。

PeneBlast主要用于前处理和模型求解。通过PeneBlast，用户可以完成从建立三维实体模型、物质点离散、从材料库中选择材料、施加边界条件、设置求解参数到生成供求解器读取的输入文件这一系列前处理过程。其中材料库中包含了常用的材料参数，用户可以直接在材料库中选择使用某种材料用于仿真模拟，也可拷贝、修改或添加新的材料到材料库中。

生成输入文件后，可以在界面上调用求解器MPM3D对所建立的物质点模型进行求解，并实时监控整个求解过程；可以暂停/重启动及中断求解器的运行。

通过PeneBlast生成的供MPM3D读取的输入文件采用xml格式，易于扩充和读取。同时，PeneBlast还可导入已有的输入文件并在此基础上继续建模，PeneBlast也可对输入文件进行格式检查，并直接编辑输入文件。

后处理采用了开源的多平台数据分析和可视化应用程序ParaView，它可定性或定量分析数据并实现可视化。MPM3D^®将计算结果保存为h5格式和VTK PolyData格式，供ParaView进行数据分析与可视化处理，绘制和导出变形动画、云图，绘制时间历程曲线，提取某些点的变量值等。

 

3.2 求解器MPM3D^®的主要功能

MPM3D^®采用C++语言编写，可以使程序具有更好的可维护性和良好的可扩充性，且易于添加新类以增加新的功能。

目前，MPM3D^®中具有多种材料模型，如弹性及弹塑性模型、Moony-Rivlin超弹性模型、Johnson-Cook模型、Drucker-Prager岩土模型、HJC混凝土模型、高能炸药模型、Gurson模型以及可用于流体模拟的空材料模型等。冲击爆炸问题中涉及高温高压高应变率，相对低速下的材料响应而言，热力学效应更加明显，因此需要使用状态方程来描述压强、密度和内能之间的关系。MPM3D^®中包含用于高压固体模拟的Gruneisen状态方程、理想气体状态方程、多项式状态方程，用于爆轰产物模拟的JWL状态方程和用于模拟混凝土的HJC状态方程及Palpha状态方程。

为了处理材料的断裂破坏，MPM3D^®具有最大静水拉力失效、最大等效塑性应变失效、最大主应力/剪应力失效，最大主应变/剪应变失效以及Johnson-Cook损伤失效模型。为了处理材料的随机特性导致其损伤破坏的随机性，MPM3D^®支持随机失效设置，目前支持Gauss和Weibull两种随机分布。

MPM3D^®中爆炸采用基于CJ理论的理想爆轰，起爆方式为程序起爆。用户可以定义多个起爆点，也可以定义起爆面。为了处理冲击波，添加了人工体积粘性。

MPM3D^®可以对某物体设置初速度或常速度，或者施加随时间变化的载荷，也可施加加速度场，如重力场等。

MPM3D^®的背景网格可固定在空间不动，也可以根据物体所占的区域动态变化。用户可以设置是否使用GIMP算法，是否启用接触算法。时间积分为可变时间步长的显式积分，支持USF、USL或MUSL格式^[15]。用户即可采用OpenMP并行设置实现在共享内存计算机上的并行计算，也可以采用MPI并行设置实现在计算机集群上的并行计算。

4   应用实例

我们利用MPM3D程序计算了大量的验证和工程算例。在此，我们仅给出有代表性的算例。为了节省篇幅，我们略去了计算方案、材料参数等详细内容。

 

4.1 超高速碰撞

4.1.1弹丸超高速碰撞薄板—碎片云

一半径为7.5mm、质量为20g的铅弹以6.56km/s的速度撞击厚度为6.35mm的铅靶。图2比较了30.6 ms时碎片云的数值模拟结果和实验结果^[24]。由MPM3D得到的30.6 ms时碎片云前端和靶板的距离和碎片云宽度分别为198 mm和142 mm，其实验结果分别为200 mm和145 mm。

     

图2  30.6 ms时刻的数值结果和实验结果比较

Figure 2 Comparison of debris cloud at 30.6 ms between MPM3D and experiment results

 

4.1.2弹丸超高速碰撞厚板—成坑

一质量为0.5g球形铜弹丸以6.0km/s的速度撞击厚度为40mm，长宽均为80mm的铜靶。表1比较了弹坑参数的数值模拟结果与实验结果^[25]，两者吻合较好。图3给出了32ms时刻碰撞变形结果的切片图。

表 1 32ms时刻弹坑参数计算结果与实验结果比较：

Table 1：Comparison of the size of crater at 32 ms between MPM3D and experimental results

	弹坑深度（mm）	弹坑直径（mm）	坑深/坑径
实验	14	25.4	0.55
MPM3D	14	25	0.56

 

图3  32ms时刻铜弹丸碰撞厚板形状

Figure 3  The deformed shape of thick plate and sphere at 32ms

 

4.2 侵彻问题

4.2.1    钢珠侵彻薄板

一直径为为10mm的钢珠正侵彻一厚度为1mm，直径为178mm的薄圆板。图4比较了靶体的最终变形数值模拟结果与实验结果^[26]。实验所得靶体的变形高度与靶体的穿孔直径比值h/D为0.84，MPM3D计算得到的值为0.85。

(a)实验结果

 

(b) MPM3D 计算结果

图4 靶体最终变形结果

Figure 4  Final deformed shape of the target

 

4.2.2金属靶斜侵彻

一长为88.9 mm，直径为12.9 mm的尖拱钢弹以575m/s的速度撞击倾角为30°、厚度为26.3mm的铝合金靶。图5比较了194.4ms时刻的数值模拟结果和实验结果^[27]。由MPM3D得到的剩余弹速为453.4 m/s，与实验值455 m/s吻合很好。

    

 

(a) 实验结果        (b) MPM3D 结果

图5  194.4 ms时刻尖拱弹侵彻铝靶变形结果

Figure 5 The deformed shape of ogive-nose projectile impacting aluminum target at194.4 ms

 

4.2.3钢筋混凝土靶正侵彻

一长度为144mm，直径为25.4，尖拱钢弹正侵彻厚度为178mm的钢筋混凝土靶板。图6比较了弹体以不同速度入射时穿过靶体后的剩余速度的数值模拟结果与实验结果^[28]，图7给出了当弹体以初速度749m/s打中钢筋时的钢筋混凝土及钢筋破坏情况。

图6  剩余弹速比较

Figure 6 Comparison of the residual velocity of projectile

(a)  靶体的破坏情况

(a) 钢筋的破坏情况

图 7  钢筋混凝土的损伤结果

Figure 7 The damage shape of the reinforced concrete

 

4.3 爆炸问题

4.3.1    激波管

激波管由薄膜分为两个初始静止的、具有不同密度和压力的区域，如图8所示。在时间t > 0时，薄膜破裂，激波和接触面会以不同的速度从左向右运动。通常观察t = 0.143时刻的结果，此时激波前进了大约0.25。

图8  激波管示意图

Figure 8 The sketch map of shock tube problem

MPM3D采用GIMP算法给出的t = 0.143时刻密度、速度、压力和比内能的分布曲线与解析解吻合很好，如图9所示。

图9 密度、速度、压力和比内能MPM3D结果与解析解(实线)的比较

Figure 9 Profiles of density, velocity, pressure and energy obtained by analytical solution (solid line) and MPM3D simulation

 

4.3.2    一维板条TNT爆轰

一维爆轰过程经常作为测试爆轰模拟的数值方法的基准检验问题，长为0.1m的板条TNT炸药在一端起爆，并在起爆端固定。图10给出了由0至14ms时间内的压力、密度和速度沿板条的曲线分布图，相邻曲线时间间隔为1 ms。数值模拟所得各曲线的峰值与解析解吻合较好。

图 10在起爆过程中沿着板条的压力曲线

Figure 10  Profiles of pressure along the slab during the detonation process

 

4.3.3    爆轰驱动飞片

爆轰驱动飞片在多个领域有着广泛的应用，该过程涉及爆轰波的传播、爆轰产物的急剧膨胀及其与周围材料相互作用产生的大变形。图11为非对称构型装药一维和二维模型的飞片终态速度与炸药爆轰速度比值随质量比变化的结果，数值模拟结果与Gurney公式及其修正的经验公式^[29][30]预测值吻合很好。

图11 不同质量比下的无量纲飞片终态速度

Figure 11 Dimensional terminal velocity of flyer with different mass ratio

 

4.3.4    爆炸碎片模拟

金属壳在爆炸作用下生成的碎片极具破坏性，研究其机理对公共安全有着重要的意义。我们采用Gurson模型及Weibull随机失效机制来模拟材料的断裂破碎过程的随机性。图12为模拟得到的柱壳破碎过程，图13为三维球壳在爆炸作用下产生的裂纹和碎片形状。

         

     (a) 10ms       (b) 20ms        (c) 30ms        (d) 40ms

图12 采用Gurson模型的柱壳破碎过程

Figure 12 The crushing process of cylindrical shell based on Gurson model

图13 球壳的破碎形态

Figure 13 The crushing shape of spherical shell

 

4.3.5 聚能射流模拟

聚能射流在穿甲、石油开采等领域有着重要的应用。在此，我们模拟的聚能装药的初始构型如图14所示，半锥顶角为38度。炸药为TNT，外壳和药型罩材料为铜。

图14 聚能装药的初始构型

Figure 14 The initial configuration of shape charge

射流尖端的速度和杵的速度可以由Birkhoff 等人提出的理论模型进行估计。根据理论估计^[29]，本例中射流尖端速度应为3.7 km/s，杵的速度应为0.6 km/s。图15(a)显示了模拟得到的21ms时的x方向速度分布，和理论估计基本吻合。图15 (b)显示了射流中温度的分布情况。

(a) 射流中的速度分布

(b) 射流中的温度分布

图15  21ms时刻射流中的速度分布和温度分布

Figure 15 The distribution of velocity and temperature of jet at 21ms

 

4.4    其它问题

4.4.1 边坡失效问题

在岩土工程中，常采用铝棒堆积物的流动模拟沙和土的颗粒流动。我们用Drucker-Prager模型对Bui进行的干燥铝棒堆积物的坍塌流动实验^[31]进行了模拟。图16(a)和（b）比较了数值模拟结果与实验结果。图16（c）将实验得到的失效区和计算结果进行了定量比较，表明坡面及坡体失效区的数值模拟结果与实验结果吻合。

(a) 实验结果

图 16 无粘沙土失效面MPM3D结果与实验结果对比

Figure 16 Comparison between MPM3D simulation and experiment result of no-cohesive soil failure

 

4.4.2 金属切削问题

切削加工技术是工业装备过程中最为重要和普遍的技术之一，切屑过程伴随着高速、高应变率、断裂和高温，极易产生绝热剪切带进而产生锯齿状切屑。我们模拟了一简单的切屑过程，如图17所示。在0.01ms时第1条剪切带触发，随着更多的剪切带触发，锯齿状切屑形成。

      

0.01ms                   0.03ms

      

0.05ms                 0.07ms

图17 切削过程

Figure 17 The processing of cutting

 

4.3碳纳米管复合材料

碳纳米管有着突出的力学性能,将其掺杂于基体中可形成高性能的复合材料。但碳纳米管细长比大，在基体中的均是以弯曲的形式存在，这给基于网格的建模带来极大的困难。在MPM3D中，利用塌落生成的碳管的几何模型布点，可构建出碳管的计算模型；管间正交布点，可构建出基体材料模型。

图18为建好的碳纳米管复合材料的计算模型，基体材料点未画出，碳管间黑点用于代表碳管间近距离相互作用。图19为该模型受剪力作用下的力-位移曲线，可得出碳管间相互作用对复合材料整体模量的影响。

图18 考虑碳管近距离作用的离散

Figure 18  Discrete model with considering near-range action between carbon tubes

图19 管间不同连接强度的力-位移曲线

Figure 19  Shear force and displacement relationship for different connections between tubes

 

5. 结论

本文介绍了冲击爆炸问题的显式三维物质点法数值仿真软件MPM3D及其图形用户界面系统PeneBlast，可运行于Windows、Linux和Mac OS等多种平台下。大量的超高速碰撞、侵彻爆炸数值算例说明了该软件的可靠性和准确性。该软件可以用于超高速碰撞、侵彻和爆炸等问题的数值仿真，为航天器空间碎片防护设计、常规武器研发与防护工程设计等领域提供有效的数字化设计工具。

致  谢

感谢已毕业博士生马上、黄鹏及出站博士后马志涛对本课题所做出的出色工作。

参考文献

1. 宁建国，王成，马天宝. 爆炸与冲击动力学[M]. 北京：国防工业出版社，2010.
2. 刘军，何长江，梁仙红. 三维弹塑性流体力学自适应欧拉方法研究 [J]. 高压物理学报，2008,22(1): 155-178.
3. 何长江，于志鲁，冯其京. 高速碰撞的三维欧拉数值模拟方法 [J]. 爆炸与冲击，1999,19(3): 216-221.
4. 王景焘，张德良，刘凯欣. 基于CE/SE方法的二维Euler型多物质流体弹塑性问题计算 [J]. 计算物理, 2007,24(4): 395-401.
5. Wang J, Liu K, Zhang D. An improved CE/SE scheme for multi-material elastic-plastic flows and its applications [J]. Computers & Fluids, 2009, 38(3): 544-551.
6. 刘桂荣，刘谋斌. 光滑粒子流体动力学 — 一种无网格粒子法[M].韩旭 等译. 长沙：湖南大学出版社,2005.
7. 杨秀敏. 爆炸冲击现象数值模拟[M]. 合肥：中国科学技术大学，2010.
8. 王宇新，陈震，孙明. 滑移爆轰问题无网格MPM法模拟 [J]. 力学与实践, 2006,29(3): 20-25.
9. 王宇新，顾元宪，孙明. 无网格MPM法在冲击载荷问题中的应用 [J]. 工程力学, 2006,23(5): 46-51.
10. Brackbill, J U, Kothe D B, and Ruppel H M. FLIP: a low-dissipation, particle-in-cell method for fluid flow [J]. Computer Physics Communications, 1988, 48: 25-38.
11. Sulsky D, Chen Z, Schreyer H L. A Particle Method for History-Dependent Materials [J]. Computer Methods in Applied Mechanics and Engineering, 1994, 118(1-2): 179–196.
12. Bardenhagen S G, Kober E M. The generalized interpolation material point method [J]. CMES-Computer Modeling in Engineering & Sciences, 2004, 5(6): 477–495.
13. Hu W, Chen Z. A multi-mesh MPM for simulating the meshing process of spur gears [J]. Computers and Structures 2003, 81:1991-2002.
14. Bardenhagen S G, Brackbill J U, Sulsky D. The material point method for granular materials [J]. Computer Methods in Applied Mechanics Engineering 2000, 187:529-541.
15. Nairn J A. Material point method calculations with explicit cracks [J]. CMES-Computer Modeling in Engineering & Sciences, 2003, 4(6):649–663.
16. 马上, 张雄, 邱信明. 超高速碰撞问题的三维物质点法[J]. 爆炸与冲击, 2006, 26(3): 273–278.
17. Zhang X, Sze K Y, Ma S. An explicit material point finite element method for hyper-velocity impact [J]. International Journal for Numerical Methods in Engineering, 2006, 66: 689-706
18. Lian YP, Zhang X, Zhou X. Ma ZT. A FE-MP method and its application in modeling dynamic response of reinforced concrete subjected to impact loading [J]. Computer Methods in Applied Mechanics and Engineering (2011). In press.
19. Ma S, Zhang X, Lian Y, Zhou X. Simulation of high explosive explosion using adaptive material point method [J]. CMES-Computer Modeling In Engineering & Sciences, 2009, 39: 101-123
20. Huang P, Zhang X, Ma S, et al. Contact Algorithms for the Material Point Method in Impact and Penetration Simulation [J]. International Journal for Numerical Methods in Engineering, 2011, 85(4):498-517.
21. Ma ZT, Zhang X, Huang P. An Object-Oriented MPM Framework for Simulation of Large Deformation and Contact of Numerous Grains [J]. CMES: Computer Modeling in Engineering & Sciences, 2010, 55(1): 61-88.
22. Huang P, Zhang X, Ma S, et al. Shared Memory OpenMP Parallelization of Explicit MPM and Its Application to Hypervelocity Impact [J]. CMES: Computer Modeling in Engineering & Science, 2008, 38(2): 119–147.
23. Zhang YT, Zhang X. An alternated grid updating parallel algorithm for material point method using OpenMP [J].  CMES: Computer Modeling in Engineering & Sciences, in press.
24. Anderson C E J, Trucano T G, Mullin S A. Debris cloud dynamics. International Journal of Impact Engineering, 1990, 9: 89-113.
25. C.A. Wingate, R.F. Stellingwerf, R.F. Davidson et al., Models of high velocity impact phenomena [J]. International journal of impact engineering, 1993, 14: 819-830.
26. Seoa S W, Min O K, Lee J H. Application of an improved contact algorithm for penetration analysis in SPH [J]. International Journal of Impact Engineering, 2008, 35: 578-588.
27. Piekutowski A J, Forrestal M J, Poormon K L, Warren T L. Perforation of aluminium plates with ogive-nose steel rods at normal and oblique impacts [J]. International Journal of Impact Engineering 1996, 18:877- 887.
28. Hanchak S J, Forrestal M J, Young E R, et al. Perforation of concrete slabs with 48 MPa and 140 MPa and 140 MPa unconfined compressive strengths [J]. Int. J. Impact Eng, 1992, 12(1): 127.
29. Mayers M A. Dynamic Behavior of Materials [M]. John Wiley & Sons, 1992
30. Zukas J A, Walters W P, (eds.). Explosive Effects and Applications [M]. Springer, 1998.
31. Ha H. Bui, Ryoichi Fukagawa, et al. Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model [J]. International Journal for Numerical and Analytical Methods in Geomechanics. 2008; 32:1537-1570.