Impact force analysis using the B-splinematerial point method

Use your smartphone to scan this QR code and download this article ABSTRACT In the MPM algorithm, all the particles are formulated in a single-valued velocity field hence the non-slip contact can be satisfied without any contact treatment. However, in some impact and penetration problems, the non-slip contact condition is not appropriate andmay even yield unreasonable results, so it is important to overcome this drawback by using a contact algorithm in the MPM. In this paper, the variation of contact force with respect to time caused by the impact is investigated. The MPM using the Lagrange basis function, so causing the cell-crossing phenomenon when a particle moves from one cell to another. The essence of this phenomenon is due to the discontinuity of the gradient of the linear basis function. The accuracy of the results is therefore also affected. The high order B-spline MPM is used in this study to overcome the cell-crossing error. The BSMPMuses higher-order B-spline functions tomake sure the derivatives of the shape functions are continuous, so that alleviate the error. The algorithm of MPM and BSMPM has some differences in defining the computational grid. Hence, the original contact algorithm in MPM needs to be modified to be suitable in order to use in the BSMPM. The purpose of this study is to construct a suitable contact algorithm for BSMPM and then use it to investigate the contact force caused by impact. Some numerical examples are presented in this paper, the impact of two circular elastic disks and the impact of a soft circular disk into a stiffer rectangular block. All the results of contact force obtained from this study are compared with finite element results and perform a good agreement, the energy conservation is also considered.


INTRODUCTION
The material point method (MPM) was first developed in 1994 by Sulsky and his colleagues 1 . Over 25 years of development, the number of researchers working on it is increasing more and more. Many universities and institutes around the world have investigated this method, such as Delft University of Technology 2 , Stuttgart University 3 , Cardiff University 4 . The MPM uses both Lagrangian description and Eulerian description 1 so it has the advantages of both descriptions. MPM has been widely used to simulate high-velocity problems such as impact 5 and explosion 6 , large deformation problems 7 , fracture 8 and also Fluid-Structure Interaction 9 . However, the original MPM has a major shortcoming that affects the simulation results. When a particle moves across a cell boundary, it will lead to numerical errors due to the discontinuity of the gradient of the basis functions 1 . This is called the "cellcrossing error" 2 . In order to alleviate the effect of this phenomenon, different methods were proposed. Bardenhagen et al. proposed the Generalized Interpolation Material Point method (GIMP) 10 . Variants in the GIMP branch were also introduced, Steffen et al. proposed the Uniform GIMP (uGIMP) 11 , the Convected Particle Domain Interpolation (CPDI) was introduced by Sadeghirad et al. 12 14 by applying the high order B-spline function into MPM algorithm. The BSMPM is then further improved by Tielen et al. 2 , Gan et al. 15 , Wobbes et al. 16 .
In the MPM algorithm, a single-valued velocity field is used for all particles so the non-slip contact condition between two bodies is satisfied automatically 1 . However, in some impact and penetration problems, the non-slip contact condition is not appropriate, so it is important to develop a contact algorithm for MPM. York et al. proposed a simple contact algorithm for MPM 17 , Bardenhagen et al. proposed an algorithm for multi-velocity field 18 , and many other improvements can be mentioned as Hu and Chen 19 , Huang et al. 20 , Nairn 21 , Ma et al. 22 .
This study using the BSMPM to mitigate the cellcrossing error. The BSMPM and MPM have differences in computational grid definition. Therefore, the contact algorithm for MPM cannot be directly applied to BSMPM. In this paper, the contact algorithm is modified to a suitable form to the BSMPM. The implementation steps are mentioned in Section 2.3. The contact force obtained from impact of two elastic objects are compared with the result from FEM, a slight difference between FEM and MPM (and BSMPM) results is observed and explained in Section 3.
The i-th B-spline basis function of order d (N i,d ) is defined by using Cox-de Boor recursion formula 15 . Firsly, the zeroth order (d=0) basis function must be defined the non-zero intervals [ξ i , ξ i+1 ) are called knot spans 2 . After obtaining N i,0 , higher order (d ≥ i ) basis functions are defined as the formula below in which the fraction 0/0 is assumed to be zero.
In two dimensions, the bivariate B-spline functions can be built from the tensor product of the univariate ones 8

B-spline Material Point Method
In 2D BSMPM, the computational domain is discretized by a parametric grid 15 . This grid is defined by two open knot vectors on two orthogonal directions Ξ = { ξ 1 , ξ 2 , ..., ξ n+p , ξ n+p+1 } and I = { η 1 , η 2 , ..., η m+q , η m+q+1 } as shown in Figure 2. The numbers of basis functions in ξ and η direction are n and m, respectively, so the total number of basis functions is n × m. A tensor product grid with the total of n × m nodes is constructed as shown in Figure 3, each node of this grid corresponds to one B-spline basis function as defined in Eq. (4). For example, the node with the position (1, 3) on the grid corresponds to the basis function N 1,3 (ξ , η) = N 1,p (ξ ) N 3,q (η). All the nodes on this tensor product grid are defined as control points in BSMPM (the same role for grid node in MPM), and in practice these control points are arbitrary distribution 15 . As shown in Figure 4 a second order (quadratic) BSMPM grid. The cell is made from 3 knot spans in x direction and 2 knot spans in y direction, so the number of knots in knot vectors are 4 and 3, respectively. The number of control points in x direction is 3 (knot spans) + 2 (order) = 5 and in y direction is equal 4. These control points play the role of grid nodes in the original MPM, the knots from knot vectors are only used for creating a computational grid. At can be seen in Figure 4, each cell has 9 control points, for example, the lower-left cell related to [1,2,3,6,7,8,11,12,13].
where (x min , y min ) is the lower-left control point and (x max , y max ) is the upper-right control point. This is the formula for mapping between the parameter space to the physical space. The derivatives of the B-spline basis functions are given as below 23 where J is the Jacobian matrix and defined by and the components are computed as where P denotes the coordinates of the control points and A is the global index of control point 23 .
In the BSMPM, for convenient the knot vector for an interval [0, L] is defined by where △x denotes the length of knot span 15 . And note that the knot vector must be normalized before a parametric grid is created, so the knot vector is rewritten as the following form , this is also a difference in the parameter space between B-spline basis functions and other functions. The control points are arbitrary distributed and they are in the physical coordinate (x, y).

Contact algorithm
This section presents the algorithm proposed by Bardenhagen et al. 18 and makes appropriate modification to apply into the BSMPM. When two bodies are approaching each other, there is a region where they have some of the same control points. These control points are viewed as the contact points, the contact algorithm is applied on these points only. In the contact region, the following equation 18 is used as a condition to check if two bodies are in contact or release ( where i denotes the i-th body in the computational domain, v cm I is the center-of-mass velocity 23 of the control point I-th for each pair in contact In Eq. (9), n i I is the normal vector of control point I-th of body i-th and computed as following steps. Firstly, the density ρ c for each cell in contact state is computed as below 23 where V e is volume of cell e-th, x c is the center of cell e-th. Remember that in the BSMPM each cell is made of knot spans (see Figure 4). In 1D, the function S x (x) is given by the following definition 23 The function S 2 (x, y) in Eq. (11) is obtained by multiplying two 1D functions S 2 (x, y) = S x (x) S y (y) . Finally, the normal vector of control point I-th is obtained 23 where G I is the derivatives of the B-spline basis functions. Before applying into Eq. (9) for checking contact, the normal vector in Eq. (13) must be normalized 23 The implementation of contact algorithm into the BSMPM algorithm can be summarized as following steps: Step 1: Mapping data from particles to control points 1. Compute the mass of I-th control point from the

RESULTS
Two numerical examples are presented in this section, particularly: • Collision of two circular disks. • Collision of a circular disk onto a rectangular block.
The first example investigates the contact of two circular surface with the same material. The second example studies the contact of a soft circular surface and hard flat surface.
To validate the results from these two examples, corresponding FEM models are created from ABAQUS software. FEM model is prepared with very fine mesh and set up with the same parameters and initial conditions as MPM (and BSMPM) model.

Collision of two circular disks
The problem is shown in Figure 5, two elastic disks with the same radius R = 0. The kinetic and strain energy obtained from BSMPM and FEM is shown in Figure 6. Kinetic energy in BSMPM decreases earlier than the result from FEM and strain energy in BSMPM increases earlier. This is reasonable for the contact in BSMPM algorithm and will be explained in the comment of Figure 7. The value of kinetic energy in both case are the same, while the strain energy in BSMPM is lower than FEM. Both case are in frictionless contact, so there is no energy  loss from friction, the strain energy loss in BSMPM is caused by other error factors. The variation of contact force during the impact process is shown in Figure 7. The FEM model used to simulate this problem has 3288 nodes. The results from MPM and BSMPM show that the impact of two bodies occurs earlier than the result in FEM as mentioned before. This is because the contact force in MPM is computed in the node of the computational grid (or control point in BSMPM), not in the particle of the body, so when two bodies approach the contact region and have the same control points, the contact is detected immediately although two bodies have not touched each other yet. In FEM, the contact is only detected when two bodies touch each other, so the contact force obtained in FEM is later than MPM. The contact force obtained from BSMPM using higher order B-spline functions also shows the smooth curve compared to the MPM and FEM. Figure 8 shows the von-Mises stress field during the impact process of two disks using the BSMPM. In detail, two disks approaching each other in Figure 8 (a), then two disks touch each other as shown in Figure 8 (b), two disks deform during the impact as shown in Figure 8 (c) and then bounce back in Figure 8 (d). After impact, two disks move far away as shown inFigure 8 (e).

Collision of a circular disk onto a rectangular block
In this example, a circular disk collides onto a stiffer rectangular block, as shown in Figure 9. The contact force obtained in this example also shows the similarity to the conclusions from the previous example. Figure 10 also shows that the impact occurs earlier in BSMPM, because BSMPM has more control points (nodes) than MPM so the contact is detected earlier. Similarly to the previous example, the contact force in BSMPM is smoother than the curve from  FEM and MPM. Figure 11 shows the collision of two objects, the von-Mises stress field and maximum stress field are presented.
To investigate the convergence of BSMPM and MPM, the computational domain with a set of 60 × 60 knot spans is retained. Different numbers of particles per cell (PPC) 4, 9 and 16 are analyzed. Figure 12 shows the total energy of the system respect to time. From the initial conditions, the total energy can be computed as ρπR 2 tv 2 /2 = 2.512 J and plotted by the black line in the figure. As shown in Figure 12, the case of MPM with PPC = 4 gives a very large deviation, and when PPC = 9, the result is significantly improved. In the case of BSMPM, there are no significant deviations and the results are slightly improved when increasing PPC.

DISCUSSIONS
As present in Section 3, there is a slight difference in the results of MPM, BSMPM and FEM. The mag-  nitude of the contact force in the three methods is the same, but the impact occurs earlier in MPM and BSMPM. This is because the contact force in MPM is computed in the node of the computational grid (or control point in BSMPM), not in the discrete particle of the object, so when two objects approach the contact region and have same the grid node (control points), the contact is detected immediately although two objects have not touched each other yet, there is still a gap. And because BSMPM uses higher-order shape functions, so the BSMPM has more control points (grid node) than MPM, the contact is therefore detected earlier. Moreover, if the contact algorithm is not used in MPM, the non-slip contact can still be determined automatically when two objects have the same grid node. In FEM, the contact is only detected when two objects touch each other (or even penetrate into each other), so the contact force obtained in FEM is later than MPM and BSMPM.

CONCLUSIONS
The contact algorithm has been successfully modified and applied into the BSMPM. The contact force obtained from this research is compared to FEM. A slight difference in the result is observed, this is because the contact force is calculated at the control point in the computational grid instead of the discrete