MPM Methods and Simulation Timing

From OSUPDOCS
Jump to navigation Jump to search

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[1]), the shape function (Sip) and shape function gradient (Gip) 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 Ni(x) is grid shape function for node i. The various MPM methods depend on the technique used for evaluating these 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.

The first decision is to choose χ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 grid shape function (e.g., Sip = Ni(xp)), which recovers "Classic MPM" or the original derivation of MPM.[2][3][4] 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 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 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). 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[1]), it is too inefficient for arbitrary 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.[5] 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 Sip and Gip 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 cj are nodes of the deformed particle domain (e.g., usually corners of the domain) and Mj(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 Sip and Gip in terms for Ni(cj).

Two versions of CPDI have appeared. The first version finds the corners of the deformed domain using the deformation gradient on the particle.[5] If rj0 is a vector from the particle to any edge location j in the undeformed domain, then rj = Fprj0, where Fp 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[6] 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 adds 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, 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 uGIMP, cpGIMP may handle uniaxial deformation better, but may not be any better for general deformation modes.

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 loops need to track stresses and strains on the particles. This tracking is done by calculating velocities on the grid (from the gird momenta). Once these velocities are found, it is possible to get velocity gradients on the particle, which are sufficient to update particle strains and through constitutive laws to update particle stresses. In the first MPM work,[2] 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 whenever an edge particle is near a node. It is unstable because of a potential division by a near-zero number. This numerical challenge is eliminated (in a careful algorithm) by updating strains and stresses right after extraplated momenta rather than after updating momenta. This observation led to two variants of MPM for stress and strain updates:

  • USF - to update strains and stresses first or to update them after the initial extrapolation in the above time step tasks.[7]
  • SZS - to follow the recommendation of Sulsky, Zhou, and Schreyer[3] 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.

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

A unique option in NairnMPM combines these two methods:

  • USAVG - to update strains both after the first extrapolation and after a second extrapolation of updated momenta to the gird.

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 SZS using half that time step. In other words, compared to USF and SZS, the USAVG is analogous to using half the time step without the cost of doubling the number of calculations. In one impact simulation, the USAVG method also exhibited improved energy conservation[8] (although this result was with classic MPM methods and may not be as dramatic with newer, improved MPM methods).

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[9] 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[9]). 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 greatly reduce errors caused by 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 are sometimes needed and commonly are needed to resolve instabilities around fixed-velocity boundary conditions on the grid or created by rigid particles.

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)

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>
<GIMP type='(MPM_method)'/>
<MatlPtsPerElement>(totalNumber)</MatlPtsPerElement>
<TimeStep>(timeStep)</TimeStep>
<TimeFactor>(timeFactor)</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 2) - update first and last. This option is the default option.
    • 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 or if no <GIMP> command is used in an XML input file.
    • uGIMP (or GIMP or 1)1 - undeformed GIMP. If <GIMP> command is used in XML files without a type attribute, this option is used.
    • lCPDI (or CPDI or 2)1 - CPDI using linear shape functions.
    • qCPDI (or 3)1 - CPDI using quadratic shape functions (available in 2D only).3
  • (axisNumber) (scripted files only) - the number of material points along each axis in the grid.2 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.2 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 ms (in XML files, the time is in sec or determined by a units attribute)
  • (timeFactor) - the Courant factor to use when determining the MPM time step. The time step used in the analysis will be the minimum of the one determined using this factor and the one entered in (timeStep). The default Courant factor is 0.5.
  • (maxTime) - the time that the calculations will stop in ms (in XML files, the time is in sec 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).
  • (kmax) - the number of terms to use when evaluating the incremental deformation gradient.4 The default is 2.

Notes

  1. uGIMP, lCPDI, and qCPDI all require use of a regular grid (i.e., a grid with constant size and orthogonal elements).
  2. 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.
  3. 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.
  4. 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
    

References

  1. 1.0 1.1 S. G. Bardenhagen and E. M. Kober, "The Generalized Interpolation Material Point Method," Computer Modeling in Engineering & Sciences, 5, 477-496 (2004).
  2. 2.0 2.1 D. Sulsky, Z. Chen, and H. L. Schreyer, "A Particle Method for History-Dependent Materials," Comput. Methods Appl. Mech. Engrg., 118, 179-186 (1994).
  3. 3.0 3.1 D. Sulsky, S. -J. Zhou, and H. L. Schreyer, "Application of a Particle-in-Cell Method to Solid Mechanics," Comput. Phys. Commun., 87, 236-252 (1995).
  4. S. -J. Zhou, "The Numerical Prediction of Material Failure Based on the Material Point Method," Ph.D. Thesis, University of Mexico (1998).
  5. 5.0 5.1 A. Sadeghirad, R. M. Brannon, and J. Burghardt, "A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations," Int. J. Numer. Meth. Engng., 86, 1435–1456 (2011)
  6. A. Sadeghirad, R. M. Brannon,and J. E. Guilkey, "Second-order convected particle domain interpolation (CPDI2) with enrichment for weak discontinuities at material interfaces," Int. J. Numer. Meth. Engng., 95, 928–952 (2013)
  7. S. G. Bardenhagen, "Energy Conservation Error in the Material Point Method," J. Comp. Phys., 180, 383-403 (2002).
  8. J. A. Nairn, "Material Point Method Calculations with Explicit Cracks," Computer Modeling in Engineering & Sciences, 4, 649-664 (2003). (See PDF)
  9. 9.0 9.1 C. Moler and C. Van Loan, "Nineteen Dubious Ways to Compute the Exponential of a Matrix. Twenty Five Years Later." SIAM Review, 46, (2003)