Material point method

The material point method (MPM) is a numerical technique used to simulate the behavior of solids, liquids, gases, and any other continuum material. Especially, it is a robust spatial discretization method for simulating multi-phase (solid-fluid-gas) interactions. In the MPM, a continuum body is described by a number of small Lagrangian elements referred to as 'material points'. These material points are surrounded by a background mesh/grid that is used to calculate terms such as the deformation gradient. Unlike other mesh-based methods like the finite element method, finite volume method or finite difference method, the MPM is not a mesh based method and is instead categorized as a meshless/meshfree or continuum-based particle method, examples of which are smoothed particle hydrodynamics and peridynamics. Despite the presence of a background mesh, the MPM does not encounter the drawbacks of mesh-based methods (high deformation tangling, advection errors etc.) which makes it a promising and powerful tool in computational mechanics.

The MPM was originally proposed, as an extension of a similar method known as FLIP (a further extension of a method called PIC) to computational solid dynamics, in the early 1990 by Professors Deborah L. Sulsky, Zhen Chen and Howard L. Schreyer at University of New Mexico. After this initial development, the MPM has been further developed both in the national labs as well as the University of New Mexico, Oregon State University, University of Utah and more across the US and the world. Recently the number of institutions researching the MPM has been growing with added popularity and awareness coming from various sources such as the MPM's use in the Disney film Frozen.

An MPM simulation consists of the following stages:

(Prior to the time integration phase)

1. Initialization of grid and material points.
   1. A geometry is discretized into a collection of material points, each with its own material properties and initial conditions (velocity, stress, temperature, etc.)
   2. The grid, being only used to provide a place for gradient calculations is normally made to cover an area large enough to fill the expected extent of computational domain needed for the simulation.

(During the time integration phase - explicit formulation)

1. Material point quantities are extrapolated to grid nodes.
   1. Material point mass ( m m p {\textstyle m_{mp}} ), momenta ( P m p → {\textstyle {\vec {P_{mp}}}} ), stresses ( σ ¯ ¯ m p {\displaystyle {\boldsymbol {\bar {\bar {\sigma }}}}_{mp}} ), and external forces ( b → {\displaystyle {\vec {b}}} ) are extrapolated to the nodes at the corners of the cells within which the material points reside. This is most commonly done using standard linear shape functions ( N n d − m p {\textstyle N_{nd-mp}} ), the same used in FEM.
   2. The grid use the material point values to create the masses ( M n o d e {\textstyle M_{node}} ), velocities ( V n o d e → {\textstyle {\vec {V_{node}}}} ), internal and external force vectors ( F n o d e i n t e r n a l → {\textstyle {\vec {F_{node}^{\mathsf {internal}}}}} , F n o d e e x t e r n a l → {\textstyle {\vec {F_{node}^{\mathsf {external}}}}} ) for the nodes: M n o d e = ∑ m p m m p     N m p − n d {\displaystyle M_{node}=\sum _{mp}m_{mp}~~N_{mp-nd}} V n o d e → = 1 M n o d e     ∑ m p P m p →     N m p − n d {\displaystyle {\vec {V_{node}}}={1 \over M_{node}}~~\sum _{mp}{\vec {P_{mp}}}~~N_{mp-nd}} F n o d e i n t e r n a l → = ∑ m p     σ ¯ ¯ m p     ∇ N m p − n d {\displaystyle {\vec {F_{node}^{\mathsf {internal}}}}=\sum _{mp}~~{\bar {\bar {\sigma }}}_{mp}~~\nabla N_{mp-nd}} F n o d e e x t e r n a l → = ∑ m p b →     N m p − n d {\displaystyle {\vec {F_{node}^{\mathsf {external}}}}=\sum _{mp}{\vec {b}}~~N_{mp-nd}}
2. Equations of motion are solved on the grid.
   1. Newton's 2nd Law is solved to obtain the nodal acceleration ( A n o d e → {\displaystyle {\vec {A_{node}}}} ) A n o d e → = F n o d e e x t e r n a l + F n o d e i n t e r n a l → M n o d e → {\displaystyle {\vec {A_{node}}}={\vec {F_{node}^{external}+{\vec {F_{node}^{internal}}} \over M_{node}}}}
   2. New nodal velocities are found ( V n o d e → ~ {\displaystyle {\tilde {\vec {V_{node}}}}} ). V n o d e → ~ = V n o d e → + A n o d e → d t {\displaystyle {\tilde {\vec {V_{node}}}}={\vec {V_{node}}}+{\vec {A_{node}}}\mathrm {d} t}
3. Derivative terms are extrapolated back to material points
   1. Material point acceleration ( a m p → {\displaystyle {\vec {a_{mp}}}} ), deformation gradient ( F m p ¯ ¯ {\displaystyle {\mathcal {\bar {\bar {F_{mp}}}}}} ) (or strain rate ( ε ˙ m p ¯ ¯ {\displaystyle {\bar {\bar {{\dot {\varepsilon }}_{mp}}}}} ) depending on the strain theory used) is extrapolated from the surrounding nodes using similar shape functions to before ( N n d − m p {\displaystyle N_{nd-mp}} ). a m p → = ∑ n d A n o d e →     N n d − m p {\displaystyle {\vec {a_{mp}}}=\sum _{nd}{\vec {A_{node}}}~~N_{nd-mp}} ε ˙ m p ¯ ¯ = ∑ n d   1 2     [ V n o d e → ∇ N n d − m p + ( V n o d e ∇ N n d − m p ) T ] {\displaystyle {\bar {\bar {{\dot {\varepsilon }}_{mp}}}}=\sum _{nd}~{1 \over 2}~~[{\vec {V_{node}}}\nabla N_{nd-mp}+(V_{node}\nabla N_{nd-mp})^{T}]}
   2. Variables on the material points: positions, velocities, strains, stresses etc. are then updated with these rates depending on integration scheme of choice and a suitable constitutive model.
4. Resetting of grid.

   Now that the material points are fully updated at the next time step, the grid is reset to allow for the next time step to begin.

The PIC was originally conceived to solve problems in fluid dynamics, and developed by Harlow at Los Alamos National Laboratory in 1957.^[1] One of the first PIC codes was the Fluid-Implicit Particle (FLIP) program, which was created by Brackbill in 1986^[2] and has been constantly in development ever since. Until the 1990s, the PIC method was used principally in fluid dynamics.

Motivated by the need for better simulating penetration problems in solid dynamics, Sulsky, Chen and Schreyer started in 1993 to reformulate the PIC and develop the MPM, with funding from Sandia National Laboratories.^[3] The original MPM was then further extended by Bardenhagen et al.. to include frictional contact,^[4] which enabled the simulation of granular flow,^[5] and by Nairn to include explicit cracks^[6] and crack propagation (known as CRAMP).

Recently, an MPM implementation based on a micro-polar Cosserat continuum^[7] has been used to simulate high-shear granular flow, such as silo discharge. MPM's uses were further extended into Geotechnical engineering with the recent development of a quasi-static, implicit MPM solver which provides numerically stable analyses of large-deformation problems in Soil mechanics.^[8]

Annual workshops on the use of MPM are held at various locations in the United States. The Fifth MPM Workshop was held at Oregon State University, in Corvallis, OR, on April 2 and 3, 2009.

The uses of the PIC or MPM method can be divided into two broad categories: firstly, there are many applications involving fluid dynamics, plasma physics, magnetohydrodynamics, and multiphase applications. The second category of applications comprises problems in solid mechanics.

The PIC method has been used to simulate a wide range of fluid-solid interactions, including sea ice dynamics,^[9] penetration of biological soft tissues,^[10] fragmentation of gas-filled canisters,^[11] dispersion of atmospheric pollutants,^[12] multiscale simulations coupling molecular dynamics with MPM,^[13][14] and fluid-membrane interactions.^[15] In addition, the PIC-based FLIP code has been applied in magnetohydrodynamics and plasma processing tools, and simulations in astrophysics and free-surface flow.^[16]

As a result of a joint effort between UCLA's mathematics department and Walt Disney Animation Studios, MPM was successfully used to simulate snow in the 2013 animated film Frozen.^[17][18][19]

MPM has also been used extensively in solid mechanics, to simulate impact, penetration, collision and rebound, as well as crack propagation.^[20][21] MPM has also become a widely used method within the field of soil mechanics: it has been used to simulate granular flow, quickness test of sensitive clays,^[22] landslides,^[23][24][25] silo discharge, pile driving, fall-cone test,^{[26][27][28][29]} bucket filling, and material failure; and to model soil stress distribution,^[30] compaction, and hardening. It is now being used in wood mechanics problems such as simulations of transverse compression on the cellular level including cell wall contact.^[31] The work also received the George Marra Award for paper of the year from the Society of Wood Science and Technology.^[32]

One subset of numerical methods are Meshfree methods, which are defined as methods for which "a predefined mesh is not necessary, at least in field variable interpolation". Ideally, a meshfree method does not make use of a mesh "throughout the process of solving the problem governed by partial differential equations, on a given arbitrary domain, subject to all kinds of boundary conditions," although existing methods are not ideal and fail in at least one of these respects. Meshless methods, which are also sometimes called particle methods, share a "common feature that the history of state variables is traced at points (particles) which are not connected with any element mesh, the distortion of which is a source of numerical difficulties." As can be seen by these varying interpretations, some scientists consider MPM to be a meshless method, while others do not. All agree, however, that MPM is a particle method.

The Arbitrary Lagrangian Eulerian (ALE) methods form another subset of numerical methods which includes MPM. Purely Lagrangian methods employ a framework in which a space is discretised into initial subvolumes, whose flowpaths are then charted over time. Purely Eulerian methods, on the other hand, employ a framework in which the motion of material is described relative to a mesh that remains fixed in space throughout the calculation. As the name indicates, ALE methods combine Lagrangian and Eulerian frames of reference.

PIC methods may be based on either the strong form collocation or a weak form discretisation of the underlying partial differential equation (PDE). Those based on the strong form are properly referred to as finite-volume PIC methods. Those based on the weak form discretisation of PDEs may be called either PIC or MPM.

MPM solvers can model problems in one, two, or three spatial dimensions, and can also model axisymmetric problems. MPM can be implemented to solve either quasi-static or dynamic equations of motion, depending on the type of problem that is to be modeled. Several versions of MPM include Generalized Interpolation Material Point Method ^[33];Convected Particle Domain Interpolation Method;^[34] Convected Particle Least Squares Interpolation Method.^[35]

The time-integration used for MPM may be either explicit or implicit. The advantage to implicit integration is guaranteed stability, even for large timesteps. On the other hand, explicit integration runs much faster and is easier to implement.

Unlike FEM, MPM does not require periodical remeshing steps and remapping of state variables, and is therefore better suited to the modeling of large material deformations. In MPM, particles and not the mesh points store all the information on the state of the calculation. Therefore, no numerical error results from the mesh returning to its original position after each calculation cycle, and no remeshing algorithm is required.

The particle basis of MPM allows it to treat crack propagation and other discontinuities better than FEM, which is known to impose the mesh orientation on crack propagation in a material. Also, particle methods are better at handling history-dependent constitutive models.

Because in MPM nodes remain fixed on a regular grid, the calculation of gradients is trivial.

In simulations with two or more phases it is rather easy to detect contact between entities, as particles can interact via the grid with other particles in the same body, with other solid bodies, and with fluids.

MPM is more expensive in terms of storage than other methods, as MPM makes use of mesh as well as particle data. MPM is more computationally expensive than FEM, as the grid must be reset at the end of each MPM calculation step and reinitialised at the beginning of the following step. Spurious oscillation may occur as particles cross the boundaries of the mesh in MPM, although this effect can be minimized by using generalized interpolation methods (GIMP). In MPM as in FEM, the size and orientation of the mesh can impact the results of a calculation: for example, in MPM, strain localisation is known to be particularly sensitive to mesh refinement. One stability problem in MPM that does not occur in FEM is the cell-crossing errors and null-space errors^[36] because the number of integration points (material points) does not remain constant in a cell.

• Center for Simulation of Accidental Fires and Explosions – MPM code available
• NairnMPM – open source
• MPM3D - open source (MPM3D-F90) and free trial version (MPM3D)
• Taichi - Physically Based Computer Graphics Library – open source MPM code available
• Anura3D open source – software for geotechnical problems and soil-water-structure interactions by Anura3D MPM Research Community