MPM Methods and Simulation Timing

These command select the MPM method to use and control time step and total time for the simulation.

Theory: Shape Functions

Many tasks in MPM involve extrapolations from particles to the grid or from the grid to the particles. These extrapolations are controlled by "Shape Functions," and the various MPM methods correspond to various methods for evaluating the shape functions. In the most generalized description of MPM (called GIMP for General Interpolation Material Point), the shape function (S_ip) and shape function gradient (G_ip) for the node i/particle p pair are:

      [math]\displaystyle{ S_{ip} = {\int_{\Omega_p} \chi_p(\vec x)N_i(\vec x) dV\over \int_{\Omega_p} \chi_p(\vec x) dV}\qquad{\rm and}\qquad G_{ip} = {\int_{\Omega_p} \chi_p(\vec x)\nabla N_i(\vec x) dV\over \int_{\Omega_p} \chi_p(\vec x) dV} }[/math]

where Ω_p is the domain for particle p, χ_p(x) is the particle basis shape function for particle p, and N_i(x) is grid shape function for node i. The various MPM methods depend on the technique used for evaluating the shape function integrals. The chart on the right shows an MPM "genealogy" or shows how the various MPM methods descend from GIMP. The methods colored in blue are the ones supported in OSParticulas and in NairnMPM.

File:MPMTree.png

The first decision is to chose χ_p(x); the chart shows three options. The choice of the Dirac delta function (δ(x)) leads to classic MPM methods. The choice of a Heaviside function (i.e., χ_p(x) = 1 in Ω_p and 0 elsewhere) leads to a family of Heaviside GIMP methods. Finally the choice of some other function (e.g., a Gaussian centered on the particle) could lead to other MPM methods, but is currently an unexplored area.

Delta Function GIMP Methods - Classic MPM

If χ_p(x) = δ(x), the shape functions reduce to the the grid shape function (e.g., S_ip = N_i(x_p)), which recovers "Classic MPM" or the original derivation of MPM. Note that "Classic GIMP predates GIMP even though it descends from GIMP in the MPM "genealogy." Although Classic MPM is the most efficient method, it is rarely used because it is not very robust for large displacement problems. Whenever particle cross cell boundaries, they can cause noise that soon leads to unsatisfactory results. These cell-crossing artifacts can be eliminated (or greatly reduced) by using different χ_p(x).

Heaviside GIMP Methods

MPM results are greatly improved by choosing χ_p(x) to be a Heaviside function. With this function, the shape function integrals reduce to integrals over the particle domain of the grid-based shape function, which are normalized to particle volume (see chart). All current MPM methods descend from this selection. The specific descendant methods depend on how they integrate over the particle domain.

uGIMP

The simplest approach is to assume the particle domain has the same shape as the initial domain and merely translates with particle motion. This method, which is called uGIMP (for undeformed GIMP), has the greatest efficiency because shape functions can be explicitly calculated for all possible particle locations relative to a regular grid. uGIMP works well from small to fairly large deformation and greatly reduces cell-crossing artifacts. Very large tension strains (>50%), however, can cause undeformed particle domains to become numerically separated, which can lead to non-physical fracture of the object. When this occurs, one option is to use a different method that accounts for deformed particle domain.

Finite GIMP

In finite GIMP, the integral is over the deformed particle domain. Although this has been shown to give improved results in 1D (where it is practical), it is too inefficient for arbitarary deformation of 2D and 3D particle domains. As a result, exactly-integrated, finite GIMP is not used.

CPDI

A useful improvement for large to massive deformation is CPDI or the Convected Partical Domain Integration approach. In CPDI, the integral over the deformed particle domain is made practical by approximating the integrand. In brief:

1. The deformed particle domain is discretized by isoparametric finite elements. In practice, the discretization uses a single element, but it could use more.
2. The Integrands for S_ip and G_ip are then approximated by expanding them in the finite element shape functions:
         [math]\displaystyle{ N_i(\vec x) \approx \sum_j N_i(c_j) M_j(\vec x)\qquad{\rm and}\qquad \nabla N_i(\vec x) \approx \sum_j N_i(c_j) \nabla M_j(\vec x) }[/math]
   where c_j are nodes of the deformed particle domain (e.g., usually corners of the domain) and M_j(x) are isoparametric finite element shape functions for the elements used to discretize the deformed particle domain.
3. Integrals over the deformed domain can then be evaluated to give S_ip and G_ip in terms for N_i(c_j).

Two versions of CPDI have appeared. The first version finds the corners of the deformed domain using the deformation gradient on the particle. If r_j0 is a vector from the particle to any edge location j in the undeformed domain, then r_j = F_pr₀, where F_p is the deformation gradient on the particle, is the vector to that location in the deformed domain. When using this approach, an initial 2D rectangle domain deforms to a parallelogram and an initial 3D orthogonal box deforms to a parallelepiped. The subclasses of this method depend on the specific finite element types used to map that particle domain:

• lCPDI for "linear" CPDI: this method maps the corners of the domain to a single linear element (4 node element in 2D and 8 node element in 3D) resulting in an isoparametric element with linear shape functions.
• qCPDI for "quadratic" CPDI: this method additionally maps midpoints of edges to nodes (e.g., 8 node element in 2D) resulting in an isoparametric element with quadratic shape functions.
• other...: this option leaves open future development such as different shape functions or discretization of each domain into more than one element.

In CPDI2 the corner locations are separately tracked during the analysis and their current locations are used to define the deformed particle domain. For example, a 2D domain deforms to an arbitrary quadrilateral element and a 3D domain to an arbitrary 8-cornered brick element.

CPDI methods are less efficient the uGIMP because of the extra overhead required to locate corner nodes needed for evaluating shape functions. CPDI2 add the additional overheard of tracking and updating corner particles. Nevertheless, the use of CPDI is crucial for handling problems with deformation beyond the capabilities of uGIMP.

cpGIMP

In cpGIMP, the deformed domain is stretched in orthogonal directions but is never sheared. In other words, particle domains that start as squares (in 2D) or cubes (in 3D) will deform to rectangles (in 2D) or orthogonal boxes (in 3D). Compared to uGIMP, cpGIMP will handle uniaxial deformation better, but it does not account for arbitrary deformation that includes shear deformation.

Scripted Commands

In scripted input files, the MPM method, the strain updating method, and the number of points per element are set with the following commands:

MPMMethod (strain update),(MPM method)
PtsPerElement (number)
TimeStep (time step)
MaximumTime (time)

where

• (strain update) - the method used to update stresses and strains on the particles. The options are:
  • USF (or 0) - update first
  • USAVG (or 2) - update first and last
  • SZS (or 3) - modified method to update last
• (MPM method) - the MPM method used for finding shape functions. The options are:
  • Classic (or Dirac or 0) - classic MPM. This method is the default if no option is specified.
  • uGIMP (or GIMP or 1) - undeformed GIMP.
  • lCPDI (or CPDI or 2) - CPDI using linear shape functions.
  • qCPDI (or 3) - CPDI using quadratic shape functions (available is 2D only)
• (number) - The number of material points in each background grid cell
• (time step) - The MPM time step in ms.
• (time) - The time that the calculations will stop in ms/

• DefGradTerms Command