# Numerical Investigation on Improving the Computational Efficiency of the Material Point Method

Authors: Jingxin Ma; Shengzhuo Lu; Yaqin Shi; Keyi Chen; Weidong Chen
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Based on the basic theory of the material point method (MPM), the factors affecting the computational efficiency are analyzed and discussed, and the problem of improving calculation efficiency is studied. This paper introduces a mirror reflection boundary condition to the MPM to solve axisymmetric problems; to improve the computational efficiency of solving large deformation problems, the concept of "dynamic background domain (DBD)" is also proposed in this paper. Taking the explosion and/or shock problems as an example, the numerical simulation are calculated, and the typical characteristic parameters and the CPU time are compared. The results show that the processing method introducing mirror reflection boundary condition and MPM with DBD can improve the calculation efficiency of the corresponding problems, which, under the premise of ensuring its calculation accuracy, provide useful reference for further promoting the engineering application of this method.


---

## Body

## 1. Introduction

In this paper, we want to discuss how to improve the computational efficiency of axisymmetric problems and large deformation problems when the meshfree method material point method (MPM) is used. Meshfree/meshless methods have been progressed significantly and achieved remarkable progress in recent years [1, 2]. Smoothed particle hydrodynamics (SPH) is the first meshless method [3], and it can solve problems of solid mechanics [4], fluid dynamics [5, 6], machining [7], heat conduction [8], and so on. To remove the instabilities in SPH methods, a reproducing kernel particle method (RKPM) [9] and Meshless local-Petrov Galerkin (MLPG) [10] method were developed. The extended finite-element method (XFEM) was introduced by Belytschko and Black [11], as an extension of the standard FEM, and it has been used to deal with crack problems and so on [12–14]. Rabczuk and Belytschko [15] proposed the cracking particles method (CPM), and it has been used to model the dynamic crack propagation [16–18]. Also, isogeometric analysis (IGA) [14, 19] has recently surfaced as a potential technique in the field of solid mechanics analysis. In 1994, the EFG method [20] was developed and used a global weak form to solve a variety of problems whether under various loadings or conditions immaculately [21–23]. To deal with the problem of crack tracking, the concept of dual-horizon is introduced to Peridynamics (DH-PD) to consider the unbalanced interactions between the particles with different horizon sizes, which can solve the “ghost force” issue [24, 25]. The phase field modelling was also used to deal with the crack propagation [26–29]. A numerical manifold method was developed to deal with the Kirchhoff plate [30–32].The material point method is an extension of the fluid implicit particle (FLIP) [33]/particle in cell (PIC) [34] method, proposed by Sulsky, Chen, and Schreyer [35]. In comparison with other methods, the MPM takes advantage of both Lagrangian and Eulerian descriptions [36, 37], which have proven useful for solving impact and contact problems [38, 39], penetration [40, 41], explosion problems [42–44], combustion reaction [45], cracks and fracture [46–49], stochastic and random problems [50, 51], and so on. In the framework of the MPM, the background area is discretized by Eulerian grid and grid nodes, the material domain is discretized by a finite number of particles, and the calculation time is comprised of multiple time steps. There are two phases in the single time step: In the first phase, the physical information on material points are mapped to grid nodes by shape functions; after obtaining the kinematic solution on the grid nodes, the physical information is mapped back to material points to update their characteristic information such as positions and velocities. Hence, the MPM avoids the mesh distortion in the Lagrangian method and the convection problems in Eulerian methods [33, 34]. Also, the MPM shows many advantages over other meshless methods in tension stability and efficiency [36–38].Many interesting scenarios from the preceding categories lend themselves to being treated as axisymmetric. These include, for example, impact and cratering events, rod penetration, and shaped charge jet formation. The material point method proposed earlier treated particles as point masses, and it has long been recognized that this original MPM (i.e., particles represented spatially as Dirac delta functions) has serious element-crossing artifacts. In addition to reformulation of axisymmetry using GIMP methods, numerous other extensions to the MPM have developed since the original work on axisymmetric MPM [52, 53]. Nairn and Guilkey have established a series of new shape functions which differ mostly from planar shape functions near the origin where r = 0 [36]. In this paper, on the basis of not changing the right-angled coordinate system, the axisymmetric boundary conditions are introduced, and the three-dimensional problem is simplified to a two-dimensional plane problem for processing.The authors have been extensively researching the MPM as it appears to resolve the large deformation problems. However, like all novel techniques, it suffers from a few deficiencies that require further research. A common problem is that when using the MPM to calculate the problem, the physical values of the particles should be updated strictly according to the dimension characteristics of the problem; it means that the two-dimensional problems and three-dimensional problems are handled in different ways. Therefore this paper introduces the mirror reflection boundary condition on the basis of the theory of the material point method, which can simplify the three-dimensional axisymmetric problem into a two-dimensional problem and can solve the three-dimensional axisymmetric problems more efficiently. Furthermore, when dealing with the large deformation problems, the physical information of all material points requires two mappings from the material points to the background grid and the background grid to the material points within each calculation step in the original MPM [54–56]. In order to improve the computational efficiency, the concept of dynamic background domain (DBD) is presented in this paper, and it allows the method to skip the mapping process of background grids and material points which not are updated in the physical information during the calculation process. The algorithm, thereby, allows a more thorough search of all the materials point and seems to offer a greater ability to find global minima.This paper is organized as follows. The algorithm of the explicit MPM is reviewed briefly in Section2. The detailed formulation of the axisymmetric problems solution is presented in Section 3, which also describes the concept of dynamic background domain. The numerical simulations and several conclusions are given in Section 4. Finally, Section 5 draws a summary and some concluding remarks. The source codes of the mirror reflection boundary conditions and dynamic background domain are given in Algorithms 1 and 2.Algorithm 1: Source code for mirror reflection boundary condition.
module specular_reflectionuse param, only: parts, dy, partstotal, round, symmetry, reptotalimplicit nonecontainsopen(1,file = "evaluate.dat")subroutine evaluate()implicit noneinteger i, j, k, tif (round == 1) thenj = 0do i = 1, partstotalif (parts(i)%x(2) < = 3∗dy)thenj  =  j + 1k = partstotal + jsymmetry(j)%p(1) = isymmetry(j)%p(2) = kwrite(1, ∗) symmetry(j)%p(1), symmetry(j)%p(2)end ifreptotal = jt = jend doelset = reptotalend ifdo i = 1, tparts(symmetry(i)%p(2))%x(1) = parts(symmetry(i)%p(1))%x(1)parts(symmetry(i)%p(2))%x(2) = -parts(symmetry(i)%p(1))%x(2)parts(symmetry(i)%p(2))%vx(1) = parts(symmetry(i)%p(1))%vx(1)parts(symmetry(i)%p(2))%vx(2) = -parts(symmetry(i)%p(1))%vx(2)end doend subroutineclose(1)end moduleAlgorithm 2: Source code for dynamic background domain.
fxmax = 0.0fymax = 0.0zxmax = 0.0zymax = 0.0do i = 1, ptotalif(parts(i)%pres/= 0.0.or. parts(i)%v/= 0.0) thenif(fxmax.ge.parts(i)%x(1)) fxmax = parts(i)%x(1)if(fymax.ge.parts(i)%x(2)) fymax = parts(i)%x(2)if(zxmax.le.parts(i)%x(1)) zxmax = parts(i)%x(1)if(zymax.le.parts(i)%x(2)) zymax = parts(i)%x(2)end ifend dofxmax = fxmax-3.0∗dxfymax = fymax-3.0∗dyzxmax = zxmax + 3.0∗dxzymax = zymax + 3.0∗dyif(round==2) thendo i = 1, ptotalparts(i)%sign = 1end doend ifdo i = 1, ptotaldeallocate(parts(i)%nodes,s(i)%pnode,s(i)%s,s(i)%sx)if(parts(i)%sign==0) thenelseif(parts(i)%x(1)<=zxmax.and.parts(i)%x(1)>=fxmax.and.parts(i)%x(2)<=zymax.and.parts(i)%x(2)>=fymax) thenparts(i)%sign=0end ifend ifend do


## 2. Material Point Method

### 2.1. Governing Equations

In the MPM, the background grid is divided to calculate the momentum equation, and the material domainΩ is discretized into a collection of particles, as shown in Figure 1.Figure 1
Discrete schematic diagram of the MPM.For continuous materials, they are governed by the conservation of momentum equations:(1)ρdvdt=∇⋅σ+ρb,where, ρ=ρx,t is the particle density, v=vx,t is the velocity vector, σ=σx,t is the stress tensor, and b=bx,t is the body force.The mass conservation equation is(2)dρdt+ρ∇⋅v=0,and the conservation of energy equation is(3)ρdedt=σ:dεdt+ρv⋅b,where ε=εx,t is the strain tensor and dε/dt is the strain rate tensor.


### 2.2. MPM Solution Scheme

The first step in the single time step of the MPM is to extrapolate the particle momenta and masses to the background. The equations are(4)pi=∑p=1npmpvpNi,p,mij=∑p=1npmpNi,p,where pi is the nodal momentum, vp is the particle velocity, mp is the particle mass, Ni,p is the shape function for node i evaluated at the current position of particle p, and mij is the nodal mass in the lumped (or diagonal) mass matrix.In this representation, the acceleration vector at each grid point is found by(5)∑j=1nmijaj=fiint+fiext,where fiint and fiext are internal and external force vectors associated with node i, and they are defined by(6)fiint=−∑p=1npmpρpσxpt,t⋅∇Nixp,fiext=∑p=1npmpNixpbxp,t+∫Γττ¯xp,t⋅Nixpds,where Γτ is the stress boundary and τ¯i is the surface traction associated with node i.Calculating the nodal acceleration, the nodal velocities can be updated, and after mapping from nodes to particles, the acceleration and density of particles can be updated and combined with the constitutive equation and EOS:(7)ρn+1=ρn1+trΔεpn,ρn+1=ρn1+trΔεpn.After these steps, the now-deformed background grid is discarded and the original grid is used for the next time step [52, 55, 57].To begin the next time step, the nodal velocities must be calculated by extrapolation from the material points to the grid, and it can be expressed as(8)∑j=1nmijvj=∑p=1npmpvpNi,pxp.


## 2.1. Governing Equations

In the MPM, the background grid is divided to calculate the momentum equation, and the material domainΩ is discretized into a collection of particles, as shown in Figure 1.Figure 1
Discrete schematic diagram of the MPM.For continuous materials, they are governed by the conservation of momentum equations:(1)ρdvdt=∇⋅σ+ρb,where, ρ=ρx,t is the particle density, v=vx,t is the velocity vector, σ=σx,t is the stress tensor, and b=bx,t is the body force.The mass conservation equation is(2)dρdt+ρ∇⋅v=0,and the conservation of energy equation is(3)ρdedt=σ:dεdt+ρv⋅b,where ε=εx,t is the strain tensor and dε/dt is the strain rate tensor.


## 2.2. MPM Solution Scheme

The first step in the single time step of the MPM is to extrapolate the particle momenta and masses to the background. The equations are(4)pi=∑p=1npmpvpNi,p,mij=∑p=1npmpNi,p,where pi is the nodal momentum, vp is the particle velocity, mp is the particle mass, Ni,p is the shape function for node i evaluated at the current position of particle p, and mij is the nodal mass in the lumped (or diagonal) mass matrix.In this representation, the acceleration vector at each grid point is found by(5)∑j=1nmijaj=fiint+fiext,where fiint and fiext are internal and external force vectors associated with node i, and they are defined by(6)fiint=−∑p=1npmpρpσxpt,t⋅∇Nixp,fiext=∑p=1npmpNixpbxp,t+∫Γττ¯xp,t⋅Nixpds,where Γτ is the stress boundary and τ¯i is the surface traction associated with node i.Calculating the nodal acceleration, the nodal velocities can be updated, and after mapping from nodes to particles, the acceleration and density of particles can be updated and combined with the constitutive equation and EOS:(7)ρn+1=ρn1+trΔεpn,ρn+1=ρn1+trΔεpn.After these steps, the now-deformed background grid is discarded and the original grid is used for the next time step [52, 55, 57].To begin the next time step, the nodal velocities must be calculated by extrapolation from the material points to the grid, and it can be expressed as(8)∑j=1nmijvj=∑p=1npmpvpNi,pxp.


## 3. Axisymmetric Problems Solution and Dynamic Background Domain

As the MPM has matured, various extensions have been developed. The two core issues in this paper are that the material point method, which introduces the mirror reflection boundary condition, is used to calculate the axisymmetric problems and the concept of “dynamic background domain” is put forward to calculate the large deformation problems. It is worth noting that when calculating the axisymmetric problems useing the MPM with a mirror reflection boundary condition, the rectangular coordinate system, instead of the column coordinate system mentioned in [36], was used.


### 3.1. Mirror Reflection Boundary Condition

Many interesting scenarios lend themselves to being treated as axisymmetric. These include, for example, rod penetration, explosion, and impact events. Many scholars have also conducted research on formulation for an axisymmetric MPM. The particles were treated as point masses with the axisymmetric formulation presented by Sulsky and Schreyer [52], and they use column coordinate systems to deal with the Taylor impact problems. A new shape function and gradients specific for axisymmetry were presented by Nairn and Guilkey [36]; the transformation and extensions of shape function included traction boundary conditions, thermal conduction, explicit cracks, solvent diffusion, and convected particle-domain integration (CPDI).When solving the specific axisymmetric problem, it is only necessary to analyze the motion of the research object on the plane containing the symmetric axes [36]. However, we want to use rectangular coordinate system to solve the explosion problems, and the simple application of traditional symmetric boundaries cannot be used to solve those axisymmetric problems within the framework of the MPM. In the specific calculation process, we need to solve the following two problems:(1)
The influence domain of a single particle is not limited to the initial background grid in which the particle is located, so it is necessary to establish the virtual particles of multiple rows of meshes when solving the axisymmetric problem by using MPM. In the theory of the GIMP (Generalized Interpolation Material Point Method), each material point affects the nodes of its upper and lower three-row grids during the first mapping process. Meanwhile, each material point affects the nodes of the only grid which is located during the first mapping process based on the theory of the standard MPM. In order to broaden the applicability of the “Mirror reflection boundary condition,” we choose to generate three rows of virtual grids and generate corresponding virtual particles within the virtual grid on the other side of the symmetric axis, as shown in Figure2.(2)
After the arrangement of virtual mesh and virtual particles is completed, the constraints that the particles should meet need to be considered. According to the characteristics of axisymmetric problems, the physical information on the two particles relative to the symmetrical axis should be identical. At the end of each time step, considering the calculation process of the material point method in a single time step, the displacement and velocity of the virtual particles should be related to the displacement and velocity of the corresponding particles: in the direction parallel to the symmetrical axis, the size is equal and the direction is the same. In the direction perpendicular to the symmetrical axis, the size is equal and the direction is opposite. It can be expressed as(9)xx′=xx,xy′=−xy,vx′=vx,vy′=−vy.Figure 2
Virtual grid schematic.At this point, after the physical information is mapped from the particles to the grid node in the next time step, the physical information, such as the quality of the virtual nodes and the momentum of the virtual particles, are exactly the same as the physical information of its corresponding nodes and particles, and the calculated results will be identical, so that the axisymmetric problems can be resolved with the MPM introducing the boundary conditions; meanwhile, the number of the material point can be greatly reduced during the process of calculation. The source code for this section is detailed in Algorithm1.


### 3.2. Dynamic Background Domain (DBD)

According to the basic theory of the material point method, there are two mapping processes between the particles and the grids in a single time step; the calculation area needs to be larger than the arrangement area of particles in case of the large deformation problems, such as explosion and impact engineering problems. At the beginning of the numerical calculation, the transfer of the pressure, temperature, and other physical indicators of the detonation product exists only in the small area near the energy-containing material, and the physical information of all particles in the calculated area needs to be updated at each time step when the conventional material point method is used. If we can determine the influenced area of the current time step and only update the particles in this area, the calculation time can be shortened and the computational efficiency of the numerical method can be improved under the premise of ensuring the accuracy of the calculation results.The key to the method is to search for the affected domain in the current time step. Although particles carry a lot of physical information in a single time step, the most primitive physical quantities are velocity and pressure when the feature values of particles are updated. The velocity and pressure are given by(10)u˙pn+1=u˙pn+u¨pΔt,σ˙p=T:ε˙p,where T is a four-order tensor.Therefore, we can treat the velocity and pressure of one particle as a sign of whether it reacts or not, that is, whether it is in the affected domain or not. When the speed or pressure is not the same as the initial value, it can be determined to be “active” particles. After getting the coordinate extremes of these “active” points, the affected domain is determined by the coordinate extremes, which can contain all the “active” particles. After the affected domain is determined, we only need to update the physical information of the “active” particles and the affected domain.Also, for the sake of insurance, we can still expand the domain identified by “active” particles into other three meshes, thus treating them as a new affected domain.At this point, we only need to update the physical information for the particles which are in the affected domain.The affected domain is called the “dynamic background domain, DBD,” and the specific solution steps for dynamic background domains are as follows:(1)
Before the calculation begins, we define the “Zone Properties” of all particles as “1”(2)
At the timet, we determine whether physical information, such as pressure, temperature, and velocity, of a single particle is an initial value; if yes, we go to Type. 5; if not, we go to Type. 3(3)
All the particles are treated as “active” particles, and the maximum and minimum values of coordinates of such “active” particles are determined(4)
The boundary determined by the extremum by three meshes, respectively is expanded, to determine the dynamic background domain(5)
We determine if the material point is in the dynamic background domain; if yes, we go to Type 6; if not, we go to Type 9(6)
We update the zone properties of particles as “0”(7)
We calculate the particles by the MPM, whose zone property is “0”(8)
We determine the time incrementdt(9)
We update the time step by “t=t + dt” and determine whether the defined calculation time is reached; if yes, we finish this calculation process; if not, we go to Type 3 again, until the end of the calculation programIn Algorithm2, the source code of the abovementioned expressions is given, and the flowchart for this processing method is shown in Figure 3.Figure 3
Processing methods of the dynamic background domain.


## 3.1. Mirror Reflection Boundary Condition

Many interesting scenarios lend themselves to being treated as axisymmetric. These include, for example, rod penetration, explosion, and impact events. Many scholars have also conducted research on formulation for an axisymmetric MPM. The particles were treated as point masses with the axisymmetric formulation presented by Sulsky and Schreyer [52], and they use column coordinate systems to deal with the Taylor impact problems. A new shape function and gradients specific for axisymmetry were presented by Nairn and Guilkey [36]; the transformation and extensions of shape function included traction boundary conditions, thermal conduction, explicit cracks, solvent diffusion, and convected particle-domain integration (CPDI).When solving the specific axisymmetric problem, it is only necessary to analyze the motion of the research object on the plane containing the symmetric axes [36]. However, we want to use rectangular coordinate system to solve the explosion problems, and the simple application of traditional symmetric boundaries cannot be used to solve those axisymmetric problems within the framework of the MPM. In the specific calculation process, we need to solve the following two problems:(1)
The influence domain of a single particle is not limited to the initial background grid in which the particle is located, so it is necessary to establish the virtual particles of multiple rows of meshes when solving the axisymmetric problem by using MPM. In the theory of the GIMP (Generalized Interpolation Material Point Method), each material point affects the nodes of its upper and lower three-row grids during the first mapping process. Meanwhile, each material point affects the nodes of the only grid which is located during the first mapping process based on the theory of the standard MPM. In order to broaden the applicability of the “Mirror reflection boundary condition,” we choose to generate three rows of virtual grids and generate corresponding virtual particles within the virtual grid on the other side of the symmetric axis, as shown in Figure2.(2)
After the arrangement of virtual mesh and virtual particles is completed, the constraints that the particles should meet need to be considered. According to the characteristics of axisymmetric problems, the physical information on the two particles relative to the symmetrical axis should be identical. At the end of each time step, considering the calculation process of the material point method in a single time step, the displacement and velocity of the virtual particles should be related to the displacement and velocity of the corresponding particles: in the direction parallel to the symmetrical axis, the size is equal and the direction is the same. In the direction perpendicular to the symmetrical axis, the size is equal and the direction is opposite. It can be expressed as(9)xx′=xx,xy′=−xy,vx′=vx,vy′=−vy.Figure 2
Virtual grid schematic.At this point, after the physical information is mapped from the particles to the grid node in the next time step, the physical information, such as the quality of the virtual nodes and the momentum of the virtual particles, are exactly the same as the physical information of its corresponding nodes and particles, and the calculated results will be identical, so that the axisymmetric problems can be resolved with the MPM introducing the boundary conditions; meanwhile, the number of the material point can be greatly reduced during the process of calculation. The source code for this section is detailed in Algorithm1.


## 3.2. Dynamic Background Domain (DBD)

According to the basic theory of the material point method, there are two mapping processes between the particles and the grids in a single time step; the calculation area needs to be larger than the arrangement area of particles in case of the large deformation problems, such as explosion and impact engineering problems. At the beginning of the numerical calculation, the transfer of the pressure, temperature, and other physical indicators of the detonation product exists only in the small area near the energy-containing material, and the physical information of all particles in the calculated area needs to be updated at each time step when the conventional material point method is used. If we can determine the influenced area of the current time step and only update the particles in this area, the calculation time can be shortened and the computational efficiency of the numerical method can be improved under the premise of ensuring the accuracy of the calculation results.The key to the method is to search for the affected domain in the current time step. Although particles carry a lot of physical information in a single time step, the most primitive physical quantities are velocity and pressure when the feature values of particles are updated. The velocity and pressure are given by(10)u˙pn+1=u˙pn+u¨pΔt,σ˙p=T:ε˙p,where T is a four-order tensor.Therefore, we can treat the velocity and pressure of one particle as a sign of whether it reacts or not, that is, whether it is in the affected domain or not. When the speed or pressure is not the same as the initial value, it can be determined to be “active” particles. After getting the coordinate extremes of these “active” points, the affected domain is determined by the coordinate extremes, which can contain all the “active” particles. After the affected domain is determined, we only need to update the physical information of the “active” particles and the affected domain.Also, for the sake of insurance, we can still expand the domain identified by “active” particles into other three meshes, thus treating them as a new affected domain.At this point, we only need to update the physical information for the particles which are in the affected domain.The affected domain is called the “dynamic background domain, DBD,” and the specific solution steps for dynamic background domains are as follows:(1)
Before the calculation begins, we define the “Zone Properties” of all particles as “1”(2)
At the timet, we determine whether physical information, such as pressure, temperature, and velocity, of a single particle is an initial value; if yes, we go to Type. 5; if not, we go to Type. 3(3)
All the particles are treated as “active” particles, and the maximum and minimum values of coordinates of such “active” particles are determined(4)
The boundary determined by the extremum by three meshes, respectively is expanded, to determine the dynamic background domain(5)
We determine if the material point is in the dynamic background domain; if yes, we go to Type 6; if not, we go to Type 9(6)
We update the zone properties of particles as “0”(7)
We calculate the particles by the MPM, whose zone property is “0”(8)
We determine the time incrementdt(9)
We update the time step by “t=t + dt” and determine whether the defined calculation time is reached; if yes, we finish this calculation process; if not, we go to Type 3 again, until the end of the calculation programIn Algorithm2, the source code of the abovementioned expressions is given, and the flowchart for this processing method is shown in Figure 3.Figure 3
Processing methods of the dynamic background domain.


## 4. Numerical Example

Two numerical examples based on steel under explosion and impact problems are presented, in which these examples involve not only large deformations but also multiple-field coupling. In the first example, both the original MPM and the MPM with a mirror reflection boundary condition are used to calculate an axisymmetric problem. In the second example, the original MPM and MPM with DBD are compared to solve the large deformation problem. In both examples, numerical simulations are all occurring in the water medium, and the equation descriptions and the determined values of the material are employed according to [57]. All the simulations presented in this section are performed with a FORTRAN-based material point method. The calculations are processed on a personal computer with an Intel core Q6600 processor and 4 GB of RAM. The calculated results are described as follows.


### 4.1. Mirror Reflection Boundary Condition (Axisymmetric Problem)

In this numerical simulation, the explosive impact dynamic response of explosives in waters with approximate infinite size is numerically simulated, and the explosives are still TNT  materials. In the symmetrical surface of the three-dimensional problem, the computational domain is full of water with 2 m × 2  m dimension. TNT is a cuboid of 0.02  m × 0.02  m dimension and located in the middle of the domain. A thick steel plate of 0.01  m × 2  m dimension is mounted on one side of the TNT, with the clamped-edge condition assumed. The model of the axisymmetric problem takes the upper part of the model of the 3D problem, which is shown in Figure4.Figure 4
Model of the axisymmetric problem. (a) Model of the 3D problem (symmetrical surface). (b) Model of the axisymmetric problem.
(a)(b)In the process of calculation, the grid size is 0.005  m × 0.005  m, material point size is 0.0025  m × 0.0025  m, and calculation time is 220μs. The detonation product is described by the JWL state equation [58]:(11)p=A1−wR1Ve−R1V+B1−wR2Ve−R2V+wEV,where A=3.738 × 105 MPa, B=3.75 × 103 MPa, R1 = 4.15, R2 = 0.9, and w = 0.35.The water medium and steel plate are described by the Mie–Grüneisen equation of state (EOS) [58]:(12)p=ρ0C02μ+γ0+aμe,μ<0,ρ0C2μ1+1−γ0/2μ−a/2μ2/1−S1−1μ−S2μ2/μ+1−S3μ3/μ+12,μ>0,where the parameter values for water and steel plates are shown in Table 1.Table 1
Mie–Grüneisen equation of state parameters for water and steel plates.
ParametersC (m/s)S1S2S3γ0ae0Value for water14802.56−1.9860.22680.500.000.00Value for steel46001.490.000.002.170.460.00The stress-strain relationship of steel is described by the Johnson–Cook constitutive model [58]:(13)σy=A+Bεpn1+Clnε˙p∗1−T∗m,where the parameter values for water and steel plates are shown in Table 2.Table 2
Johnson–Cook material model parameters for steel.
ParametersA (MPa)B (MPa)C0nmTm/KTr/KValue7925100.0140.261.031793293The damage failure of steel is described by the Johnson–Cook failure model [58].(14)D=∑Δεpεf,where εf is the destroying strain of materials, which can be described by(15)εf=D1+D2expD3σ∗×1+D4lnε˙p∗×1+D5lnT∗,where the values of D1–D5 are 0.05, 3.44, -2.12, 0.002, and 1.61.Particles of diverse positions and materials are selected as observation points, and the pressure and internal energy of the observation points are recorded. The detonation point is located in the center of TNT and takes the detonation point coordinates as (0,0), and the eigenvalues of the observation points are shown in Table3. Pressure-time curves and energy-time curves of the gauges are shown in Figures 5 and 6. Pressure information from different locations is also recorded. The coordinates of locations are shown in Table 4, and the pressure-time curves are shown in Figure 7.Table 3
Information of Lagrange particles.
NumberCoordinate value (cm)NotesGauge 1(0,0)Ignition pointGauge 2(1,0)/Gauge 3(5,0)/Gauge 4(10,0)/Gauge 5(0,5)/Gauge 6(1,5)/Gauge 7(5,5)/Gauge 8(10,5)/Figure 5
Pressure-time curves of Lagrange particles. (a) Pressure-time curve of gauge 1. (b) Pressure-time curve of gauge 2. (c) Pressure-time curve of gauge 3. (d) Pressure-time curve of gauge 4. (e) Pressure-time curve of gauge 5. (f) Pressure-time curve of gauge 6. (g) Pressure-time curve of gauge 7. (h) Pressure-time curve of gauge 8.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 6
Energy-time curves of Lagrange particles. (a) Energy-time curve of gauge 1. (b) Energy-time curve of gauge 2. (c) Energy-time curve of gauge 3. (d) Energy-time curve of gauge 4. (e) Energy-time curve of gauge 5. (f) Energy-time curve of gauge 6. (g) Energy-time curve of gauge 7. (h) Energy-time curve of gauge 8.
(a)(b)(c)(d)(e)(f)(g)(h)Table 4
Information of Euler coordinates.
NumberCoordinate value (cm)NotesPosition 1(1.125,0.125)/Position 2(1.875,0.125)/Position 3(1.125,5.125)/Position 4(1.875,5.125)/Figure 7
Pressure-time curves of the Euler coordinate. (a) Pressure-time curve of position 1. (b) Pressure-time curve of position 2. (c) Pressure-time curve of position 3. (d) Pressure-time curve of position 4.
(a)(b)(c)(d)The results of the 3D problem and axisymmetric problem are basically exactly the same, and the peak and changing trends of pressure at the observation point are exactly the same. Because the explicit integration algorithm was used during the process of calculation, the updated time step will be different from 3D problems and 2D problems. As the calculation progresses, the results of the pressure-time curve and energy-time curve will also be different. The pressure contours of the whole model at different typical times are also shown in Figures8–11.Figure 8
Pressure contours at 10μs. (a) Pressure contour of the 3D problem. (b) Pressure contour of the axisymmetric problem.
(a)(b)Figure 9
Pressure contours at 50μs. (a) Pressure contour of the 3D problem. (b) Pressure contour of the axisymmetric problem.
(a)(b)Figure 10
Pressure contours at 100μs. (a) Pressure contour of the 3D problem. (b) Pressure contour of the axisymmetric problem.
(a)(b)Figure 11
Pressure contours at 220μs. (a) Pressure contour of the 3D problem. (b) Pressure contour of the axisymmetric problem.
(a)(b)Table5 records the computing CPU time of the original 3D MPM and MPM with a mirror reflection boundary. The original MPM costs about 1560s when it is used to solve this 3D problem, while the MPM with a mirror reflection boundary costs about 208s when it is used to solve the same problem, which reduces the computing time significantly.Table 5
Comparison of CPU time.
MethodCPU time/s3D MPM1560.38MPM with a mirror reflection boundary208.17


### 4.2. Dynamic Background Domain (Explosion and Impact Dynamic Response Problem)

The problem of explosion impact dynamic response of TNT  explosive to steel plate in a water medium is numerically simulated. The computational domain is full of water with 1 m  ×  1  m × 1m dimension. TNT is a cuboid of 0.02  m  ×  0.02  m × 1m dimension and located in the middle of the domain. A thick steel plate of 0.01  m  ×  1 m × 1 m dimension is mounted on one side of the TNT, with the clamped edge condition assumed. Because of symmetry, the analogy model can be further simplified as a 2D planar model, as shown in Figure12.Figure 12
Model of explosion and impact dynamic response problem.In the process of calculation, the grid size is 0.005  m  ×  0.005  m, material point size is 0.0025  m  ×  0.0025  m, and the calculation time is 200μs.The layout of coordinate systems, constitutive equations, and equation of state are similar with example 1, and the eigenvalues of the observation points are shown in Table6.Table 6
Viewpoint eigenvalues of Section4.2.
NumberCoordinate value (cm)NotesGauge 1(0,0)Ignition pointGauge 2(1,0)/The results of pressure-time curves are shown in Figure13 and Table 7. Also, the pressure-time change relationship at the same observation point is basically the same. The pressure contours of the whole model are shown in Figure 14. The result shows the characteristics of large deformation and multifield coupling, and to these two results, the radius of the detonation wave is about 43 cm and vacuum areas with a radius of about 6 cm are both formed in the center area of the explosive.Figure 13
The results of pressure at viewpoints. (a) Pressure-time curve of gauge 1. (b) Pressure-time curve of gauge 2.
(a)(b)Table 7
The peak pressure of viewpoints in Section4.2.
ViewpointsMethodsPeak pressure (Mbar)Error value (Mbar)Relative errors (%)Gauge 1MPM with DBD0.0690.0045.48MPM without DBD0.073Gauge 2MPM with DBD0.0600.0069.09MPM without DBD0.066Figure 14
Pressure contours at 200μs. (a) MPM with DBD (50 μs). (b) MPM without DBD (50 μs). (c) MPM with DBD (100 μs). (d) MPM without DBD (100 μs). (e) MPM with DBD (200 μs). (f) MPM without DBD (200 μs).
(a)(b)(c)(d)(e)(f)Table8 records the computing time of the MPM with the dynamic background domain and the original MPM without the dynamic background domain. The original method takes 712.69s, and the proposed method takes 182.51s. The speed is a quarter of times faster than that of the original MPM simulation, and the computation efficiency of the method with DBD is greatly increased. Hence, the MPM with DBD makes it more efficient and convenient to solve the large deformation and large-scale structure.Table 8
Comparison of CPU time.
MethodCPU time/sMPM without DBD712.69MPM with DBD182.51


## 4.1. Mirror Reflection Boundary Condition (Axisymmetric Problem)

In this numerical simulation, the explosive impact dynamic response of explosives in waters with approximate infinite size is numerically simulated, and the explosives are still TNT  materials. In the symmetrical surface of the three-dimensional problem, the computational domain is full of water with 2 m × 2  m dimension. TNT is a cuboid of 0.02  m × 0.02  m dimension and located in the middle of the domain. A thick steel plate of 0.01  m × 2  m dimension is mounted on one side of the TNT, with the clamped-edge condition assumed. The model of the axisymmetric problem takes the upper part of the model of the 3D problem, which is shown in Figure4.Figure 4
Model of the axisymmetric problem. (a) Model of the 3D problem (symmetrical surface). (b) Model of the axisymmetric problem.
(a)(b)In the process of calculation, the grid size is 0.005  m × 0.005  m, material point size is 0.0025  m × 0.0025  m, and calculation time is 220μs. The detonation product is described by the JWL state equation [58]:(11)p=A1−wR1Ve−R1V+B1−wR2Ve−R2V+wEV,where A=3.738 × 105 MPa, B=3.75 × 103 MPa, R1 = 4.15, R2 = 0.9, and w = 0.35.The water medium and steel plate are described by the Mie–Grüneisen equation of state (EOS) [58]:(12)p=ρ0C02μ+γ0+aμe,μ<0,ρ0C2μ1+1−γ0/2μ−a/2μ2/1−S1−1μ−S2μ2/μ+1−S3μ3/μ+12,μ>0,where the parameter values for water and steel plates are shown in Table 1.Table 1
Mie–Grüneisen equation of state parameters for water and steel plates.
ParametersC (m/s)S1S2S3γ0ae0Value for water14802.56−1.9860.22680.500.000.00Value for steel46001.490.000.002.170.460.00The stress-strain relationship of steel is described by the Johnson–Cook constitutive model [58]:(13)σy=A+Bεpn1+Clnε˙p∗1−T∗m,where the parameter values for water and steel plates are shown in Table 2.Table 2
Johnson–Cook material model parameters for steel.
ParametersA (MPa)B (MPa)C0nmTm/KTr/KValue7925100.0140.261.031793293The damage failure of steel is described by the Johnson–Cook failure model [58].(14)D=∑Δεpεf,where εf is the destroying strain of materials, which can be described by(15)εf=D1+D2expD3σ∗×1+D4lnε˙p∗×1+D5lnT∗,where the values of D1–D5 are 0.05, 3.44, -2.12, 0.002, and 1.61.Particles of diverse positions and materials are selected as observation points, and the pressure and internal energy of the observation points are recorded. The detonation point is located in the center of TNT and takes the detonation point coordinates as (0,0), and the eigenvalues of the observation points are shown in Table3. Pressure-time curves and energy-time curves of the gauges are shown in Figures 5 and 6. Pressure information from different locations is also recorded. The coordinates of locations are shown in Table 4, and the pressure-time curves are shown in Figure 7.Table 3
Information of Lagrange particles.
NumberCoordinate value (cm)NotesGauge 1(0,0)Ignition pointGauge 2(1,0)/Gauge 3(5,0)/Gauge 4(10,0)/Gauge 5(0,5)/Gauge 6(1,5)/Gauge 7(5,5)/Gauge 8(10,5)/Figure 5
Pressure-time curves of Lagrange particles. (a) Pressure-time curve of gauge 1. (b) Pressure-time curve of gauge 2. (c) Pressure-time curve of gauge 3. (d) Pressure-time curve of gauge 4. (e) Pressure-time curve of gauge 5. (f) Pressure-time curve of gauge 6. (g) Pressure-time curve of gauge 7. (h) Pressure-time curve of gauge 8.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 6
Energy-time curves of Lagrange particles. (a) Energy-time curve of gauge 1. (b) Energy-time curve of gauge 2. (c) Energy-time curve of gauge 3. (d) Energy-time curve of gauge 4. (e) Energy-time curve of gauge 5. (f) Energy-time curve of gauge 6. (g) Energy-time curve of gauge 7. (h) Energy-time curve of gauge 8.
(a)(b)(c)(d)(e)(f)(g)(h)Table 4
Information of Euler coordinates.
NumberCoordinate value (cm)NotesPosition 1(1.125,0.125)/Position 2(1.875,0.125)/Position 3(1.125,5.125)/Position 4(1.875,5.125)/Figure 7
Pressure-time curves of the Euler coordinate. (a) Pressure-time curve of position 1. (b) Pressure-time curve of position 2. (c) Pressure-time curve of position 3. (d) Pressure-time curve of position 4.
(a)(b)(c)(d)The results of the 3D problem and axisymmetric problem are basically exactly the same, and the peak and changing trends of pressure at the observation point are exactly the same. Because the explicit integration algorithm was used during the process of calculation, the updated time step will be different from 3D problems and 2D problems. As the calculation progresses, the results of the pressure-time curve and energy-time curve will also be different. The pressure contours of the whole model at different typical times are also shown in Figures8–11.Figure 8
Pressure contours at 10μs. (a) Pressure contour of the 3D problem. (b) Pressure contour of the axisymmetric problem.
(a)(b)Figure 9
Pressure contours at 50μs. (a) Pressure contour of the 3D problem. (b) Pressure contour of the axisymmetric problem.
(a)(b)Figure 10
Pressure contours at 100μs. (a) Pressure contour of the 3D problem. (b) Pressure contour of the axisymmetric problem.
(a)(b)Figure 11
Pressure contours at 220μs. (a) Pressure contour of the 3D problem. (b) Pressure contour of the axisymmetric problem.
(a)(b)Table5 records the computing CPU time of the original 3D MPM and MPM with a mirror reflection boundary. The original MPM costs about 1560s when it is used to solve this 3D problem, while the MPM with a mirror reflection boundary costs about 208s when it is used to solve the same problem, which reduces the computing time significantly.Table 5
Comparison of CPU time.
MethodCPU time/s3D MPM1560.38MPM with a mirror reflection boundary208.17


## 4.2. Dynamic Background Domain (Explosion and Impact Dynamic Response Problem)

The problem of explosion impact dynamic response of TNT  explosive to steel plate in a water medium is numerically simulated. The computational domain is full of water with 1 m  ×  1  m × 1m dimension. TNT is a cuboid of 0.02  m  ×  0.02  m × 1m dimension and located in the middle of the domain. A thick steel plate of 0.01  m  ×  1 m × 1 m dimension is mounted on one side of the TNT, with the clamped edge condition assumed. Because of symmetry, the analogy model can be further simplified as a 2D planar model, as shown in Figure12.Figure 12
Model of explosion and impact dynamic response problem.In the process of calculation, the grid size is 0.005  m  ×  0.005  m, material point size is 0.0025  m  ×  0.0025  m, and the calculation time is 200μs.The layout of coordinate systems, constitutive equations, and equation of state are similar with example 1, and the eigenvalues of the observation points are shown in Table6.Table 6
Viewpoint eigenvalues of Section4.2.
NumberCoordinate value (cm)NotesGauge 1(0,0)Ignition pointGauge 2(1,0)/The results of pressure-time curves are shown in Figure13 and Table 7. Also, the pressure-time change relationship at the same observation point is basically the same. The pressure contours of the whole model are shown in Figure 14. The result shows the characteristics of large deformation and multifield coupling, and to these two results, the radius of the detonation wave is about 43 cm and vacuum areas with a radius of about 6 cm are both formed in the center area of the explosive.Figure 13
The results of pressure at viewpoints. (a) Pressure-time curve of gauge 1. (b) Pressure-time curve of gauge 2.
(a)(b)Table 7
The peak pressure of viewpoints in Section4.2.
ViewpointsMethodsPeak pressure (Mbar)Error value (Mbar)Relative errors (%)Gauge 1MPM with DBD0.0690.0045.48MPM without DBD0.073Gauge 2MPM with DBD0.0600.0069.09MPM without DBD0.066Figure 14
Pressure contours at 200μs. (a) MPM with DBD (50 μs). (b) MPM without DBD (50 μs). (c) MPM with DBD (100 μs). (d) MPM without DBD (100 μs). (e) MPM with DBD (200 μs). (f) MPM without DBD (200 μs).
(a)(b)(c)(d)(e)(f)Table8 records the computing time of the MPM with the dynamic background domain and the original MPM without the dynamic background domain. The original method takes 712.69s, and the proposed method takes 182.51s. The speed is a quarter of times faster than that of the original MPM simulation, and the computation efficiency of the method with DBD is greatly increased. Hence, the MPM with DBD makes it more efficient and convenient to solve the large deformation and large-scale structure.Table 8
Comparison of CPU time.
MethodCPU time/sMPM without DBD712.69MPM with DBD182.51


## 5. Conclusions

With the progress of science and technology and the development of computer technology, it is gradually becoming a trend to solve the problems of large deformation or fluid-solid coupling such as explosion and penetration by means of numerical simulation. When the Lagrangian algorithm is used to calculate the large deformation problems such as explosion detonation, there is a shortcoming of mesh distortion, which will make a great error in the calculation results, and it is difficult to calculate the fluid-solid coupling problem by the Euler algorithm, such as the material boundary and the historical deformation of the tracking material. As one of the nongrid algorithms, the material point method has high computational efficiency, and there are no mesh distortion problems, and the coupling conditions are automatically satisfied, so it has outstanding advantages in dealing with the problems of large deformation and fluid-solid coupling. On the basis of the theory of the material point method and aiming at the problem of large deformation, this paper puts forward the concept of the dynamic background domain, introduces the mirror reflection boundary condition, studies the axisymmetric problem, uses Fortran 95 language to write a set of calculation programs suitable for the axisymmetric problem, and compares and validates the concrete examples. On the premise of ensuring the accuracy of its calculation, the method with a mirror reflection boundary condition can save about 85% of the calculation time when it is used to solve axisymmetric problems; to solve the large deformation problem mentioned in this paper, the method with DBD can also save about 75% of the calculation time.(1)
In view of the axisymmetric problem, the introduction of specular reflection boundary conditions can ensure that the three-dimensional axisymmetric problem is reduced to a two-dimensional problem without changing the original right-angled coordinate system, and the calculation efficiency is improved on the premise of ensuring the accuracy of the results.(2)
This paper analyzes the characteristics of the problem of large deformation, puts forward the concept of dynamic background domain, calculates the influence area of a large deformation problem, updates only the particle information in the affected area, and improves the calculation efficiency under the premise of ensuring the accuracy of the calculation results.



---

*Source: 8854318-2021-02-13.xml*