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_pi) and shape function gradient (G_pi) for extrapolation from node i to particle p are:

      [math]\displaystyle{ S_{pi} = {\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_{pi} = {\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 characteristic function for particle p, and N_i(x) is grid shape function for node i. The various MPM methods depend on 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.

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 from 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 NairnMPM is to use 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 (do not use yet for axisymmetric MPM).

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 are reduced by using B2SPLINE (because first derivatives are continuous), but problems often still persist. The better approach to handling cell cross artifacts in to choose a different χ_p(x) (e.g., the Heaviside methods).

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 functions, 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 (do not use yet for axisymmetric MPM).

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_pi and G_pi 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_pi and G_pi 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 with the <DefGradTerms> command.

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 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. An alternative to very low C is to try restarting time steps that develop high accelerations on the grid.

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.

Restarting Time Steps

Some problems might develop high accelerations on the grid that cause simulations to break down (such as particles leaving the grid). The first steps to fixing such problems should be:

1. Use a smaller time step C factor
2. Add an AdjustTimeStep Custom Task
3. Increase number of terms used to find the incremental deformation gradient - this option works for all materials except small-strain materials when they are not using the largeRotation=1 option.

If none of these work, you can restart time steps as needed based on nodal velocities and accelerations. Each MPM time step extrapolates mass (m), momentum (p), and force (f) to the grid. When time-step restarting is activated, each step checks the virtual nodal displacement on each node implied by extrapolated nodal values:

      [math]\displaystyle{ d = \| \vec p + {1\over2}\vec f\Delta t\|{\Delta t\over m} }[/math]

where Δt is the current time step. If any node's d relative to its cell size is greater than an entered rescaling factor, the time step is discarded and restarted using a smaller time step. This option is invoked using the <RestartScaling> command. To prevent the time step from permanently remaining small (and making a simulation too slow), this option should always be used in conjunction with an AdjustTimeStep Custom Task. This tasks will allow the time step to increase once nodal accelerations get under control, but should include its maxIncrease parameter to prevent returning immediately to the previous time step.

Input Commands

The MPM shape functions, 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)>
CPDIrcrit (rcrit)

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)'/>
<CPDIrcrit>(rcrit)</CPDIrcrit>
<MatlPtsPerElement>(totalNumber)</MatlPtsPerElement>
<TimeStep>(timeStep)</TimeStep>
<TimeFactor>(timeFactor)</TimeFactor>
<TransTimeFactor>(transTimeFactor)</TimeFactor>
<MaxTime>(maxTime)</MaxTime>
<DefGradTerms>(kmax)</DefGradTerms>
<RetartScaling CFL='(restartFactor)'>(restartScaling)</RestartScaling>

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.
• (rcrit) - one problem with CPDI shape functions is that if a particle deforms too much, its' corners may extend beyond ghost rows into other patches used in parallel coding. If that becomes a problem, this parameter will restrict corner deformation to keep them within a single patch. Enter maximum allowed particle radius in units of cell size. Truncation occurs when half of any diagonal of the deformed particle exceeds (rcrit). The particle domain truncation implements the methods described in Homel et al. (see Eqs. (4) to (11)).
• (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 for 2D and 1 for 3D.
• (restartFactor) and (restartScaling) - these parameters invoke the time-step restarting option. If the implied displacement of any node is greater than (restartFactor)*(cell size), the time step is changed to |(restartScaling)|**Δt where Δt is the current time step. The entered (restartFactor) should be between -1 and 1. The scaling uses absolute value of the entered value. Negative values are a "verbose" mode that prints a warning eah time the time step is restarted; postive values restart time steps without notice. Values out side the range -1 to 1 are change to 0.5 (preserving the sign). A report at the end of a simulations summarizes the number of times steps that had to be restarted.⁸

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. The default number of terms for the incremental deformation gradient is 2 for 2D simulations and 1 for 3D simulations (the latter for efficieny). Because these settings work well for most calculations, there is currently no script command to change this setting. If necessary, you use an XMLData Command:

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

8. Because restarting time steps is uncommon (and often does not fix the modeling problems), there is currently no script command to invoke it. If necessary, you can use an XMLData Command:

   XMLData "MPMHeader"
        <RetartScaling CFL='(restartFactor)'>(restartScaling)</RestartScaling>
   EndXMLData