MPM Methods and Simulation Timing

The input commands described below set details about the MPM method and control the simulation timing. Before describing the commands, the first few sections give background theory about these options.

Theory: Shape Functions

Many tasks in MPM involve extrapolations from particles to the grid or from the grid to 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 characteristic function for particle p, and N_i(x) is grid shape function for node i. The various MPM methods depend choices for χ_p(x) and N_i(x) and 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 NairnMPM.

File:MPMTree.png

Particle Shape Functions

One choice is for χ_p(x) or the particle shape functions; 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.

Spline Grid Shape Functions

A second choice is for grid shape functions can be derived Bⁿ spline basis functions, which are C^(n-1) continuous spline functions on the grid. Most MPM methods choose B¹ order spline functions or linear shape functions. In 1D, linear shape functions have the form:

      [math]\displaystyle{ N_i(x) = \left\{ \begin{array}{cc} 1+{x-x_i\over \Delta x} & {\rm for\ } x_i-\Delta x \lt x \lt x_i \\ 1-{x-x_i\over \Delta x} & {\rm for\ } x_i \le x \lt x_i+\Delta x \\ 0 & {\rm otherwise} \end{array} \right. }[/math]

where x_i is coordinate of node i and Δx is cell size (assumed same on both sides of node i).

Another choice supported in this code is to choose B² spline functions or quadratic spline functions. In 1D, quadratic spline shape functions have the form:

      [math]\displaystyle{ N_i(x) = \left\{ \begin{array}{cc} {1\over 8}\left(3+{2(x-x_i)\over \Delta x}\right)^2 & {\rm for\ } x_i-{3\Delta x\over 2} \lt x \lt x_i-{\Delta x\over 2} \\ {3\over 4} - \left({x-x_i\over \Delta x}\right)^2 & {\rm for\ } x_i-{\Delta x\over 2} \le x \le x_i+{\Delta x\over 2} \\ {1\over 8}\left(3-{2(x-x_i)\over \Delta x}\right)^2 & {\rm for\ } x_i+{\Delta x\over 2} \lt x \lt x_i+{3\Delta x\over 2} \\ 0 & {\rm otherwise} \end{array} \right. }[/math]

Classic MPM Methods

If χ_p(x) = δ(x), the shape functions reduce to the grid shape function (e.g., S_ip = N_i(x_p)), which recovers "Classic MPM" or the original derivation of MPM. (when using linear N_i(x)) Note that "Classic GIMP" predates GIMP even though it descends from GIMP in the MPM "genealogy." The "Dirac" method options are

• Classic (or Dirac) chooses "Classic MPM" with linear N_i(x) on the grid.
• B2SPLINE chooses Dirac particle functions with quadratic spline N_i(x) on the grid.

Although Classic MPM is the most efficient method, it is rarely used because it is not very robust for large displacement problems. Whenever particles 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 B2SPLINE (because first derivatives are continuous) and/or by choosing a 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). Many MPM methods descend from this selection. The specific descendant methods depend on how they integrate over the particle domain.

Undeformed GIMP

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 "Undeformed GIMP", has the greatest efficiency because shape functions can be explicitly calculated for all possible particle locations relative to a regular grid. Undeformed GIMP works well from small to fairly large deformation and greatly reduces cell-crossing artifacts. The two undeformed GIMP methods can be selected using

• uGIMP or undeformed GIMP with linear N_i(x) on the grid.
• B2GIMP or undefomed GIMP with quadratic spline N_i(x) on the grid.

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.

Deformed GIMP

The following GIMP methods attempt to account for particle deformation in evaluation of the shape function integrals. The methods differ in how they do the domain integration and in the grid shape functions that are integrated.

1. Finite GIMP

In finite GIMP, the shape function integrals are fully integrated over the deformed particle domain. This option can be selected using Finite for shape function method. This exact integration is available only in OSParticulas, only for 2D calculations, and only for linear grid shape functions. The full integration compared to approximate integral is similar in efficiency, but surprisingly, the results are not improved over other methods. The reasons are not yet clear. One hypothesis is that while finite GIMP exactly integrates the full particle domain, the particle shape function is still treating particle properties (such as velocity) as constant over the domain (due to Heaviside χ_p(x) function). These contrary views may be inconsistent.

2. Convected Partical Domain Integration

A useful improvement for large to massive deformation is Convected Particle Domain Integration approach or CPDI. In CPDI, the integral over the deformed particle domain is approximated. 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).

The version of CPDI available here 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_j0, 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 grid shape functions being used and on the finite element types used to map that particle domain:

• lCPDI for "linear" CPDI: this method uses linear grid shape functions and 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 uses linear grid shape functions and additionally maps midpoints of edges to nodes (e.g., 8 node element in 2D) resulting in an isoparametric element with quadratic shape functions. This method is only available for 2D simulations.
• B2CPDI for "B2 Spline" CPDI: this method uses B₂ spline grid shapes functions and maps the deformed domain as in lCPDI.

CPDI methods are less efficient the undeformed GIMP because of the extra overhead required to locate corner nodes needed for evaluating shape functions. Nevertheless, the use of CPDI is crucial for handling problems with deformation beyond the capabilities of undeformed GIMP.

3. Convected Particle Domain Integration 2

A second version of CPDI, known as CPDI2, has been developed. 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. The efficiency should be similar to CPDI, but CPDI2 adds the additional overheard of tracking and updating corner particles. CPDI2 options are currently not implemented in either NairnMPM or OSParticulas.

4. cpGIMP

In cpGIMP, or convected particle GIMP, 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). In 1D, cpGIMP is a true finite GIMP method, but in 2D or 3D, it does not account for arbitrary deformation that includes shear deformation. Compared to undeformed GIMP, cpGIMP may handle uniaxial deformation better, but may not be any better for general deformation modes. This approach is currently not implemented in either NairnMPM or OSParticulas.

Theory: Stress and Strain Updates

The basic tasks in an MPM time step are:

1. Extrapolate particle mass and momentum to the grid
2. Find internal and external forces
3. Update momenta on the grid using those forces
4. Update particle position and velocity using updated grid information

To correctly model internal forces, this loop needs to track stresses and strains on the particles. This tracking is done by calculating velocities on the grid (from the grid momenta). Once these velocities are found, they can be extrapolated to the particle using gradient shape functions to get velocity gradients on the particle. Particle velocity gradients provide the key input for updating particle strains and through constitutive laws to update particle stresses. In the first MPM work, the strain and stress update was done using the updated momenta on the grid. It is easy to show, however, that this approach is unstable if the update particle position interacts with a node having zero mass. This numerical challenge can be eliminated by two variants to MPM:

• USL - follow the recommendation of Sulsky, Zhou, and Schreyer which was to do a second extrapolation of momenta after updating momenta on the grid and use those new momenta to update strains and stresses. This approach adds step 5 to MPM time step above to reextrapolate momenta to the grid and then update stress and strain after that step.
• USF - update strains and stresses first or to update them after the initial extrapolation in step 1 of above time step tasks.

These two approaches give similar results. The USF method is more efficient than USL because it avoids the need to do a second extrapolation of momenta.

A unique option in NairnMPM combines these two methods:

• USAVG - update strains after the first extrapolation and after a second extrapolation of updated momenta on the grid.

The USAVG method seems to behave like a midpoint rule. In observations of temporal convergence, the USAVG method for a given time step gives results close to USF or USL using half that time step. In other words, compared to USF and USL, the USAVG is analogous to using half the time step without the cost of doubling the number of calculations. In most simulation, the USAVG method also exhibits improved energy conservation.

Although Sulsky, Zhou, and Schreyer recommend a second momentum extrapolation prior to update, that second extrapolation is not needed provided one is careful in calculation of shape functions. Skipping the second extrapolation is always more efficient and can result in better simulations. Whether or not to use the second extrapolation, however, seems to depend on the problem. Some problems with contact or cracks seem to require it while problems with liquids seem improved without it. The default is to include the second extrapolation, but commands below let you pick from the following options:

• USL+ - update strains after momentum update and after a second extrapolation (same as USL above).
• USL- - update strains after momentum update, but do not do a second extrapolation.
• USAVG+ - update strains after the first extrapolation and after a second extrapolation of updated momenta on the grid (same as USAVG above).
• USAVG- - update strains before and after the momentum update, but do not do a second extrapolation.

The key to skipping the second extrapolation compared to early MPM work, is to extrapolate velocity gradient from updated momenta by using gradient shape function based on the original particle position rather then the updated position.

Incremental Deformation Gradient

During strain updating, some materials (e.g., hyperelastic materials) use the full incremental deformation gradient, dF, which is defined by

      [math]\displaystyle{ {dF\over dt} = \nabla\vec v F }[/math]

where Δt is the time step and ∇v is the velocity gradient calculated from the current nodal velocities. If velocity gradient is constant for the time step, the exact solution is

      [math]\displaystyle{ F(t+\Delta t) = \exp(\Delta t \nabla\vec v) F(t) = dF\cdot F(t) }[/math]

where

      [math]\displaystyle{ dF = \exp(\Delta t \nabla\vec v) }[/math]

is the incremental deformation gradient. An interesting review article on finding the exponential of a matrix, which is needed here, concludes no single way is stable and efficient for all problems, but in MPM with sufficiently small time steps, an expansion method should work well

      [math]\displaystyle{ \exp(\Delta t \nabla\vec v) = \sum_{k=0}^{k_{max}} {(\nabla u)^k \over k!} = I + \nabla u + \sum_{k=2}^{k_{max}} {(\nabla u)^k \over k!} }[/math]

where ∇u = Δt ∇v is the incremental displacement gradient for the current time step. Many computational mechanics methods truncate after the first term (the linear term), because it avoids the need to take powers of a matrix. But, using just the linear term can have issues for problems with lots of rotational motion. Fortunately the extra terms can all be evaluated with no matrix multiplications (in 2D) and with only the square of a matrix (in 3D) by using the "Cayley-Hamilton" theory (see method #8 in the review article). You can select the number of terms to use in the expansion. The default is 2, which seems to work well for almost all simulations and to reduce errors caused by large rotations.

Theory: MPM Time Step

In explicit computation mechanics, the Courant–Friedrichs–Lewy condition (CFL condition) is a necessary condition for convergence. It states that the time step, Δt, must satisfy:

      [math]\displaystyle{ \Delta t \le { \Delta x\over v} }[/math]

where Δx is the minimum dimension of any cell in the background grid and v is the maximum wave speed in the material (i.e., the highest wave speed in anisotropic materials). It is usually better to be below this minimum by a input factor of less than one. The factor is also known as the "Courant" factor:

      [math]\displaystyle{ C = { v\Delta t\over \Delta x} \qquad {\rm or} \qquad\Delta t =C { \Delta x\over v} }[/math]

The entered C must be less then 1.0. The default in NairnMPM is C = 0.5. Lower values (sometimes much loweer) may be are needed to resolve instabilities around fixed-velocity boundary conditions on the grid or created by rigid particles when modeling contact, and when modeling brittle materials using damage mechanics.

When transport tasks are activated (conduction, diffusion, or poroelasticity), convergence of transport calculations requires the time step to satisfy

      [math]\displaystyle{ \Delta t_T \le { \Delta x^2\over 2}{C_T\over \max(\kappa)} }[/math]

where C_T and [math]\displaystyle{ \kappa }[/math] depend on the type of transport being modeled. For conduction, [math]\displaystyle{ C_T=\rho C_v }[/math] and [math]\displaystyle{ \kappa }[/math] is thermal conductivity. For diffusion, C_T=1 and [math]\displaystyle{ \kappa }[/math] is diffusion coefficient. For poroelasticity, C_T=1/Q and [math]\displaystyle{ \kappa }[/math] is Darcy's law permittivity divided by fluid viscosity. It is usually better to be below this minimum by a input factor of less than one. The factor is also known as the "Transport Courant" factor or transport time step is:

      [math]\displaystyle{ \Delta t_T = F { \Delta x^2\over 2}{C_T\over \max(\kappa)} }[/math]

The entered F must be less then 1.0. The default in NairnMPM is F = 0.5.

When running coupled calculations, the simulation time step must use the minimum of mechanics and transport time steps or:

      [math]\displaystyle{ \Delta t_{final} = \min(\Delta t,\Delta t_T) }[/math]

For realistic properties, the minimum is usually [math]\displaystyle{ \Delta t }[/math] by several orders of magnitude, which means entered value for F has no effect. Notice, however, that transport time step scales with [math]\displaystyle{ \Delta x^2 }[/math] while mechanics time step is linear in [math]\displaystyle{ \Delta x }[/math]. The minimum value could change if using very small cells.

Input Commands

The MPM method, the strain updating method, the number of points per element, and simulation timing are set with one or more of the commands in this section. In script input files, the commands are:

MPMMethod (strainUpdate),<(MPM_method)>
PtsPerElement (axisNumber)
TimeStep (timeStep),<(maxTime)>,<(timeFactor)>
MaximumTime (maxTime)
CFLFactor (timeFactor),<(transTimeFactor)>

In XML input files, the settings are made with one or more of the following commands (all of which must be within the <MPMHeader> element):

<MPMMethod>(strainUpdate number)</MPMMethod>
<SkipPostExtrapolation/>
<GIMP type='(MPM_method)'/>
<MatlPtsPerElement>(totalNumber)</MatlPtsPerElement>
<TimeStep>(timeStep)</TimeStep>
<TimeFactor>(timeFactor)</TimeFactor>
<TransTimeFactor>(transTimeFactor)</TimeFactor>
<MaxTime>(maxTime)</MaxTime>
<DefGradTerms>(kmax)</DefGradTerms>

where

• (strainUpdate) - the method used to update stresses and strains on the particles. The options are (scripted files can use name or number while XML files must use number):
  • USF (or 0) - update first
  • USAVG or USAVG+ (or 2) - update strain first and last. This option is the default option.
  • USAVG- (no number) - update strain first and last, but last is modified¹. Because there is no number option, XML input files should add the <SkipPostExtraolation> command to get this method.
  • USL or USL+ (or 3) - update strain last (SZS is old name for this method that is still accepted)
  • USL- (no number) - update strain last, using modified¹ method. Because there is no number option, XML input files should add the <SkipPostExtraolation> command to get this method.
• (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 or if no <GIMP> command is used in an XML input file.
  • uGIMP (or GIMP or 1)^2,3 - undeformed GIMP. If <GIMP> command is used in XML files without a type attribute, this option is used.
  • lCPDI (or CPDI or 2)³ - CPDI using linear shape functions.
  • qCPDI (or 3)³ - CPDI using quadratic shape functions (available in 2D only).⁵
  • Finite (or 4)^2,3 - exactly integrated finite GIMP (available in 2D only and only in OSParticulas).
  • B2GIMP (or 5)^2,4 - undeformed GIMP methods using quadratic spline grid shape functions.
  • B2SPLINE (or 6)^2,4 - Classic method using quadratic spline grid shape functions.
  • B2CPDI (or 7)^2,4 - CPDI methods using quadratic spline grid shape functions.
• (axisNumber) (scripted files only) - the number of material points along each axis in the grid.⁶ The options are
  • 1 to 5 for 2D - these settings results in 1, 4, 9, 16, or 25 points per element. The default is 2 or 4 points per element.
  • 1 to 3 for 3D - these settings results in 1, 8 or 27 points per element. The default is 2 or 8 points per element.
• (totalNumber) (XML files only) - the number of material points in each element of the background grid.⁶ The options are
  • 1, 4, 9, 16, or 25 for 2D. The default is 4 points per element.
  • 1, 8 or 27 for 3D. The default is 8 points per element.
• (timeStep) - the MPM time step in alt time units (in XML files, the time is in time units or determined by a units attribute). This input will be ignored unless it is less that time step calculated for input (or default) values for (timeFactor) and (transTimeFactor).
• (maxTime) - the time that the calculations will stop in alt time units (in XML files, the time is in time units or determined by a units attribute). In scripted files, this value can be entered as second argument to the TimeStep command or in a separate MaximumTime command (whichever is last will be the one that is used) (it can be an entity).
• (timeFactor) and (transTimeFactor) - the Courant factors to use when determining the MPM time step. The time step used in the analysis will be the minimum of the one determined using these factors and the one entered in (timeStep). The default Courant factors are 0.5. In scripted files, The (timeFactor) can be entered as third argument to the TimeStep command or first argument in a CFLFactor command
• (kmax) - the number of terms to use when evaluating the incremental deformation gradient.⁷ The default is 2.

Notes

1. The USAVG+ and USL+ methods re-extrapolate momenta to the grid and re-impose contact conditions prior to updating the strain after the momenta update. The USAVG- and USL- methods skip the second extrapolation. The method to use may depend on the problem, but some problems with contact seem to require the "+" methods. If both are equal, the "-" methods are more efficient.
2. uGIMP, Finite, B2GIMP, B2SPLINE, and B2CPDI all require use of a regular grid (i.e., a grid with constant size and orthogonal elements).
3. These methods are based on linear N_i(x) grid shape functions.
4. These methods are based on quadratic spline N_i(x) grid shape functions.
5. qCPDI is less efficient than lCPDI and in trial runs appears to have little or no benefit. It is available for comparison and potential future research into MPM shape function methods.
6. The number of points per element must be set before defining material points in your analysis. Once you start defining material points, it cannot be changed. In other words elements with different numbers of material points cannot be used in a single analysis. Note that script files and XML files set this number differently. Script files set the number for each axis (which is squared for 2D or cubed for 3D) while XML files gives the total number for the cells.
7. Because two terms for the incremental deformation gradient seems good for all calculations, there is currently no script command to change this setting. If necessary, you can change it using an XMLData Command:

   XMLData "MPMHeader"
       <DefGradTerms>(kmax)</DefGradTerms>
   EndXMLData