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

Based on the basic theory of the material point method (MPM), the factors aﬀecting the computational eﬃciency are analyzed and discussed, and the problem of improving calculation eﬃciency is studied. This paper introduces a mirror reﬂection boundary condition to the MPM to solve axisymmetric problems; to improve the computational eﬃciency 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 reﬂection boundary condition and MPM with DBD can improve the calculation eﬃciency 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.


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. e 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][13][14]. Rabczuk and Belytschko [15] proposed the cracking particles method (CPM), and it has been used to model the dynamic crack propagation [16][17][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][22][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]. e phase field modelling was also used to deal with the crack propagation [26][27][28][29]. A numerical manifold method was developed to deal with the Kirchhoff plate [30][31][32]. e 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][43][44], combustion reaction [45], cracks and fracture [46][47][48][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. ere 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][37][38].
Many interesting scenarios from the preceding categories lend themselves to being treated as axisymmetric. ese include, for example, impact and cratering events, rod penetration, and shaped charge jet formation. e 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. e 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. erefore 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][55][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. e algorithm, thereby, allows a more thorough search of all the materials point and seems to offer a greater ability to find global minima.
is paper is organized as follows. e algorithm of the explicit MPM is reviewed briefly in Section 2. e detailed formulation of the axisymmetric problems solution is presented in Section 3, which also describes the concept of dynamic background domain. e numerical simulations and several conclusions are given in Section 4. Finally, Section 5 draws a summary and some concluding remarks. e source codes of the mirror reflection boundary conditions and dynamic background domain are given in Algorithms 1 and 2.

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.
For continuous materials, they are governed by the conservation of momentum equations: is the stress tensor, and b � b(x, t) is the body force. e mass conservation equation is and the conservation of energy equation is where ε � ε(x, t) is the strain tensor and dε/dt is the strain rate tensor.

MPM Solution
Scheme. e first step in the single time step of the MPM is to extrapolate the particle momenta and masses to the background. e equations are where p i is the nodal momentum, v p is the particle velocity, m p is the particle mass, N i, p is the shape function for node i evaluated at the current position of particle p, and m ij is the nodal mass in the lumped (or diagonal) mass matrix.
In this representation, the acceleration vector at each grid point is found by where f int i and f ext i are internal and external force vectors associated with node i, and they are defined by module specular_reflection use param, only: parts, dy, partstotal, round, symmetry, reptotal implicit none contains open(1,file � "evaluate.dat") subroutine evaluate() implicit none integer i, j, k, t if (round �� 1) then (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 do end subroutine close (1) end module ALGORITHM 1: Source code for mirror reflection boundary condition. Mathematical Problems in Engineering 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: 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

Axisymmetric Problems Solution and Dynamic Background Domain
As the MPM has matured, various extensions have been developed. e 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.

Mirror Reflection Boundary Condition.
Many interesting scenarios lend themselves to being treated as axisymmetric. ese include, for example, rod penetration, explosion, and impact events. Many scholars have also conducted research on formulation for an axisymmetric MPM. e 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) e 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 deallocate(parts(i)%nodes,s(i)%pnode,s(i)%s,s(i)%sx) if(parts(i)%sign��0) then else if(parts(i)%x(1)<�zxmax.and.parts(i)%x(1)>�fxmax.and.parts(i)%x(2)<�zymax.and.parts(i)%x (2)  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 Figure 2. (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 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. e source code for this section is detailed in Algorithm 1.

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.
e 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. e velocity and pressure are given by where T is a four-order tensor. erefore, 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. e 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 time t, 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) e 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 increment dt Mathematical Problems in Engineering (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 program In Algorithm 2, the source code of the abovementioned expressions is given, and the flowchart for this processing method is shown in Figure 3.

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. e calculations are processed on a personal computer with an Intel core Q6600 processor and 4 GB of RAM. e calculated results are described as follows.

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. e model of the axisymmetric problem takes the upper part of the model of the 3D problem, which is shown in Figure 4.
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. e detonation product is described by the JWL state equation [58]: where A�3.738 × 10 5 MPa, B�3.75 × 10 3 MPa, R 1 � 4.15, R 2 � 0.9, and w � 0.35. e water medium and steel plate are described by the Mie-Grüneisen equation of state (EOS) [58]: where the parameter values for water and steel plates are shown in Table 1. e stress-strain relationship of steel is described by the Johnson-Cook constitutive model [58]: where the parameter values for water and steel plates are shown in Table 2. e damage failure of steel is described by the Johnson-Cook failure model [58].
where ε f is the destroying strain of materials, which can be described by where the values of D 1 -D 5 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. e 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 Table 3. 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. e coordinates of locations are shown in Table 4, and the pressure-time curves are shown in Figure 7.
e 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. e pressure contours of the whole model at different typical times are also shown in Figures 8-11.

Dynamic Background Domain (Explosion and Impact Dynamic Response Problem).
e problem of explosion impact dynamic response of TNT explosive to steel plate in a water medium is numerically simulated. e 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 Figure 12.
In  Table 6. e results of pressure-time curves are shown in Figure 13 and Table 7. Also, the pressure-time change relationship at the same observation point is basically the same. e pressure contours of the whole model are shown in Figure 14. e 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. Table 8           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.

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 fluidsolid 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) is 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.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.