Impact force analysis using the B-spline material point method

artphone to scan this penetration problems, the non-slip contact condition is not appropriate and may even yield unrea- QR code and download this article sonable 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 in- vestigated. 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 BSMPM uses higher-order B-spline functions to make sure the derivatives of the shape functions are 1Department of Engineering Mechanics, continuous, so that alleviate the error. The algorithm of MPM and BSMPM has some differences in Faculty of Applied Sciences, Ho Chi defining the computational grid. Hence, the original contact algorithm in MPM needs tobemodi- Minh City University of Technology, fied to be suitable in order to use in the BSMPM. The purpose of this study is to construct asuitable Vietnam. contact algorithm for BSMPM and then use it to investigate the contact force caused by impact. 2 Some numerical examples are presented in this paper, the impact of two circular elastic disks and Vietnam National University Ho Chi the impact of a soft circular disk into a stiffer rectangular block. All the results of contact force ob- Minh City, Vietnam. tained from this study are compared with finite element results and perform a good agreement, Correspondence the energy conservation is also considered. Key words: BSMPM, contact algorithm, contact force, impact, MPM Vay Siu Lo, Department of Engineering Mechanics, Faculty of Applied Sciences, Ho Chi Minh City University of Technology, Vietnam. the Convected Particle Domain Interpolation (CPDI) Vietnam National University Ho Chi Minh INTRODUCTION 12 City, Vietnam. The material point method (MPM) was first devel- was introduced by Sadeghirad et al. . Zhang et Email: losiuvay@hcmut.edu.vn oped in 1994 by Sulsky and his colleagues 1. Over al. modified the gradient of shape functions to en- hance the MPM 13. Steffen et al. introduced the B- Correspondence 25 years of development, the number of researchers spline MPM (BSMPM) 14 by applying the high order Thien Tich Truong, Department of working on it is increasing more and more. Many uni- Engineering Mechanics, Faculty of versities and institutes around the world have investi- B-spline function into MPM algorithm. The BSMPM Applied Sciences, Ho Chi Minh City gated this method, such as Delft University of Tech- is then further improved by Tielen et al. 2, Gan et al. 15, University of Technology, Vietnam. nology 2, Stuttgart University 3, Cardiff University 4. Wobbes et al. 16. Vietnam National University Ho Chi Minh City, Vietnam. The MPM uses both Lagrangian description and Eu- In the MPM algorithm, a single-valued velocity field is lerian description 1 so it has the advantages of both Email: tttruong@hcmut.edu.vn used for all particles so the non-slip contact condition descriptions. MPM has been widely used to simu- between two bodies is satisfied automatically 1. How- History late high-velocity problems such as impact 5 and ex- • Received: 18-11-2020 ever, in some impact and penetration problems, the plosion 6, large deformation problems 7, fracture 8 and • Accepted: 11-3-2021 non-slip contact condition is not appropriate, so it is also Fluid-Structure Interaction 9. • Published: 30-3-2021 important to develop a contact algorithm for MPM. However, the original MPM has a major shortcom- DOI : 10.32508/stdjet.v4i1.794 York et al. proposed a simple contact algorithm for ing that affects the simulation results. When a par- MPM 17, Bardenhagen et al. proposed an algorithm ticle moves across a cell boundary, it will lead to nu- for multi-velocity field 18, and many other improve- merical errors due to the discontinuity of the gradi- ments can be mentioned as Hu and Chen 19, Huang et ent of the basis functions 1. This is called the “cell- Copyright 20 21 22 2 al. , Nairn , Ma et al. . © VNU-HCM Press. This is an open- crossing error” . In order to alleviate the effect of access article distributed under the this phenomenon, different methods were proposed. This study using the BSMPM to mitigate the cell- terms of the Creative Commons Bardenhagen et al. proposed the Generalized In- crossing error. The BSMPM and MPM have differ- Attribution 4.0 International license. terpolation Material Point method (GIMP) 10. Vari- ences in computational grid definition. Therefore, the ants in the GIMP branch were also introduced, Stef- contact algorithm for MPM cannot be directly applied fen et al. proposed the Uniform GIMP (uGIMP) 11, to BSMPM. In this paper, the contact algorithm is Cite this article : Lo V S, Nguyen N T, Nguyen M N, Truong T T. Impact force analysis using the B-spline material point method. Sci. Tech. Dev. J. – Engineering and Technology; 4(1):721-729. 721 Science & Technology Development Journal – Engineering and Technology, 4(1):721-729 modified to a suitable form to the BSMPM. The im- plementation steps are mentioned in Section 2.3. The contact force obtained from impact of two elastic ob- jects are compared with the result from FEM, a slight difference between FEM and MPM (and BSMPM) re- sults is observed and explained in Section 3. METHODOLOGY B-spline basis functions Considering a vector containing non-decrease val- ues X=fx1;x2;:::;xn+d;xn+d+1g, where n is the num- ber of basis functions, d is the polynomial order. Each value in this vector is called knot and satis- fies the relation x1 ≤ x2 ≤ ::: ≤ xn+d ≤ xn+d+1. Vector X contains a sequence of knots is called the knot vector 2. The B-spline basis functions are con- structed by a knot vector. A uniform knot vector is a knot vector containing equally distributed knots, e.g X=f0;1;2;3;4;5g is a uniform knot vector. From the relation of the knots sequence, one notices that the value of adjacent knots can be repeated, if x1 and xn+d+1 are repeated d+1 times, it is an open knot vec- tor 2, e.g X=f0;0;0;1;2;3;4;5;5;5g is an open knot vector with n = 7 and d = 2. The i-th B-spline basis function of order d (Ni;d) is defined by using Cox-de Boor recursion formula 15. Firsly, the zeroth order (d=0) basis function must be Figure 1: (a) Quadratic B-spline basis functions (d defined X = ( = 2) built from an open uniform knot vector f0;0;0;0:5;1;1;1g and (b) Cubic B-spline basis func- 1 i f xi ≤ x ≤ x N = i+1 (1) tions (d = 3) defined by X = f0;0;0;0:5;1;1;1g i;0 0 otherwise the non-zero intervals [xi; xi+1) are called knot 2 ≥ spans .After obtaining Ni;0, higher order (d i ) ba- where p and q are the order of the univariate basis sis functions are defined as the formula below function. x − xi Two important properties of B-spline basis functions Ni;d (x) = Ni;d−1 (x) xi+d − xi are: they are non-negative for all x and the functions x − x (2) n i+d+1 x have the partition of unity property, i.e. ∑ Ni;d = + x − x Ni+1;d−1 ( ) i=1 i+d+1 i+1 1 15. in which the fraction 0/0 is assumed to be zero. Fig- ure 1 shows the high order B-spline basis functions B-spline Material Point Method (d=2, d=3). In 2D BSMPM, the computational domain is dis- 15 The derivatives of basis function Ni;d fxg are calcu- cretized by a parametric grid . This grid is de- lated as following fined by two open{ knot vectors on two} orthogo- X = x ;x ;:::;x ;x I = dN (x) {nal directions 1 2 } n+p n+p+1 and i;d d x = Ni;d−1 ( ) h1;h2;:::;hm+q;hm+q+1 as shown in Figure 2. The dx xi+d − xi d (3) numbers of basis functions in x and h direction are n − Ni+1;d−1 (x) xi+d+1 − xi+1 and m, respectively, so the total number of basis func- tions is n × m. A tensor product grid with the total In two dimensions, the bivariate B-spline functions of n × m nodes is constructed as shown in Figure 3, can be built from the tensor product of the univari- each node of this grid corresponds to one B-spline ba- ate ones 8 sis function as defined in Eq. (4). For example, the Ni; j (x;h) = Ni;p (x)Nj;q (h) (4) node with the position (1, 3) on the grid corresponds 722 Science & Technology Development Journal – Engineering and Technology, 4(1):721-729 to the basis function N1;3 (x;h) = N1;p (x)N3;q (h). The figure also shows a particle located in the upper- All the nodes on this tensor product grid are defined middle cell, so this particle is mapped to [7, 8, 9, 12, as control points in BSMPM (the same role for grid 13, 14, 17, 18, 19]. node in MPM), and in practice these control points are arbitrary distribution 15. Figure 4: A quadratic (d=2) BSMPM grid Unlike the original MPM, the particles in BSMPM are considered in the whole discretized domain, instead Figure 2: A 2D parametric grid constructed from of a specific cell, as shows in the equation below 23 two open knot vectors X and I x − x x = min ; xmax − xmin y − y (5) h = min ymax − ymin where (xmin, ymin) is the lower-left control point and (xmax, ymax ) 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 2 3 [ ] [ ] ¶x ¶x [ ] ¶N ¶N 6 ¶x ¶y 7 ¶N − = 4 5 = J 1 (6) ¶x ¶x ¶h ¶h ¶x x ¶x ¶y where J is the Jacobian matrix and defined by 2 3 ¶x ¶x 6 ¶x ¶h 7 = 6 7 Figure 3: A tensor product grid containing n × m J 4 ¶y ¶y 5 (7) nodes (control points) ¶x ¶h and the components are computed as As shown in Figure 4 a second order (quadratic) ¶x ¶N (x) BSMPM grid. The cell is made from 3 knot spans in x ¶x = ∑ PA ¶x (8) A=1 direction and 2 knot spans in y direction, so the num- ber of knots in knot vectors are 4 and 3, respectively. where P denotes the coordinates of the control points The number of control points in x direction is 3 (knot and A is the global index of control point 23. spans) + 2 (order) = 5 and in y direction is equal 4. In the BSMPM, for convenient the knot These control points play the role of grid nodes in the vector for an interval [0;L] is defined by original MPM, the knots from knot vectors are only X = f0;:::;0;4x;24x;:::;L − 4x;L;:::;Lg, used for creating a computational grid. At can be seen where 4x denotes the length of knot span 15. in Figure 4, each cell has 9 control points, for example, And note that the knot vector must be normal- the lower-left cell related to [1, 2, 3, 6, 7, 8, 11, 12, 13]. ized before a parametric grid is created, so the 723 Science & Technology Development Journal – Engineering and Technology, 4(1):721-729 knotn vector is rewritten as the followingo form where GI is the derivatives of the B-spline basis func- X0 4x 24x L−4x = 0;:::; L ; L ;:::; L ;1;:::;1 , this is also tions. a difference in the parameter space between B-spline Before applying into Eq. (9) for checking contact, the basis functions and other functions. The control normal vector in Eq. (13) must be normalized 23 points are arbitrary distributed and they are in the i i nI physical coordinate (x, y). nI =  i  (14) nI Contact algorithm The implementation of contact algorithm into the This section presents the algorithm proposed by Bar- BSMPM algorithm can be summarized as following denhagen et al. 18 and makes appropriate modifica- steps: tion to apply into the BSMPM. Step 1: Mapping data from particles to control points When two bodies are approaching each other, there 1. Compute the mass of( I-th) control point from the is a region where they have some of the same control i;t ∑ i;t i i-th body: mI = p NI xp Mp points. These control points are viewed as the con- 2. Compute the momentum( of I-th) control point from tact points, the contact algorithm is applied on these i;t ∑ i;t i the i-th body: pI = p NI xp (Mv)p points only. 18 3. Compute external force at control point I from i-th In the contact region, the following equation is used ext;i;t body: f as a condition to check if two bodies are in contact or I 4. Compute internal force at control( point) I from i-th release int;i;t −∑ i;t s i;t ∇ i;t ( body: fI = p Vp p NI xp ( ) ≥ i;t − cm i 0 contact i;t vI vI :nI = (9) 5. Compute the total force at control point I: fI = < 0 release ext;i;t int;i;t fI + fI where i denotes the i-th body in the computational Step 2: Update the control point momentums: cm 23 i;t+4t i;t i;t 4 domain, vI is the center-of-mass velocity of the pI = pI + fI t control point I-th for each pair in contact Step 3: Imposed boundary conditions at specific con- 1;t+4t 2;t+4t trol points (if needed) cm pI + pI vI = (10) Step 4: Contact force calculating (for contact points m1;t + m2;t I I only). i In Eq. (9), nI is the normal vector of control point I-th 1. Calculate the normal vector from Eq. (14) of body i-th and computed as following steps. 2. Calculate the center-of-mass velocity using Eq. Firstly, the density rc for each cell in contact state is (10) computed as below 23 3. Check the contact condition in Eq. (9) np ( ) If two body are in contact, continue sub-step 4 and 5. ri 1 i 2 i − i c = ∑ mpS xp xc (11) If not, move to Step 5. Ve p=1 4. Compute contact( force at contact) control points I: contract;i;t mi;t cm;t i;t where Ve is volume of cell e-th, xc is the center of cell I − fI = 4ti vI vI e-th. Remember that in the BSMPM each cell is made 5. Correct the control point momentums: of knot spans (see Figure 4). correct;i;t i;t+4t contract;i;t fI = fI + fI x In 1D, the function S (x) is given by the following def- Step 5: Mapping data from control points to particles inition 23 i;t+4t i;t 1. Update particle( velocities:) vp = vp + 4 ( ) Sx (x) t ∑ N xt f i;t + f contact;i;t 8 mi;t I I p I I > 1 3 9 3 1 I > x2 + x + ; − ≤ x ≤ − i;t+4t i;t > 2h2 2h 8 2h 2h 2. Update( ) particle positions: xp = xp + 1 3 1 1 4t ∑ t correct;i;t+4t − x2 + ; − ≤ x ≤ (12) i;t I NI xp pI 2 mI = > h 4 2h 2h > 1 2 3 9 1 3 3. For MUSL only, get control point velocities: > − ≤ ≤ i;t+4t correct;i;t+4t i;t > 2 x x + ; x :> 2h 2h 8 2h 2h vp = pI =mI ; i;t+4t 0 otherwise 4. Compute( ) particle gradient velocity: Lp = ∑ ∇ i;t i;t+4t The function S2 (x;y) in Eq. (11) is obtained by mul- I NI xp vI tiplying two 1D functions S2 (x;y) = Sx (x)Sy (y): 5. Update( particle gradient) deformation tensor: i;t+4t i;t+4t i;t Finally, the normal vector of control point I-th is ob- Fp = I + Lp 4t Fp tained 23 i;t+4t 6.( Update) particle volume: Vp = ( ) 4 i ∑ i ri det Fi;t+ t V i;0 nI = c GI xc c (13) p p 724 Science & Technology Development Journal – Engineering and Technology, 4(1):721-729 7. ( Compute) strain increment: 4ep = i;t+4t sym Lp 4t, then compute the stress in- crement 4sp. i;t+4t i;t 8. Update particle stresses: sp = sp + 4sp Then, reset the computational grid and move to the next time step. RESULTS Two numerical examples are presented in this section, particularly: • Collision of two circular disks. • Collision of a circular disk onto a rectangular block. Figure 5: Impact of two circular disks. The first example investigates the contact of two cir- cular surface with the same material. The second ex- ample studies the contact of a soft circular surface and hard flat surface. To validate the results from these two examples, cor- responding FEM models are created from ABAQUS software. FEM model is prepared with very fine mesh and set up with the same parameters and initial con- ditions 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:2 m and the thickness is one unit. The material properties used in this prob- Figure 6: Kinetic and strain energy. lem are: Young’s modulus E = 1000 Pa, Poisson ra- 3 tio v = 0:3, and the mass density r0 = 1000 kg=m . The coordinate of the center of the lower-left disk is loss from friction, the strain energy loss in BSMPM is (0:2; 0:2), the upper-right disk is (0:7; 0:7), two disks caused by other error factors. × 2 are in a square domain of size 0:9 0:9 m . The initial The variation of contact force during the impact pro- velocities of the particles v = (0:1; 0:1) m=s, for the cess is shown in Figure 7. The FEM model used to upper-right disk, the velocities of the particles are set simulate this problem has 3288 nodes. The results − to vp = v and for the lower-left vp = v. from MPM and BSMPM show that the impact of two The computational domain is discretized into 40×40 bodies occurs earlier than the result in FEM as men- knot spans. Each computational cell has 9 particles. tioned before. This is because the contact force in The original MPM with Lagrange basis and quadratic MPM is computed in the node of the computational BSMPM (d=2) are concerned in this example. grid (or control point in BSMPM), not in the particle The time step for this simulation is chosen as 4t = of the body, so when two bodies approach the contact 0:001 s, the total simulation time is 3 s. So, there is region and have the same control points, the contact 3000 steps in this simulation. is detected immediately although two bodies have not The kinetic and strain energy obtained from BSMPM touched each other yet. In FEM, the contact is only and FEM is shown in Figure 6. Kinetic energy in detected when two bodies touch each other, so the BSMPM decreases earlier than the result from FEM contact force obtained in FEM is later than MPM. and strain energy in BSMPM increases earlier. This is The contact force obtained from BSMPM using higher reasonable for the contact in BSMPM algorithm and order B-spline functions also shows the smooth curve will be explained in the comment of Figure 7. The compared to the MPM and FEM. value of kinetic energy in both case are the same, while Figure 8 shows the von-Mises stress field during the the strain energy in BSMPM is lower than FEM. Both impact process of two disks using the BSMPM. In de- case are in frictionless contact, so there is no energy tail, two disks approaching each other in Figure 8 (a), 725 Science & Technology Development Journal – Engineering and Technology, 4(1):721-729 Figure 7: Impact force obtained from FEM, MPM and BSMPM (d=2) Figure 9: A circular disk collides with a rectangular block. 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). Af- ter impact, two disks move far away as shown inFig- ure 8 (e). Collision of a circular disk onto a rectangu- lar block In this example, a circular disk collides onto a stiffer rectangular block, as shown in Figure 9. The radius of the circular disk is R = 0:2 m, and the thickness is one unit. The material properties used for circular disk are: Young’s modulus E1 = 1000 Pa, Figure 10: Impact force of example 5.2 obtained Poisson ratio v1 = 0:3, and the mass density r1 = from FEM, MPM and BSMPM (d=2) 1000 kg=m3. The rectangular block is made from 6 stiffer material with Young’s modulus E2 = 10 Pa, r Poisson ratio v2 = 0:3, and the mass density 1 = FEM and MPM. 3 × 2 5000 kg=m , the rectangular size is 1 0:2m . Dis- Figure 11 shows the collision of two objects, the von- tance between the center of the circular disk to the top Mises stress field and maximum stress field are pre- of rectangular block is 0:3 m. The computational do- sented. × 2 main is a square with dimension of 1:2 1:2 m . The To investigate the convergence of BSMPM and MPM, − initial velocity of the disk is v = (0; 0:2) m=s. In this the computational domain with a set of 60 × 60 knot simulation, the gravitational acceleration is ignored. spans is retained. Different numbers of particles per The computational domain is discretized into a set of cell (PPC) 4, 9 and 16 are analyzed. Figure 12 shows × 60 60 knot spans. Each cell has 9 particles. The the total energy of the system respect to time. From nodes (or control points) on the bottom line of the the initial conditions, the total energy can be com- rectangular is fixed in two direction x and y. puted as rpR2tv2=2 = 2:512 J and plotted by the The time step size is chosen as 4t = 0:001 s, and the black line in the figure. As shown in Figure 12, the total simulation time is 2 s. So, there is 2000 steps in case of MPM with PPC = 4 gives a very large devi- this simulation. ation, and when PPC = 9, the result is significantly The contact force obtained in this example also shows improved. In the case of BSMPM, there are no signif- the similarity to the conclusions from the previous ex- icant deviations and the results are slightly improved ample. Figure 10 also shows that the impact occurs when increasing PPC. earlier in BSMPM, because BSMPM has more control points (nodes) than MPM so the contact is detected DISCUSSIONS earlier. Similarly to the previous example, the con- As present in Section 3, there is a slight difference tact force in BSMPM is smoother than the curve from in the results of MPM, BSMPM and FEM. The mag- 726 Science & Technology Development Journal – Engineering and Technology, 4(1):721-729 Figure 8: Impact of two circular disks. Figure 12: Energy of the system respect to time with different number of particles per cell. 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 Figure 11: Circular disk deforms during the impact determined automatically when two objects have the process onto a stiffer surface. (a) von-Mises stress same grid node. In FEM, the contact is only detected and (b) Maximum stress. 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. nitude of the contact force in the three methods is the same, but the impact occurs earlier in MPM and CONCLUSIONS BSMPM. This is because the contact force in MPM is The contact algorithm has been successfully modi- computed in the node of the computational grid (or fied and applied into the BSMPM. The contact force control point in BSMPM), not in the discrete parti- obtained from this research is compared to FEM. A cle of the object, so when two objects approach the slight difference in the result is observed, this is be- contact region and have same the grid node (control cause the contact force is calculated at the control points), the contact is detected immediately although point in the computational grid instead of the discrete 727 Science & Technology Development Journal – Engineering and Technology, 4(1):721-729 particles. This is an inherent property of the MPM al- 8. Nairn JA. Material point method calculations with ex- gorithm, so it is inevitable. More study on the contact plicit cracks. Computer Modeling in Engineering & Sciences. 2003;4:649–663. algorithm need to be done to overcome this disadvan- 9. II ARY, Sulsky D, Schreyer HL. Fluid-membrane interac- tage and improve the accuracy of the method. tion based on the material point method. International Jour- nal for Numerical Methods in Engineering, vol. 48, pp. 901- ACKNOWLEDGMENT 924, 2000.;Available from: https://doi.org/10.1002/(SICI)1097- 0207(20000630)48:63.0.CO;2-T. This research is funded by Ho Chi Minh City Univer- 10. Badenhagen SG, Kober EM. The generalized interpolation ma- sity of Technology – VNU-HCM, under grant num- terial point method. Computer Modeling in Engineering & Sci- ence. 2004;5(6):477–495. ber T-KHUD-2020-47. We acknowledge the support 11. Steffen M, et al. Examination and analysis of implementation of time and facilities from Ho Chi Minh City Uni- choices within the material point method (MPM). Computer versity of Technology (HCMUT), VNU-HCM for this Modeling in Engineering & Science. 2008;31(2):107–127. 12. Sadeghirad A, Brannon RM, Burghardt J. A convected par- study. ticle domain interpolation technique to extend applicability of the material point method for problems involving mas- ABBREVIATIONS sive deformations. International Journal for Numerical Meth- ods in Engineering. 2011;86(12):1435–1456. Available from: BSMPM: B-spline Material Point Method https://doi.org/10.1002/nme.3110. FEM: Finite Element Method 13. Zhang DZ, Ma X, Giguere P. Material Point Method enhanced MPM: Material Point Method by modified gradienet of shape function. Journal of Computa- tional Physics. 2011;230(16):6379–6398. Available from: https: MUSL: Modified Update Stress Last //doi.org/10.1016/j.jcp.2011.04.032. 14. Steffen M. Examination and Analysis of Implementation CONFLICT OF INTEREST Choices within the Material Point Method (MPM). Computer Modellin in Engineering and Sciences. 2008;31(2):107–127. Group of authors declare that this manuscript is origi- 15. Gan Y, et al. Enhancement of the material point method using nal, has not been published before and there is no con- B-spline basis functions. Numerical Methods in Engineering. flict of interest in publishing the paper. 2017;113(3):411–431. Available from: https://doi.org/10.1002/ nme.5620. 16. Wobbes E, Moller M, Galavi V, Vuik C. Conservative Tay- AUTHOR CONTRIBUTION lor Least Squares reconstruction with application to material Vay Siu Lo is work as the chief developer of the point methods. International Journal for Numerical Methods in Engineering. 2018;117(3):271–290. Available from: https: method and the manuscript editor. //doi.org/10.1002/nme.5956. Nha Thanh Nguyen and Minh Ngoc Nguyen take part 17. II ARY, Sulsky D, Schreyer HL. The material point method in the work of gathering data and checking the numer- for simulation of thin membranes. International Journal ical results. for Numerical Methods in Engineering. 1999;44(10):1429– 1456. Available from: https://doi.org/10.1002/(SICI)1097- Thien Tich Truong is the supervisor of the group, he 0207(19990410)44:103.0.CO;2-4. also contributes ideas for the proposed method. 18. Bardenhagen SG, et al. An improved contact algorithm for the material point method and application to stress propagation REFERENCES in granular material. Computer Modeling in Engineering and Sciences. 2001;2(4):509–522. 1. Zhang X, Chen Z, Liu Y. The Material Point Method 19. Hu W, Chen Z. A multi-mesh MPM for simulating the A Continuum-Based Particle Method for Extreme Loading meshing process of spur gears. Computers and Structures. Cases. Elsevier. 2017;Available from: https://doi.org/10.1016/ 2003;81(20):1991–2002. Available from: https://doi.org/10. B978-0-12-407716-4.00003-X. 1016/S0045-7949(03)00260-8. 2. Tielen R. High order material point method. master disserta- 20. Huang P, Zhang X, Ma S, Huang X. Contact algorithms for the tion, Delft University of Technology, Delft. 2016;. material point method in impact and penetration simulation. 3. Jassim I, Hamad F, Vermeer P. Dynamic material point method International journal for numerical methods in engineering. with applications in geomechanics. Stuttgart University, Ger- 2010;85(4):498–517. Available from: https://doi.org/10.1002/ many. 2011;. nme.2981. 4. Nguyen VP. Material point method: basics and applications. 21. Nairn JA. Modeling imperfect interfaces in the material point Cardiff University, Department of Civil Engineering. 2014;. method using multimaterial methods. Computer Modeling in 5. Wang YX, et al. Response of multi-layered structure due to im- Engineering and Science. 2013;92(3):271–299. Available from: pact load using material point method. Engineering Mechan- https://doi.org/10.32604/cmes.2013.092.271. ics. 2007;24(12):186–192. 22. Ma J, Wang D, Randolph MF. A new contact algorithm in the 6. Hu WQ, Chen Z. Model-based simulation of the synergis- material point method for geotechnical simulations. Interna- tic effects of blast and fragmentation on a concrete wall us- tional Journal for Numerical and Analytical Methods in Ge- ing the MPM,” International Journal of Impact Engineering. omechanics. 2014;38(11):1197–1210. Available from: https: 2006;32(12):2066–2096. Available from: https://doi.org/10. //doi.org/10.1002/nag.2266. 1016/j.ijimpeng.2005.05.004. 23. Nguyen VP. Material point method: basis and application. 7. Andersen S, Andersen L. Analysis of spatial interpolation 2014;Available from: https://researchgate.net/publication/ in the material-point method. Computers and Structures. 262415477_Material_point_method_basics_and_applications. 2010;p. 506–518. Available from: https://doi.org/10.1016/j. compstruc.2010.01.004. 728 Tạp chí Phát triển Khoa học và Công nghệ – Kĩ thuật và Công nghệ, 4(1):721-729 Open Access Full Text Article Bài Nghiên cứu Phân tích lực va chạm bằng phương pháp Điểm vật liệu sử dụng hàm dạng B-spline Lồ Sìu Vẫy1,2,*, Nguyễn Thanh Nhã1,2, Nguyễn Ngọc Minh1,2, Trương Tích Thiện1,2,* TÓM TẮT Trong giải thuật MPM, các điểm vật liệu được xây dựng trong một trường vận tốc đơn trị nên sự tương tác/tiếp xúc không trượt giữa