=x*x*x/{\\x\\? (19) where x(r) is the piecewise-constant basis function: X(x) \\x\\ < -jl otherwise (20) and I is width of %. The B-spline basis function associated with node i, is then 0,-(jc) = lp and gradient weighting function V0,- for a piecewise-constant xP with a characteristic length of I. Notice that that while j looks smooth in Fig. 4(a), the dashed lines show locations of breaks in continuity which become apparent in Vtpjp in Fig. 4(b). These breaks in continuity will become important in later analysis. (a) GIMP Weighting Function (b) GIMP Gradient Weighting Function Figure 4: Example GIMP weighting function (f>jp, and gradient weighting function V0,- centered at 0 using piecewise linear grid basis functions and piecewise constant particle characteristic functions Xp- Dotted lines denote breaks in the continuity of the functions. 3.4 Kinematic Boundary Conditions One of the conveniences afforded by the use of a Cartesian background grid is the ease of applica Implementation Choices within the Material Point Method 115 tion of kinematic boundary conditions. That is, Dirichlet or Neumann conditions, or a combination, on the velocity field. Note that if no treatment is given to the boundary nodes, then particles are able to freely advect from the computational domain in what could be considered a zero gradient Neumann condition. This important part of the algorithm has received scant treatment in the literature (although it is very relevant when one is actually implementing MPM), so we turn attention to it here. 3.4.1 Traditional MPM In traditional MPM, boundary conditions only need be applied on those nodes which coincide with the extents of the computational domain. As illustrated in Fig. 2(a) nodes beyond those boundaries are not influenced by particles within the domain. This can be considered a result of the zero width of the Dirac delta characteristic functions. Boundary conditions must be applied to the velocity that has been projected to the nodes (Eq. 2), the time advanced velocity (Eq. 6), and the acceleration (Eq. 5). For Dirichlet conditions, this simply means overwriting the calculated values for the velocities with the prescribed values. For the acceleration, some debate exists regarding the proper means of treatment. The usual approach has been to assume that a Dirichlet condition for velocity implies that the acceleration should be zero on those boundary nodes. However, it is also possible that if Eq. 6 were solved for acceleration: a,- = (vf+1-vf)/At, (30) the proper value for the acceleration at the boundary nodes would be computed based on the difference between the time advanced velocity (after application of boundary conditions) and the projected velocity (without the application of boundary conditions). Put another way, acceleration on the boundary nodes can be considered to reflect the force required to bring the velocity at those nodes from the projected value to the prescribed value. While these two approaches to the acceleration seem substantially different, the difference in simulation results is very subtle. Indeed, when both approaches are tested with the manufactured solution described in Sec. 5.1, the superiority of either is not apparent. Currently, the UCF implementation uses the boundary treatment given in Eq. 30. In addition to prescribed velocity boundary conditions, \"symmetry\" boundary conditions are also frequently useful. Symmetry BCs are used to represent a plane of symmetry, which allows the use of a reduced computational domain, or a friction- less surface. They are achieved by simply applying a zero velocity Dirichlet condition on the component of velocity normal to a boundary, while allowing the other components to remain at their computed values. Acceleration is handled in the same manner, with the normal component either zeroed out, or computed as in Eq. 30. 3.4.2 GIMP and B-Spline MPM When using GIMP or B-Spline MPM, there are additional considerations in the applications of the boundary conditions. Namely, because of their increased extents, it is possible for particles to influence, and be influenced by, nodes that lie outside of the simulation domain, (see Figures 2(b) and 2(c)). In the UCF, these are referred to as \"extra\" nodes, but may also be called \"ghost\" nodes by other investigators. Boundary condition treatment of these nodes for Dirichlet conditions is the same as for the regular boundary nodes, namely, their computed values are replaced by prescribed values as described above. In treating symmetry boundaries, the extra nodes require special care. In particular, the normal component of velocity for these nodes is no longer set to zero, but rather should be set to the negative of the value of the node opposite the boundary. The need to do so is apparent if one considers two objects approaching a collision plane symmetrically. The normal component of velocity on the opposite sides of that plane will have opposing signs. 4 Analysis and Interpretation MPM is a fairly new method and thus there has been a recent push by the MPM community to provide a more formal analysis of errors in the116 Copyright © 2008 Tech Science Press CMES, vol.31, no.2, pp.107-127, 2008 method. Bardenhagen (2002) looked at energy conservation errors in MPM, focusing on the effects of the choice of two time-stepping algorithms. Recently, Wallstedt and Guilkey (2008) expanded on the analysis of those time-stepping algorithms. Love and Sulsky (2006a) and Love and Sulsky (2006b) analyze an energy consistent implementation of MPM, the second of these showing an implicit implementation to be unconditionally stable and energy-momentum consistent. Wallstedt and Guilkey (2007) focus on velocity projection errors and present a scheme which helps ameliorate these errors. Steffen, Kirby, and Berzins (2008) perform an analysis on some of the spatial integration errors present within MPM. In this section, we continue adding a few more pieces to the error analysis of MPM. Specifically we will look at integration errors which are affected by the smoothing of the piecewise-linear basis functions. 4.1 The Relationship between GIMP and B- Splines Taking a closer look at the weighting function (Eq. 26), we see that the construction is essentially a convolution of the grid basis functions and the particle basis function %p. Since a standard piecewise-linear tf>j can also be represented as the convolution of piecewise-constant basis functions, we can rewrite the GIMP weighting function as: 4>=Xg*Xg*Xp/(\\Xg\\\\Xp\\) (31> where the width of Xg is h (the grid spacing), and the width of %p is lp, as described in the GIMP methods. The equivalent GIMP basis function would then come from evaluating Eq. 31: Q.(x) = 4>(x-Xi) (32) with the GIMP weighting function equivalent to evaluating Eq. 32 at the particle position, xp. The reason for rewriting the GIMP basis functions in this way is to demonstrate the similarities between the construction of GIMP basis functions and the construction of B-spline basis functions as in Eq. 19. Both basis functions are constructed by convolving piecewise-constant basis functions with themselves; however all of the % in the B- spline basis are of width h while one of the % functions used in the GIMP method has width lp. In cpGIMP, the particle characteristic length lp (of which there may be different lengths for different directions) is updated in time, meaning the cpGIMP weighting function (Eq. 29) is time dependent, and is different for each particle p. The ideal case would be that the updating of lp in time will cause the set of particle characteristic functions xp to perfectly tile space, but due to the rectilinear constraints of %p, this is not possible in general multi-D simulations. Because of this inability to tile space, and the recognition that the major benefit of GIMP is the smoother equivalent basis functions, a simplified standard GIMP is used in which lp = P for all time. Furthermore, if lp = I is constant for all particles p in a simulation (uGIMP), the effect is truly equivalent to using standard MPM with a smoother set of basis functions. In fact, if one were to disassociate the smoothing length, /, from the particles in a uGIMP formulation and instead leave / as a free parameter, the effect is to create quadratic B-spline-like basis functions, with / determining the maximum extent of the functions. Choosing / = 1^ would give standard GIMP. Choosing I = h would give quadratic B-spline basis functions. Choosing I = 0 would lead to the degenerate case of Xp = S(x - xp), leaving us with the standard piecewise-linear basis functions. It has been our decision to leave the smoothing length I as a free parameter the UCF, allowing for various options when running simulations. We will explore various choices of I in the sections to follow. 4.2 Smoothing Length Dependent Integration Errors Spatial integrals within MPM are performed using nodal integration - an approximation which takes the form: f f{x)da^f{xp)Vp. (33) JQ pImplementation Choices within the Material Point Method 117 An analysis of errors in the above approximation within the MPM framework was performed by Steffen, Kirby, and Berzins (2008). There, the nodal integration approximations in MPM were equated to non-uniformly-spaced midpoint integration of functions with discontinuities in various derivatives. In particular, the analysis focused on the errors when calculating the internal force (Eq. 4), which involves the following approximation: pnt __ Q (*) • VipVp (34) The main result from that analysis showed that if the particle arrangements did not respect the discontinuities which arise from the basis functions (i.e. a particle's voxel overhangs node boundaries), an extra integration, or \"jump\" error can arise in the above approximation which is of the order C{{f'p+^]]Axp+2, where the function, /. being integrated is Cp continuous (with p = - 1 for discontinuous functions). Here, [[•]] represents the jump in the p + 1 derivative of / at the discontinuity and Ax is the particle spacing. Note that the function V^( in Eq. 34 is discontinuous, thus a jump error of ff(Ax) can arise in MPM when using standard piecewise-linear basis functions, depending on particle spacing. Numerical examples of this error were shown in Steffen, Kirby, and Berzins (2008). The jump error for a single particle is calculated as Ejump = fQp f(x) dx - f(xp)Ax, where Qp spans a discontinuity, or jump, and Ax is the width of the particle. This consists of measuring the midpoint integration error for the single interval spanning the jumps. Integration approximations, such as in Eq. 33, involve integration over the whole domain, using multiple intervals, leading to a composite midpoint rule integration error which is ff(Ax2). These two errors are additive, giving a total error of the form Etotul = j fix) dQ £ f(xp) Vp = Emp + EJump, (35) where Emp is the composite midpoint error and Ejump is any errors arising from integrating across jumps. Note that if we assume particles are nonoverlapping and fill space, this equation can also be written as Etotai-'y. P Q, f(x)dn-f(xp)Vp (36) Since the errors are additive and since Emp is always ff(Ax2), Co and higher continuous functions exhibit an overall integration error which is ff(Ax2), while discontinuous functions have an error which is ff(Ax). Again, this is important because the nodal integration for the internal force calculation in Eq. 4 involves the gradients of the basis functions, which are discontinuous at grid cell boundaries when the standard piecewise- linear basis functions are used. In uGIMP, we have the choice of a smoothing parameter I (the width of our general particle characteristic function %), independent of the individual particle sizes, which ensures us C\\ continuous basis functions but leads to a situation which was not analyzed in Steffen, Kirby, and Berzins (2008) - the case where the width of the particle is greater than the smoothing length. In such cases (as illustrated in Fig. 5), a particle can span two jumps in the continuity of the basis functions. Consider a general case from Eq. 36 where a single particle, or Qp, spans three regions of a piecewise linear function. The first region (R\\) is defined by the equation y\\ = a\\x + b\\, the second (/?2) will be >'2 = c/2X ■ l>2- and the third (R3) is >'3 = a^x + Z>3 with the particle located a distance 8 inside the second region (see Fig. 5). For this118 Copyright © 2008 Tech Science Press CMES, vol.31, no.2, pp.101-121, 2008 a^=0 [ + a2 = -2/(hl) ] 83 - 0 + Ax/2-5 5 1-5 Ax/2-1+6 (b) GIMP Specific Two-Jump Situation Figure 5: Example cases of a particle's volume (the area between the square brackets) spanning two jumps in a piecewise-linear function: (a) shows a particle spanning two jumps in a general piecewise-linear function with a\\ and i = 1.2.3 the parameters describing the linear segments, while (b) shows the specific piecewise- linear uGIMP gradient function V0 and a situation where the particle size Av is greater than the smoothing length /. case, the integration error is given by: EjU,„P= I f(x)dx-f(xp)Vp (37) Jiip = / \\'\\dx+ / \\'->dx+ JQrnRt' ./Q,,n/?2' / >'3 dx - }'2 (xp) A v (38) J Q,,r 1R} 1 7 1 = - («3 - Cl2 ) + - («2 - «3 ) / Av+ 1 , («2 - «3)/<5 + \"(«3 - (11)8\" + 1 1 7 -(«! - 2«2+«3)5Av+ -(«3 - «[)Av\" 2 8 (39) where I is the width of the center region and 8 is the particle offset into the center region. Since I and 8 are both less than Av, this whole expression appears to be &'(A.\\2). However, when we consider the specific case of measuring the integration error in internal force (Eq. 34 with a = 1) when a particle spans the center region of V0 as in uGIMP (see Fig. 4(b) and Fig. 5(b)), the slopes of the left and right regions in Fig. 5(b) are zero, while the slope of the center region is dependent on the smoothing length I and the grid spacing h. Specifically, for these regions, ci\\ = 0, r/3 = 0 and ci2 = -2/(hi). When we substitute these parameters into Eq. 39, we are left with -1 p r= ___ Gjump j Av - / + 2 | (40) Here, it is clear that the jump error in this case is £>'(Av). If instead, Av < / and the particle only spans one of the uGIMP jumps, the error then takes the form: jump = -77- [<52 - <5Av + 7 Av2]. hi 4 (41) Since / no longer depends on Av, this is now f/(A\\2). To test this analysis, we calculate the force on a single node for a set of particles with constant stress a. This is the same test performed by Steffen, Kirby, and Berzins (2008) when looking at a particle spanning a single jump instead of the two-jumps analyzed above. In this case, the internal force on a node is calculated as fi = -I V0, P ip ■gpVp = -o'2V*,(xp)Vp. (42) For a constant stress, internal force should be zero, so any errors are from integrating V0,-. Fig. 6 shows the errors for various particle spac- ings when the smoothing length / = 1/10. As expected, when Av < / the error converges as I>'(A\\2). When Av becomes greater than /, the error tends towards £>'(Av). Here, we have shown errors in the internal force which are either £>'(Av) or I>'(A\\2), depending on the relationship between the smoothing length / and the particle widths Av. For stability, in addition to the typical CFL constraints one needs when smooth forces exist, we need to consider further time step restrictions when force kicks arise from these integration errors. These time step restrictions would be similar to those required, as shown by Tran, Kim, and BerzinsImplementation Choices within the Material Point Method 119 1/PPC Figure 6: Errors in internal force vs. particle spacing for constant <7=1, grid spacing h = 1, and smoothing length / = 1/10. The particles have a global uniform density, however they have a locally non-uniform spacing. Otherwise, super- convergent results are observed. Note that when particle spacing is less than the smoothing length /, the error converges as second-order. As the particle spacing becomes greater than the smoothing length, the error tends towards first-order. (2007), due to force kicks arising from grid crossing errors. While a time step of Ati may be sufficient when we are in the 6 (Ax1) error region, a smaller Atj may be required to control stability when we are in the 6(Ax) error region. 4.3 Impact of Boundary Treatments In MPM, the union of the particles' voxels are assumed to fill space and define the material of interest. However, many calculations are not performed directly on the particles, but rather on the background grid to which the particle information is projected. This projection of particle information leads to a set of \"active\" basis functions and grid cells (those which have particles in their support) which, in general, will differ geometrically than the union of particles' voxels. This can, and does, lead to a further errors in many MPM simulations. In standard MPM with piecewise-linear basis functions, the active grid cells are those which contain a particle. One could argue that a grid cell which contains no particles but still overlaps with a particle voxel (from a particle in a neighboring cell) should also be active, but is not considered so in the current MPM framework. In either case, when considering active cells on the grid, there may be a geometric error of up to h in each direction. When moving to uGIMP, or B-splines, this geometric error can become worse since the support of these basis functions are larger. Cubic B- splines, quadratic B-splines, and uGIMP can experience geometric errors of up to 2h, 3h/2 and h I 1/2. respectively. All of these errors are 6(h); however it is important to note that these geometric errors are not only a function of how well the object of interest is aligned with grid cells, but they are also a function of basis function choice. Some work has been performed on MPM background grids which more closely represent the material of interest. For example, Wallstedt (2008) has worked on an MLS representation of a material boundary and incorporates this boundary into the MPM integration routines. Here, we sidestep part of the issue by developing test problems in Section 5 whose boundaries are perfectly aligned with the grid boundaries (such as a fixed- fixed elastic bar). Even with these aligned test problems, geometric errors can still exist since information is projected to extra nodes outside the domain, as is shown in Fig. 2(b) and Fig. 2(c); information which is still used in standard kinematic boundary treatments. To illustrate this geometric error, Fig. 7 shows an example of the density field resulting from projecting particles with constant mass (a discretization of a constant density field) to the grid. The density field is calculated as p(x) = Iip4i(x) with pi = mi/(Jjp^ipVp). In this example, the constant density field spans the region [0.1] which is embedded in a grid covering the region [-0.2.1.2]. Since the deformed configuration of the material with respect to the grid is effectively the support of the fields of interest, we can see from Fig. 7 how implementing boundary conditions and modeling contact can present a challenge when wider basis functions are used. We postulate that, in general, all of the methods here can suffer from G(li) geometric errors. In the special case of boundary aligned problems,120 Copyright © 2008 Tech Science Press CMES, vol.31, no.2, pp.107-127, 2008 1.5 1 ™ n c c 0.5 Q 0 -0.5 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x Figure 7: Density fields resulting from projecting particle mass to the background grid. The true density field is shown, along with density fields calculated with piecewise-linear and quadratic B-spline basis functions. Here we can see that the geometric extent of information projected to the grid not only depends on which grid cells contain material, but also on basis function choice. -■□OOOOOOOOOOOOOOOOOOOOOOOOOOOOOir- ------True Density ------Piecewise-Linear Basis - Quadratic B-spline Basis o Particle Locations methods where information is projected to extra nodes, such as uGIMP and standard B-splines, will still be affected by this &(h) geometric error when Neumann or Dirichlet boundary conditions are applied. These methods should not be affected by this error when periodic boundary conditions are used. Methods which do not require the use of extra nodes, such as standard MPM with piecewise-linear basis functions (Fig. 2(a)), modified boundary B-splines (Fig. 3), and cpGIMP will not be affected by this geometric error. It is worth noting that while cpGIMP is implemented in the UCF using extra boundary nodes, information is not projected to these extra nodes in well-behaved boundary aligned simulations. This is because the particle p that is closest to the boundary has width lp, and is located at a position of lp/2 to the inside of the boundary, and the closest extra node is at a distance of h + lpf 2, which is the exact location where the extra node's basis function goes to zero (see Fig. 4). 5 Test Problem Development Code verification has gained renewed importance in recent decades as costly projects rely more heavily on computer simulations. Full time- dependent test problems with analytical solutions are desired so that simulation errors can be assessed. The Method of Manufactured solutions [Schwer (2002); Knupp and Salari (2003); Banerjee (2006)] begins with an assumed solution to the model equations and analytically determines the external force required to achieve that solution. This allows the user to verify the accuracy of numerical implementations, understand the effects of parameter choices within the code, and to find where bugs may exist or improvements can be made. The critical advantage afforded by MMS is the ability to test codes with boundaries or nonlinearities for which exact solutions will never be known. It is argued [Knupp and Salari (2003)] that MMS is sufficient to verify a code, not merely necessary. Since full transient mechanics solutions are often difficult to find in the literature, we will first present an overview of the method of manufactured solutions with which we will then develop both 1-D and 3-D test problems. 5.1 Method of Manufactured Solutions Overview For this paper we define several non-linear dynamic manufactured solutions and use them for subsequent testing. The solutions exercise the mathematical and numerical capabilities of the code and provide reliable test problems for ascertaining a simulation's accuracy and stability properties. Finite Element Method (FEM) texts often present Total Lagrange and Updated Lagrange forms of the equations of motion. The Total Lagrange formImplementation Choices within the Material Point Method 121 is written in terms of the reference configuration of the material whereas the Updated Lagrange form is written in terms of the current configuration. Either form can be used successfully in a FEM algorithm, and solutions from Updated and Total Lagrange formulations are equivalent [Be- lytschko, Liu, and Moran (2000)]. However, within GIMP it is necessary to manufacture solutions in the Total Lagrange formulation so that zero normal stress can be applied to free surfaces as a boundary condition. This might at first appear to conflict with the fact that GIMP is always implemented in the Updated Lagrange form. The equivalence of the two forms and the ability to map back and forth between them allows a manufactured solution in the Total Lagrange form to be validly compared to a numerical solution in the Updated Lagrange form. The equations of motion are presented in Total and Updated Lagrangian forms, respectively: VP + pob = poa (43) V