See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/262415477

Material point method: basics and applications

Data · May 2014

CITATIONS READS

6 7,911

1 author:

Vinh Phu Nguyen

Monash University (Australia)
112 PUBLICATIONS   3,333 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Computational modelling of crack propagation in solids View project

Hydraulic fracturing modelling View project

All content following this page was uploaded by Vinh Phu Nguyen on 20 May 2014.

The user has requested enhancement of the downloaded file.

 Material point method: basics and applications
Vinh Phu Nguyen
Cardiff Unitversity
Department of Civil Engineering
Institute of Advanced Mechanics and Materials

Contents
1 Lagrangian, Eulerian and Arbitrary Lagrangian-Eulerian descriptions 2

2 Introduction 3

3 Governing equations 5

4 Weak form and discretisation 5

5 Explicit time integration 7

6 Implementation 10
6.1 Eulerian grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
6.2 Shape functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
6.2.1 B-splines basis functions . . . . . . . . . . . . . . . . . . . . . . . . . 16
6.3 Initial particle distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
6.4 Visualisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

7 Examples 18
7.1 1D examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
7.2 Axial vibration of a continuum bar . . . . . . . . . . . . . . . . . . . . . . . . 22
7.3 Impact of two elastic bodies . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
7.4 Impact of two elastic rings . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

8 Generalized Interpolation Material Point Method-GIMP 32

9 Implicit time integration 34

10 Fluid-structure interaction 36

1
 1 LAGRANGIAN, EULERIAN AND ARBITRARY LAGRANGIAN-EULERIAN
DESCRIPTIONS

Figure 1: Eulerian vs Lagrangian description.

1 Lagrangian, Eulerian and Arbitrary Lagrangian-Eulerian

descriptions
In Eulerian formulations the governing equations are solved on the fixed grid, and material
moves through the mesh. In Lagrangian methods, the computational grid or mesh is attached to
the material being simulated. The Lagrangian mesh moves and distorts with the material.
From Fig. 1 it is obvious that the most distinct advantage of the Eulerian description is
the ability to deal with highly distorted motion since the grid is always fixed. However it suf-
fers from the following disadvantages (1) difficulties with history dependent materials because
material points are not tracked, (2) difficulties in defining material boundaries. The stress of
fluids is independent of its history, the Eulerian description is then widely adopted in describing
the motion of fluids. Furthermore, the material time derivative in an Eulerian description is
the sum of a spatial time derivative and a convective term (this convective term is vanished in
a Lagrangian formulation). This term is generally expensive to handle from a computational
perspective.
There are also some mixed methods which combine the advantages of these two descriptions
and to avoid their disadvantages, such as the arbitrary Lagrangian-Eulerian (ALE). In ALE,
the computational mesh is continuously moved independently of the material deformation to
optimize element shapes and describe the boundaries accurately. However, the convection term
still exists in the ALE formulation, which may cause numerical difficulties. In addition, it

2
2 INTRODUCTION

still remains a challenging task to design an efficient and effective mesh moving algorithm to
maintain mesh regularity for three-dimensional complicated material domain. The Material
Point Method (MPM) was born as a simple ALE method.

2 Introduction
The Material Point Method (MPM) is one of the latest developments in particle-in-cell (PIC)
methods. The first PIC technique was developed in the early 1950s and was used primarily
for applications in fluid mechanics. Early implementations suffered from excessive energy
dissipation which was overcome in 1986, by Brackbill and Ruppel and they introduced FLIP-the
Fluid Implicit Particle method. The FLIP technique was modified and tailored for applications
in solid mechanics by Sulsky and her co-workers Sulsky et al. [1994, 1995] and has since been
referred to as the material point method (MPM) Sulsky and Schreyer [1996].
The MPM discretizes a continuum body Ω with a finite set of np material points (or par-
ticles) in the original configuration that are tracked throughout the deformation process. The
terms particle and material point will be used interchangeably throughout this manuscript. The
mass and volume of subregions which the particles represent are memorised for these points,
but changes in the shape of the subregions are not traced1 . Let xtp (p = 1, 2, . . . , np ) denote the
current position of material point p at time t. Each material point has an associated mass Mp ,
density ρp , velocity vp , Cauchy stress tensor σp and any other internal state variables necessary
for the constitutive model. Thus, these material points provide a Lagrangian description of the
continuum body. Since each material point contains a fixed amount of mass for all time, mass
conservation is automatically satisfied. In MPM, a fixed Eulerian grid is used where the equa-
tion of balance of momentum is solved. This is in contrast to other meshfree/particle methods
where the momentum equations are solved on the particles. We refer to Fig 2 for a graphical
illustration of the discretisation in MPM. The grid covers the solid in its initial configuration
as well as the region in which the solid is expected to move. The particles interact with other
particles in the same body, with other solid bodies, or with fluids through a background Eu-
lerian grid. In general, after each time step the grid is discarded i.e., it is reset to the initial
configuration. This is the key difference between MPM and the updated Lagrange FEM, cf.
Fig. 3 and makes the method distortion-free. The MPM for solid mechanics conserves mass
and momentum by construction, but energy conservation is not explicitly enforced. Advantages
of the MPM algorithm include the absence of mesh-tangling problems and error-free advection
of material properties via the motion of the material points. Yet another advantage is that fully
coupled, fluid-structure interaction problems can be handled in a straightforward manner. A no
slip contact algorithm is automatic to the method. That is, it comes at no additional computa-
tional expense.
Applications of the MPM are emerging in ice-sheet models for climate simulation
Sulsky et al. [2007] as well as more traditional explosive-related simulations ?. In 2004, nu-
merical noise caused by cell-boundary crossing of material points was identified, and the gener-
alized interpolation material point (GIMP) method was introduced to avoid the noise ?. When
1
Advanced MPM formulations do track the subregion shape.

3
 2 INTRODUCTION

Figure 2: Material point method: Lagrangian material points overlaid over an Eulerian grid.
The dotted lines denote the physical domains for each particle. These domains can be found by
Voronoi tesselation.

Figure 3: Material point method: Lagrangian material points overlaid over an Eulerian grid.
Note that the material points can move to other elements. This is contrast to integration points
in FEM which always stay at the same element.

4
4 WEAK FORM AND DISCRETISATION

the choice of characteristic function is the Dirac delta function, the traditional MPM formula-
tion is recovered exactly. fluid-membrane interaction York et al. [1999, 2000]. ? comparative
study of SPH and MPM for hypervelocity impact problems. The conclusion was that the mate-
rial point method is an efficient and promising method for simulating the hypervelocity impact
problems. geo-engineering landslide.

3 Governing equations

4 Weak form and discretisation

The MPM also uses the weak formulation as in FEM which is given as

∂δui s
Z Z Z Z
ρδui ai dΩ + ρ σ dΩ = ρδui bi dΩ + ρδui t̄i dΓ (1)
Ω Ω ∂xj ij Ω Γt

where Ω denotes the current configuration, σijs is the specific stresses i.e., σijs = σij /ρ.
The whole material domain is described with a set of material points, and it is assumed that
the whole mass of a material subdomain is concentrated at the corresponding material point,
which means that the mass density field is expressed as
np
X
ρ(x, t) = Mp δ(x − xp ) (2)
p=1

where δ is the Dirac delta function with dimension of the inverse of volume. Substitution of
Equation (2) into Equation (1) results in

np np  np np
X X ∂δui  s
X X
Mp δui (xp )ai (xp )+ Mp σ (xp ) = Mp δui (xp )bi (xp )+ Mp δui (xp )t̄i (xp )
p=1 p=1
∂xj (xp ) ij p=1 p=1
(3)
At this point, MPM is identical to FEM in which the material points are served as integration
points. The background Eulerian grid serves as a finite element mesh with shape functions
{NI }nI=1
n
where nn denotes the total number of nodes. Thanks to the use of a grid, evaluation of
shape functions and derivatives are standard and does not involve neighbor search as in meshfree
methods such as SPH. The position and displacement of particle p are then given by
nn
X
xi (xp ) = NI (xp )xiI (4)
I=1

nn
X
ui (xp ) = NI (xp )uiI (5)
I=1

5
 4 WEAK FORM AND DISCRETISATION

where xiI is the i component of the position vector of node I. In 3D, one writes xI =
(x1I , x2I , x3I ). And uiI is the i component of the displacement vector of node I. Subscript
I denotes the value of grid node I, and subscript p denotes the value of particle p.
The velocity and acceleration fields are given by
nn
X
vi (xp ) = NI (xp )viI (6)
I=1
nn
X
ai (xp ) = NI (xp )aiI (7)
I=1

where viI , aiI are the i component of the velocity and acceleration vectors of node I.
Using the Bubnov-Galerkin method, the virtual displacement field is approximated as
nn
X
δui (xp ) = NI (xp )δuiI (8)
I=1

thus  nn 
∂δui  X ∂NI 
= δuiI (9)
∂xj (xp ) I=1 ∂xj (xp )
Substituting Equations (8), (7) and (9) into Equation (3) leads to

np np
"n #" nn
# "n #
n n

X X X X X ∂NI 
Mp NI (xp )δuiI NJ (xp )aiJ − Mp δuiI σijs (xp ) =
p=1 I=1 J=1 p=1 I=1
∂x 
j (xp )
np n
"n # p
" nn
#
X X n X X
Mp NI (xp )δuiI bi (xp ) + Mp NI (xp )δuiI t̄i (xp ) (10)
p=1 I=1 p=1 I=1

Since δuiI are arbitrary, we obtain the following set of equations (i equations for each node
I, I = 1, . . . , nn )

np nn np 
X X X ∂NI 
Mp NI (xp )( NJ (xp )aiJ ) − Mp σijs (xp ) =
p=1 J=1 p=1
∂x 
j (xp )
np np
X X
Mp NI (xp )bi (xp ) + Mp NI (xp )t̄i (xp ) (11)
p=1 p=1

which can be written in the following compact form

mIJ aJ = fIext + fIint (12)

where mIJ , fIext , fIint are the consistent mass matrix, the external force vector and the internal
force vector. This equation is exactly identical to FEM.

6
5 EXPLICIT TIME INTEGRATION

The consistent mass matrix is given by

np
X
mIJ = Mp NI (xp )NJ (xp ) (13)
p=1

Note that the mass matrix is not constant as in FEM but changes in time because the material
points move while the grid nodes are fixed. Therefore, for every time step, an inversion of the
mass matrix must be carried out. Thus, explicit integration is usually employed in MPM with
the use of a diagonal lumped mass matrix. The diagonal mass matrix is the standard row sum
matrix in which the diagonal terms are given by
nn np
nn X np
X X X
mI = mIJ = Mp NI (xp )NJ (xp ) = Mp NI (xp ) (14)
J=1 J=1 p=1 p=1

where Equation (13) and the partition of unity property of FE shape functions were used.
The external force vector is written as
np np
X X
fIext = Mp NI (xp )b(xp ) + Mp NI (xp )t̄(xp ) (15)
p=1 p=1

and the internal force vector as

np np
X X
fIint =− Mp /ρp σp ∇NI (xp ) = − Vp σ∇NI (xp ) (16)
p=1 p=1

where ∇NI = ( ∂N I ∂NI ∂NI T

, ,
∂x1 ∂x2 ∂x3
) denotes the gradient of the shape function; Vp is the volume
of particle p. Note that particle densities are defined as the ratio of particle mass to particle
volume. Note also that σp is a 3 × 3 matrix in 3D. This is different from the usual stress vector
written in Voigt notation in FEM.

5 Explicit time integration

If a lumped mass matrix is used, Equation (12) becomes

mtI atI = fIext + fIint ≡ fIt (17)

and after being multiplied with ∆t, one gets

mtI (vIt+∆t − vIt ) = fIt ∆t (18)

which permits the computation of the updated nodal momenta (mv)t+∆t I . Note that essential
boundary conditions must be applied before (18). It was shown in Ref. [39] that the explicit
procedure of time integration of the dynamic equations required shorter time increment when
the MPM is utilized than in the case of using the finite element method. However, before

7
 5 EXPLICIT TIME INTEGRATION

node particle

(a) node to particle (b) particle to node

Figure 4: FE mapping in the material point method: (a) nodes to particles mapping and (b)
particles to nodes mapping. The former is the standard FE mapping (from nodal displacements
to integration points’ strains). The latter is something unpopular in FEM although this is also
used for FE visualization (mapping the stresses of integration points to the nodes). Note also
that the former mapping is local and the latter is non-local (for C 0 shape functions). It should
be emphasized that some high order Lagrange basis functions are negative and therefore they
are not used in the particles to nodes mapping. High order B-splines are always positive and
hence can be used.

solving Equation (17), one needs the nodal masses and momenta which are not stored at the
nodes. Therefore, they are mapped from the particles using the shape functions

X
mtI = NI (xtp )Mp
p
X (19)
(mv)tI = NI (xtp )(Mv)tp
p

Note that this mapping, graphically illustrated in Fig. (4), is a kind of extrapolation and is a
non-local operation.
The particle velocities and positions are updated as follows

X
vpt+∆t = vpt + ∆t NI (xp )fIt /mtI (20)
I
X
xt+∆t
p = xtp + ∆t NI (xp )(mv)t+∆t
I /mtI (21)
I

which are standard FE interpolations, cf. Fig. (4). Note that this is the momentum formulation
in which momentum is used instead of velocity as much as possible, thus avoiding divisions
by potentially small nodal masses. Because the velocity field is single-valued, interpenetration
of material is precluded. It should be emphasized that this is totally an updated Lagrangian
formulation since the grid nodes are moved to new positions xt+∆t
I = vIt + ∆tvIt+∆t . However

8
5 EXPLICIT TIME INTEGRATION

node

particle
Figure 5: Troubled nodes with nearly zero mass resulting in infinite acceleration (node 2).

the grid nodes are not explicitly moved because it will be discarded anyway at the end of the
time step.
Next, one needs to compute the updated particle gradient velocities using the nodal veloc-
ities vIt+∆t . However if the velocities were computed as vIt+∆t = (mv)t+∆t I /mtI from Equa-
t
tion (18) the velocities would then be infinite if the mass mI is small. This happens when a
particle is very closed to a node having only one particle within its support, Fig. 5. There are at
least three ways to treat this issue: (i) using a cutoff value to detect small nodal masses, (ii) a
modified USL and (iii) USF. In the first approach, the nodal velocities are computed as
t+∆t
 (mv)I

if mtI > tol
vIt+∆t = mtI (22)
0 otherwise


which requires an extra parameter (a cutoff value) in the algorithm which is not clear how to
choose. Even a good cutoff value can be chosen, it produces an undesirable constraint which
should not be in the system. For example it results in non-zero stresses (albeit small) in particles
which are moving with the same speed.
The second way Sulsky et al. [1994] is to map the particle velocities back to the nodal
velocities using
X
(mv)t+∆t
I = NI (xp )(Mv)t+∆t
p (23)
p

and thus

NI (xp )(Mv)t+∆t NI (xp )(Mv)t+∆t

P P
(mv)t+∆t p p p p
vIt+∆t = I
= = (24)
mtI t
P
p NI (xp )Mp mI

The appearance of the shape functions in both numerator and denominator, in the second equal-
ity, cancel out its role and the numerical problem is thus cured. Our implementation uses the
final equality because mtI was computed in Equation (19). Note, however, that this option is
inefficient for Equation (23) is a particle-to-node map within the node-to-particle phase.
Next particle velocity gradients are computed

9
 6 IMPLEMENTATION

X
Lt+∆t
p ≡ ∇vpt+∆t = ∇NI (xp )vIt+∆t (25)
I
where Lp is a 3 × 3 matrix of which components are Lij = vi,j (in three dimensions). It is
a 2 × 2 matrix in two dimensions. From the velocity gradients one can compute the gradient
deformation using the relation Ḟ = LF. Utilizing a forward Euler method, one can write

Ft+∆t − Ft
= Lt+∆t Ft (26)
∆t
And the gradient deformation tensor is thus computed

Ft+∆t
p = (I + Lt+∆t
p ∆t)Ftp (27)
which allows to compute the updated particle volume Vpt+∆t= det Ft+∆t
p Vp0 .In the above
equation, I is the identity matrix.
Next, the particle stresses are updated. This depends on the constitutive model. For example,
one might need to compute the strain increment
∆ep = symLt+∆t


p ∆t (28)
and using it to compute the stress increment ∆σp . The updated particle stresses are given by

σpt+∆t = σpt + ∆σp (29)

In Box 5.1 the step-by-step algorithm of an explicit MPM is given. In the literature it is
referred to as USL (Update Stress Last). In Bardenhagen [2002], another MPM implementation
named USF (Update Stress First) was presented and for completeness, the USF procedure is
p there is an restriction on the time step dt < dx/c,
given in Box 5.2. For explicit dynamics,
where dx is the grid spacing and c = E/ρ is the speed of sound in the material.

6 Implementation
In this section, we present a serial implementation of MPM that reuses existing FE resources.
Before doing so, it is beneficial to recognize the resemblances and differences between MPM
and FEM:
• MPM can be seen as the updated Lagrangian finite element method in which material
points serve as integration points and uncoupled from the Eulerian grid. Practically ma-
terial points, where constitutive equations are applied, move from elements to elements
rather than remain at the Gauss points of an element. This requires to track the movement
of the material points.
• Contrary to standard FEM where the mesh is fitted to the physical domain of interest,
the Eulerian grid in MPM is dimensioned sufficiently large to cover the deformed solid.
Therefore, at a certain time step, there are elements with no material points. They are
labeled as inactive elements and skipped in the assembly process.

10
6 IMPLEMENTATION

Box 5.1 Solution procedure of an explicit MPM code: USL formulation.

1. Initialisation phase

(a) Particle distribution in the undeformed configuration

(b) Grid set up
(c) Initialise particle quantities such as mass, stress, strain etc.

2. Solution phase for time step t to t + ∆t

(a) Mapping from particles to nodes

i. Compute nodal mass mtI = p NI (xtp )Mp
P

ii. Compute nodal momentum (mv)tI = p NI (xtp )(Mv)tp

P

iii. Compute external force fIext,t

iv. Compute internal force fIint = − np=1
P p
Vp σp ∇NI (xp )
v. Compute nodal force fI = fIext + fIint
(b) Update the momenta (mv)t+∆t
I = (mv)tI + fI ∆t
(c) Mapping from nodes to particles
i. Update particle velocities vpt+∆t = vpt + ∆t I NI (xp )fIt /mtI
P

ii. Update particle positions xt+∆t = xtp + ∆t I NI (xp )(mv)t+∆t /mtI

P
p I
iii. Update nodal velocities vIt+∆t = (mv)It+∆t /mtI
Update gradient velocity Lt+∆t = I ∇NI (xp )vIt+∆t
P
iv. p
v. Updated gradient deformation tensor Fpt+∆t = (I + Lt+∆t
p ∆t)Ftp
vi. Update volume Vpt+∆t = det Ft+∆t
p Vp0 .
vii. Update stresses σpt+∆t = σpt + ∆σp

3. Reset the grid (if it was updated) and advance to the next time step.

11
 6 IMPLEMENTATION

Box 5.2 Solution procedure of an explicit MPM code: USF formulation.

1. Initialisation phase

2. Solution phase for time step t to t + ∆t

(a) Mapping from particles to nodes

i. Compute nodal mass mtI = p NI (xtp )Mp
P

ii. Compute nodal momentum (mv)tI = p NI (xtp )(Mv)tp

P

iii. Compute nodal velocities vIt = (mv)tI /mtI

iv. Compute gradient velocity Ltp = I ∇NI (xp )vIt
P

v. Compute gradient deformation tensor Ftp = (I + Ltp ∆t)Ftp

vi. Update volume Vpt = det Ft+∆t
p Vp0
vii. Update stresses σpt = σpt + ∆σp
viii. Compute external force fIext,t
ix. Compute internal force fIint = − np=1
P p
Vp σp ∇NI (xp )
x. Compute nodal force fI = fIext + fIint
(b) Update the momenta (mv)t+∆t
I = (mv)tI + fI ∆t
(c) Mapping from nodes to particles
i. Update particle velocities vpt+∆t = vpt + ∆t I NI (xp )fIt /mtI
P

ii. Update particle positions xt+∆t = xtp + ∆t I NI (xp )(mv)t+∆t /mtI

P
p I

3. Reset the grid (if it was updated) and advance to the next time step.

12
6 IMPLEMENTATION

• Unlike FEM, in MPM there is a mapping from Gauss points (particles) to nodes.

Some advantages of MPM are listed as follows

• Mesh distortion of Lagrangian FEM is eliminated.

• The problem of free surface is easy to solve.

• The problem of no-slip self-contact is solved automatically in the case of granular ma-
terial: there is no penetration of particles because the velocities of the material point are
single valued.

• Boundary conditions can be applied as easily as in the case of standard FEM.

• Adding material points during calculations is relatively easy as there is no coupling be-
tween them and the grid.

• The MPM can be implemented in a FEM program in a relatively easy way in constrast to
meshfree methods or ALE methods.

• Solving fluid-structure interaction is relatively straightforward: the background grid is

also used as an Eulerian grid for the fluid.

The MPM suffers from the following disadvantages

• Efficiency of MPM is lower than that of FEM due to the mappings between the back-
ground grid and the particles and the accuracy of particles quadrature used in MPM is
lower than that of Gauss quadrature used in FEM.
• This method, in analogy with meshless methods, also introduces new visualization chal-
lenges.
• Whereas the finite element method (FEM) has long-established basic verification stan-
dards (patch tests, convergence testing, etc.), no such standards have been universally
adopted within the particle method community.

Available MPM codes

• Uintah, http://www.uintah.utah.edu
• MPM-GIMP, http://sourceforge.net/p/mpmgimp/home/Home/

• Mechsys, http://mechsys.nongnu.org/index.shtml

• NairnMPM, http://osupdocs.forestry.oregonstate.edu/index.php/Main_Page

• CartaBlanca

All are written in C++ except the last one which is written in Java.

13
6.1 Eulerian grid 6 IMPLEMENTATION

Figure 6: A two dimensional structured grid.

6.1 Eulerian grid

In MPM a structured Eulerian grid is usually used for its advantages. Among other bene-
fits, a uniform Cartesian grid eliminates the need for computationally expensive neighborhood
searches during particle-mesh interaction. A structured mesh is illustrated in Fig. 6 for 2D
cases. For a particle p with coordinates (xp , yp ), it is straightforward to track which element it
belongs to using the following equation

e = [floor((xp − xmin )/∆x ) + 1] + nelx [floor((yp − ymin )/∆y)] (30)

where nelx denotes the number of elements along the x direction and (xmin , ymin ) are the mini-
mum node coordinates and ∆x, ∆y are the nodal spacing in the x and y directions, respectively.
In order to keep the implementation similar to FEM as much as possible, at every time step,
the elements must know which particles they are storing. Listing 1 gives a simple data structure
implemented in Matlab for this purpose.
Listing 1: Matlab data structures for mesh-particle interaction
pElems = o n e s ( pCount , 1 ) ;
m p o i n t s = c e l l ( elem Co u n t , 1 ) ;

f o r p = 1 : p Co u n t
x = xp ( p , 1 ) ;
y = xp ( p , 2 ) ;
e = f l o o r ( x / d e l t a x ) + 1 + numx2∗ f l o o r ( y / d e l t a y ) ;
pElems ( p ) = e ; % p a r t i c l e " p " s t a y s i n e l e m e n t " e "
end

f o r e = 1 : elem Co u n t
i d = f i n d ( pElems==e ) ;
m p o i n t s { e }= i d ; % m p o i n t s { e}−> i n d i c e s o f p a r t i c l e s i n " e "
end

14
6 IMPLEMENTATION 6.2 Shape functions

Figure 7: Material point method: One dimensional shape functions defined in global coordinate
system.

6.2 Shape functions

Although any grid can be used in MPM, a Cartesian grid is usually chosen for computational
convenience reasons. In order to avoid finding the natural coordinates of material points if shape
functions are defined in the parameter space, in MPM shape functions are conveniently defined
in the global coordinate system. In 1D, the shape functions are defined as
(
x 1 − |x − xI |/lx if |x − xI | ≤ lx
NI (x) = (31)
0 else
where lx denotes the nodal spacing or element size in the x direction. For 2D, the shape func-
tions are simply the tensor-product of the two shape functions along the x and y directions

NI (x, y) = NIx (x)N y (y) (32)

Remark 6.1. If shape functions are defined in the parent domain as most isoparametric elements
are i.e., NI = NI (ξ, η), then one has to perform one extra step–after getting the updated particle
positions, determine the natural coordinates (ξ, η) of them. For a 2D structured grid with linear
elements, it is straightforward to do so as

2x − (xe1 + xe2 )
ξ=
xe2 − xe1
(33)
2y − (y1e + y2e )
η=
y2e − y1e

where (xe1 , y1e ) and (xe2 , y2e ) denote the lower-left and upper-right nodes of element e, respec-
tively. Extension to 3D is straightforward.

15
6.3 Initial particle distribution 6 IMPLEMENTATION

(a) parametric mesh (b) physical mesh

Figure 8: A two dimensional structured grid with B-splines.

6.2.1 B-splines basis functions

Steffen et al. [2008] showed that for simple problems, the use of cubic splines improves the
spatial convergence properties of method as grid-crossing errors are reduced.
Given a knot vector Ξ1 = {ξ1 , ξ2 , . . . , ξn+p+1}, the associated set of B-spline basis functions
{Ni,p }ni=1 are defined recursively by the Cox-de-Boor formula, starting with the zeroth order
basis function (p = 0) (
1 if ξi ≤ ξ < ξi+1 ,
Ni,0 (ξ) = (34)
0 otherwise,
and for a polynomial order p ≥ 1
ξ − ξi ξi+p+1 − ξ
Ni,p (ξ) = Ni,p−1 (ξ) + Ni+1,p−1 (ξ). (35)
ξi+p − ξi ξi+p+1 − ξi+1
in which fractions of the form 0/0 are defined as zero.
In two dimensions, the knot vectors are Ξ1 = {ξ1, ξ2 , . . . , ξn+p+1} and Ξ2 =
{η1 , η2 , . . . , ηm+q+1 } and the bi-variate B-spline basis functions are defined as

Ni,j (ξ, η) = Ni,p (ξ)Mj,q (η) (36)

The mapping between the parameter space and the physical space is always linear regardless
of the basis order and is given by

x = (xmax − xmin )ξ + xmin , y = (ymax − ymin )η + ymin (37)

p−1 0
High order B-spline basis functions are C not C , the connectivity of elements is different
from standard finite elements. Fig. 9 illustrates this difference in case of cubic B-splines. Apart
from this B-splines MPM is exactly identical to standard MPM.

6.3 Initial particle distribution

There are numerous ways to obtain the initial particle distribution depends on the geometry of
the object and/or the available tools. For example, in order to generate particles for a circular

16
6 IMPLEMENTATION 6.3 Initial particle distribution

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.2 0.4 0.6 0.8 1
particle elements

1 2 3 4 5 6 7 8 nodes
Figure 9: One dimensional cubic (p = 3) B-spline basis functions on a open uniform knot Ξ =
{0, 0, 0, 0, 0.2, 0.4, 0.6, 0.8, 1, 1, 1, 1}. There are 8 nodes (control points in CAD terminology)
and 5 elements (or knot spans in CAD). At can be seen from the figure, at any point there are 4
(= (p + 1)) non-zero basis functions. Therefore each element has 4 nodes. The firstelement’s
connectivity is [1, 2, 3, 4].

17
6.4 Visualisation 7 EXAMPLES

Figure 10: Initial particle distribution: obtained using the available FE mesh generators.

disk, one can use a mesh generator to partition the disk into a set of triangles. The particles
can be taken as the centers of these triangles, Fig. 10. Alternatively, integration points of the
triangular elements are mapped to the global coordinate systems and then are used as material
points. The area of each triangle can be easily obtained and thus the particle volumes and
masses can be determined.

6.4 Visualisation
In FEM visualisation is typically performed on the mesh. Information stored at the integration
points are extrapolated to the nodes to this end. In MPM, the computational grid is fixed and
thus not suitable for visualisation. There are at least two approaches to visualizing MPM results
• Particle visualization: the particles are visualized directly as spheres. Particle data (posi-
tion, stresses etc.) are written to files (say VTK files) and can be processed by Paraview
or Visit.
• Particle mesh visualization: in standard MPM, one can generate a mesh from the parti-
cle positions using a Delaunay triangulation and use this mesh for visualization. Note
however that this method is expensive since the particles positions are evolving in time.

7 Examples
In this section various examples in one and two dimensions are provided to verify the imple-
mentation as well as to demonstrate the performance of the MPM.

7.1 1D examples
As the simplest MPM example, let us consider the vibration of a single material point as shown
in Fig. 11. The bar is represented by a single point initially located at Xp = L/2, which has an
initial velocity v0 . The material is linear elastic.

18
7 EXAMPLES 7.1 1D examples

Figure 11: Vibration of a single material point (solid point).

The exact solution is given by

1p
v(t) = v0 cos(ωt), ω= E/ρ (38)
L
for the velocity and
hv i
0
x(t) = x0 exp sin(ωt) (39)
Lw
for the position. The density ρ is constant and equals one. The constitutive equation is σ̇ = E ǫ̇,
ǫ̇ = dv/dx where E is the Young’s p modulus. The grid consists of one two-noded element.
The elastic wave speed is c = E/ρ = 2π. Boundary conditions are imposed on the grid
and demand that both the node velocity and the acceleration at x = 0 (the left node) be zero
throughout the simulation. The Matlab implementation is given in Listing 2. Comparison
between MPM and the exact solutions are given in Fig. 12. Kinetic, strain and total energies
are plotted in Fig. 13 from which one observes that energy is conserved. The strain and kinetic
energy are computed as
1 1
U = σp ǫp Vp , K = vp2 Mp (40)
2 2

Listing 2: Matlab implementation

% Computational grid
nodes = [0 L ] ;
elements = [1 2 ] ;
elem Co u n t = s i z e ( e l e m e n t s , 1 ) ;
n o d eCo u n t = s i z e ( n o d es , 2 ) ;
% Material points
xp = 0 . 5 ∗ L ; % position
Mp = 1 ; % mass
Vp = 1 ; % volume
vp = 0 . 1 ; % velocity
s = 0.; % stress
q = Mp∗ vp ; % momentum
% Time l o o p
time = 0 . 1 ;
dtime = 0 . 0 1 ;
t = 0;

19
7.1 1D examples 7 EXAMPLES

t a = [ ] ; va = [ ] ; xa = [ ] ;
while ( t < time )
% s h a p e f u n c t i o n s and d e r i v a t i v e s
N1 = 1 − abs ( xp−n o d e s ( 1 ) ) / L ;
N2 = 1 − abs ( xp−n o d e s ( 2 ) ) / L ;
dN1 = −1/L ;
dN2 = 1 / L ;
% p a r t i c l e mass and momentum t o n o d e
m1 = N1∗Mp;
m2 = N2∗Mp;
mv1 = N1∗ q ;
mv2 = N2∗ q ;

mv1 = 0 ; % Bo u n d ar y c o n d i t i o n v=0 a t n o d e 1
% internal force
f i n t 1 = − Vp∗ s ∗dN1 ;
f i n t 2 = − Vp∗ s ∗dN2 ;

f1 = fint1 ; % t h e r e i s no e x t e r n a l f o r c e
f2 = fint2 ;
% update n o d a l momenta
f1 = 0; % Bo u n d ar y c o n d i t i o n s f 1 = m1∗ a1 , a1 =0

mv1 = mv1 + f 1 ∗ d t i m e ;
mv2 = mv2 + f 2 ∗ d t i m e ;

% u p d a t e p a r t i c l e v e l o c i t y and p o s i t i o n
vp = vp + d t i m e ∗ ( N1∗ f 1 / m1 + N2∗ f 2 / m2 ) ;
xp = xp + d t i m e ∗ ( N1∗mv1 / m1 + N2∗mv2 / m2 ) ;

q = Mp∗ vp ; % momentum

v1 = N1∗Mp∗ vp / m1 ;
v2 = N2∗Mp∗ vp / m2 ;
v1 = 0 ; % boundary c o n d i t i o n
Lp = dN1 ∗ v1 + dN2 ∗ v2 ; % gradient velocity
dEps = d t i m e ∗ Lp ; % s t r a i n increment
s = s + E ∗ dEps ; % s t r e s s update
% s t o r e time , v e l o c t y f o r p l o t t i n g
ta = [ ta ; t ];
va = [ va ; vp ] ;
xa = [ xa ; xp ] ;
% advance to the next time st e p
t = t + dtime ;
end

20
7 EXAMPLES 7.1 1D examples

0.15
MPM
0.1 Exact

Velocity 0.05

−0.05

−0.1

−0.15

−0.2
0 2 4 6 8 10
Time

10
MPM
8 Exact
6
4
Displacement

2
0
−2
−4
−6
−8
0 2 4 6 8 10
Time

Figure 12: Vibration of a single material point.

−3
x 10
7
kinetic
strain
6 total

4
Energy

0
0 0.5 1 1.5 2
Time

Figure 13: Vibration of a single material point: kinetic, strain and total energies.

21
7.2 Axial vibration of a continuum bar 7 EXAMPLES

7.2 Axial vibration of a continuum bar

In this example, we study the axial vibration of a continuum bar with Young’s modulus E = 100
and the bar length L = 25.
Exact solutions for mode 1 are

v(x, t) = v0 cos(ω1 t) sin(β1 x) (41)

v0
u(x, t) = sin(ω1 t) sin(β1 x) (42)
ω1
π
p π
where ω1 = 2L E/ρ and β1 = 2L .
The initial velocity is given by

v(x, 0) = v0 sin(β1 x) (43)

The grid consists of 13 two-noded line elements (14 grid nodes) and 13 material points
placed at the center of the elements are used. The initialisation step implemented in Matlab is
given in Listing 3. The solution step implemented in Matlab is given in Listing 4. The results
are given in Fig. 14 and the evolution of kinetic energy, strain energy and their sum (the total
energy) is plotted in Fig. 15. The time increment is chosen as ∆t = 0.1∆x/c where ∆x denotes
the nodal spacing.
Listing 3: Vibration of a continuum bar: initialisation.
% C o m p u t a t i o n a l g r i d : two−n o d ed e l e m e n t s
elem Co u n t = 1 3 ;
nodes = l i n s p a c e ( 0 , L , elem Co u n t + 1 ) ;
e l e m e n t s = z e r o s ( elem Co u n t , 2 ) ;
f o r i e = 1 : elem Co u n t
elements ( ie , : ) = i e : i e +1;
end
n o d eCo u n t = s i z e ( n o d es , 2 ) ;

v e l o = z e r o s ( n o d eCo u n t , 1 ) ;
mId = f l o o r ( ( elem Co u n t ) / 2 ) + 1 ; % i d o f t h e c e n t e r mass p a r t i c l e
% M a t e r i a l p o i n t s : c e n t e r s o f t h e e l e m e n t s one p a r t i c l e p e r e l e m e n t
c = sqrt (E/ rho ) ;
beta1 = pi / 2 / L ;
omega1 = b e t a 1 ∗ c ;
d e l t a x = L / elem Co u n t ; % n o d a l s p a c i n g
% m at e ri a l p o i n t s at the c e n t e r of elements
xp = z e r o s ( elem Co u n t , 1 ) ;
f o r p = 1 : elem Co u n t −1
xp ( p ) = 0 . 5 ∗ ( n o d e s ( p ) + n o d e s ( p + 1 ) ) ;
end
p Co u n t = l e n g t h ( xp ) ;

Mp = d e l t a x ∗ o n e s ( pCount , 1 ) ; % mass
Vp = d e l t a x ∗ o n e s ( pCount , 1 ) ; % volume

22
7 EXAMPLES 7.2 Axial vibration of a continuum bar

Fp = o n e s ( pCount , 1 ) ; % gradient deformation

Vp0 = Vp ; % i n i t i a l volume
s p = z e r o s ( pCount , 1 ) ; % stress
vp = z e r o s ( pCount , 1 ) ; % velocity
% initial velocities
f o r p = 1 : p Co u n t
vp ( p ) = 0 . 1 ∗ s i n ( b e t a 1 ∗ xp ( p ) ) ;
end

Listing 4: Vibration of a continuum bar: processing step.

dtime = 0.1∗ d e l t a x / c ;
time = 100;
t = 0;

t a = [ ] ; va = [ ] ; xa = [ ] ;

nmass = z e r o s ( n o d eCo u n t ,1); % nodal mass v e c t o r

nmomentum = z e r o s ( n o d eCo u n t ,1); % nodal momentum v e c t o r
niforce = z e r o s ( n o d eCo u n t ,1); % nodal i nternal force vector
neforce = z e r o s ( n o d eCo u n t ,1); % nodal external force vector

while ( t < time )

nmass ( : ) = 0;
nmomentum ( : ) = 0 ;
niforce (:) = 0;
% loop over computational c e l l s or elements
f o r e = 1 : elem Co u n t
e s c t r = elements ( e , : ) ;
en o d e = n o d e s ( e s c t r ) ;
Le = en o d e (2) − en o d e ( 1 ) ;
m p ts = m p o i n t s { e } ;
% loop over p a r t i c l e s
f o r p = 1 : l e n g t h ( m p ts )
p i d = m p ts ( p ) ;
N1 = 1 − abs ( xp ( p i d )− en o d e ( 1 ) ) / Le ;
N2 = 1 − abs ( xp ( p i d )− en o d e ( 2 ) ) / Le ;
dN1 = −1/Le ;
dN2 = 1 / Le ;
% p a r t i c l e mass and momentum t o n o d e
nmass ( e s c t r ( 1 ) ) += N1∗Mp( p i d ) ;
nmass ( e s c t r ( 2 ) ) += N2∗Mp( p i d ) ;
nmomentum ( e s c t r ( 1 ) ) += N1∗Mp( p i d ) ∗ vp ( p i d ) ;
nmomentum ( e s c t r ( 2 ) ) += N2∗Mp( p i d ) ∗ vp ( p i d ) ;
% internal force
n i f o r c e ( e s c t r ( 1 ) ) = n i f o r c e ( e s c t r ( 1 ) ) − Vp ( p i d ) ∗ s p ( p i d ) ∗ dN1 ;
n i f o r c e ( e s c t r ( 2 ) ) = n i f o r c e ( e s c t r ( 2 ) ) − Vp ( p i d ) ∗ s p ( p i d ) ∗ dN2 ;
end
end
% u p d a t e n o d a l momenta
nmomentum ( 1 ) = 0 ; % Bo u n d ar y c o n d i t i o n s f 1 = m1∗ a1 , a1 =0
niforce (1) = 0;

23
7.3 Impact of two elastic bodies 7 EXAMPLES

f o r i = 1 : n o d eCo u n t
nmomentum ( i ) = nmomentum ( i ) + n i f o r c e ( i ) ∗ d t i m e ;
end
% u p d a t e p a r t i c l e v e l o c i t y and p o s i t i o n and s t r e s s e s
f o r e = 1 : elem Co u n t
e s c t r = elements ( e , : ) ;
en o d e = n o d e s ( e s c t r ) ;
Le = en o d e (2) − en o d e ( 1 ) ;
m p ts = m p o i n t s { e } ;
% loop over p a r t i c l e s
f o r p = 1 : l e n g t h ( m p ts )
p i d = m p ts ( p ) ;
N1 = 1 − abs ( xp ( p i d )− en o d e ( 1 ) ) / Le ;
N2 = 1 − abs ( xp ( p i d )− en o d e ( 2 ) ) / Le ;
dN1 = −1/Le ; dN2 = 1 / Le ;
i f nmass ( e s c t r ( 1 ) ) > t o l
vp ( p i d ) += d t i m e ∗ N1∗ n i f o r c e ( e s c t r ( 1 ) ) / nmass ( e s c t r ( 1 ) ) ;
end
i f nmass ( e s c t r ( 2 ) ) > t o l
vp ( p i d ) += d t i m e ∗ N2∗ n i f o r c e ( e s c t r ( 2 ) ) / nmass ( e s c t r ( 2 ) ) ;
end
xp ( p i d ) += d t i m e ∗ ( N1∗nmomentum ( e s c t r ( 1 ) ) / nmass ( e s c t r ( 1 ) ) +
N2∗nmomentum ( e s c t r ( 2 ) ) / nmass ( e s c t r ( 2 ) ) ) ;
v1 = nmomentum ( e s c t r ( 1 ) ) / nmass ( e s c t r ( 1 ) ) ;
v2 = nmomentum ( e s c t r ( 2 ) ) / nmass ( e s c t r ( 2 ) ) ;
%i f ( e s c t r ( 1 ) == 1 ) v1 = 0 ; end
Lp = dN1 ∗ v1 + dN2 ∗ v2 ; % gradient velocity
Fp ( p i d ) = ( 1 + Lp∗ d t i m e ) ∗ Fp ( p i d ) ; % gradient deformation
Vp ( p i d ) = Fp ( p i d ) ∗ Vp0 ( p i d ) ; % volume
dEps = d t i m e ∗ Lp ; % s t r a i n increment
s p ( p i d ) = s p ( p i d ) + E ∗ dEps ; % s t r e s s update
end
end
% s t o r e time , v e l o c t y f o r p l o t t i n g
t a = [ t a ; t ] ; va = [ va ; vp ( mId ) ] ; xa = [ xa ; xp ( mId ) ] ;
% advance to the next time st e p
t = t + dtime ;
end

7.3 Impact of two elastic bodies

Fig. 16 shows the problem of the impact of two identical elastic disks which move in opposite
direction towards each other Sulsky et al. [1994]. The computational domain is a square which
is discretised into 20 × 20 Q4 elements. A plane strain condition is assumed. The mesh and
the initial particle distribution are shown in Fig. 17. There are 320 particles for two disks which
are obtained by meshing (using Gmsh) the two disks and take the centers as the initial particle
positions. Initial condition for this problem is the initial velocities of the particles, vp = v for
lower-left particles and vp = −v for upper-right particles. There is no boundary conditions
in this problem since the simulation stops before the particles after impact move out of the

24
7 EXAMPLES 7.3 Impact of two elastic bodies

0.1
MPM
0.08 Exact
0.06
0.04
0.02
Velocity

0
−0.02
−0.04
−0.06
−0.08
−0.1
0 20 40 60 80 100
Time

Figure 14: Vibration of a continuum bar: comparison of the MPM velocity solution and the
exact solution.

0.08
kinetic
0.07 strain
total
0.06

0.05
Energy

0.04

0.03

0.02

0.01

0
0 20 40 60 80 100
Time

Figure 15: Vibration of a continuum bar: kinetic, strain and total energies.

25
7.3 Impact of two elastic bodies 7 EXAMPLES

Figure 16: Impact of two elastic bodies: problem statement.

computational box.
In order to check the energy conservation, the strain and kinetic energy are computed for
each time step. They are defined as
np np
X 1X
U= u p Vp , K= vp · vp Mp (44)
p=1
2 p=1

where up denotes the strain energy density of particle p, up = 1/2σp,ij ǫp,ij or explicitly as
 
1 κ+1 2 2 2
up = (σp,xx + σp,yy ) − 2(σp,xx σp,yy − σp,xy ) (45)
4µ 4

with µ and κ are the shear modulus and the Kolosov constant, respectively.
List 5 and 6 gives the code for initialising the particle data and nodal data. Lists 7, 8 and
9 presents the codes for the solution phase. Note that some C++ notations (+=) were adopted
to save space. It should be emphasized that the code has not been optimized to make it easy
to understand. One can, for example, store the shape functions and its derivatives computed in
List 8 so that they can be reused in List 9.
Listing 5: Two dimensional MPM: particle initialisation.
p Co u n t = numelem ; % # of p a r t i c l e s
Mp = o n e s ( pCount , 1 ) ; % mass
Vp = o n e s ( pCount , 1 ) ; % volume
Fp = o n e s ( pCount , 4 ) ; % gradient deformation
s = z e r o s ( pCount , 3 ) ; % stress
e p s = z e r o s ( pCount , 3 ) ; % strain
vp = z e r o s ( pCount , 2 ) ; % velocity
xp = z e r o s ( pCount , 2 ) ; % position
% i n i t i a l i s e p a r t i c l e p o s i t i o n , mass , volume , v e l o c i t y
f o r e = 1 : numelem
c o o r d = n o d e1 ( e l e m e n t 1 ( e , : ) , : ) ;
a = det ( [ coord , [ 1 ; 1 ; 1 ] ] ) / 2 ;
Vp ( e ) = a ;

26
7 EXAMPLES 7.3 Impact of two elastic bodies

Figure 17: Impact of two elastic bodies: Eulerian mesh and initial particle distribution.

Mp( e ) = a ∗ r h o ;
xp ( e , : ) = mean ( c o o r d ) ;
% t h i s i s p a r t i c u l a r f o r t h e two i m p a c t d i s k s ex am p le
i f xp ( e , 1 ) < 0 . 5
vp ( e , : ) = [ v v ] ;
else
vp ( e , : ) = [−v −v ] ;
end
Fp ( e , : ) = [ 1 0 0 1 ] ;
end
Vp0 = Vp ;

Listing 6: Two dimensional MPM: nodal initialisation.

% build the structured grid ...
% initialise nodal data
nmass = z e r o s ( n o d eCo u n t ,1); % nodal mass v e c t o r
nmomentum = z e r o s ( n o d eCo u n t ,2); % nodal momentum v e c t o r
niforce = z e r o s ( n o d eCo u n t ,2); % nodal i nternal force vector
neforce = z e r o s ( n o d eCo u n t ,2); % nodal external force vector

Listing 7: Two dimensional MPM: solution phase (explicit).

while ( t < time )
% re set grid data
nmass ( : ) = 0;
nmomentum ( : ) = 0 ;
niforce (:) = 0;
% p a r t i c l e t o g r i d nodes

27
7.3 Impact of two elastic bodies 7 EXAMPLES

See L i s t 8
% u p d a t e n o d a l momenta
nmomentum = nmomentum + n i f o r c e ∗ d t i m e ;
% nodes t o p a r t i c l e s
See L i s t 9
% s t o r e time , v e l o c t y f o r p l o t t i n g
p o s { i s t e p } = xp ;
v e l { i s t e p } = vp ;
% update the element p a r t i c l e l i s t
f o r p = 1 : p Co u n t
x = xp ( p , 1 ) ;
y = xp ( p , 2 ) ;
e = f l o o r ( x / d e l t a X ) + 1 + numx2∗ f l o o r ( y / d e l t a Y ) ;
pElems ( p ) = e ;
end
f o r e = 1 : elem Co u n t
i d = f i n d ( pElems ==e ) ;
m p o i n t s { e }= i d ;
end
% VTK o u t p u t
i f ( mod ( i s t e p , i n t e r v a l ) == 0 )
xp = p o s { i s t e p } ;
v t k F i l e = s p r i n t f ( ’../results/%s%d’ , v t k F i l e N a m e , i s t e p ) ;
V T K P a r t i c l e s ( xp , v t k F i l e , s ) ;
end
% advance to the next time st e p
t = t + dtime ; i s t e p = i s t e p + 1;
end

Listing 8: Two dimensional MPM: particles to nodes.

f o r e = 1 : elem Co u n t % loop over elements
e s c t r = element ( e , : ) ; % element co n n ec ti v it y
en o d e = n o d e ( e s c t r , : ) ; % e le m e n t node c o o r d s
m p ts = m p o i n t s { e } ; % p a r t i c l e s in s id e element e
% loop over p a r t i c l e s
f o r p = 1 : l e n g t h ( m p ts )
p i d = m p ts ( p ) ; % p a r t i c l e ID
p t ( 1 ) = ( 2 ∗ xp ( p i d , 1 ) − ( en o d e ( 1 , 1 ) + en o d e ( 2 , 1 ) ) ) / d e l t a X ;
p t ( 2 ) = ( 2 ∗ xp ( p i d , 2 ) − ( en o d e ( 2 , 2 ) + en o d e ( 3 , 2 ) ) ) / d e l t a Y ;

[N, dNdxi ] = l a g r a n g e _ b a s i s ( ’Q4’ , p t ) ; % element shape f u n c t i o n s

J0 = enode ’ ∗ dNdxi ; % element Jacobian matrix
invJ0 = inv ( J0 ) ;
dNdx = dNdxi ∗ i n v J 0 ;
% p a r t i c l e mass and momentum t o n o d e
s t r e s s = s ( pid , : ) ;
for i =1: length ( e s c t r )
id = esctr ( i ); % n o d e ID
dNIdx = dNdx ( i , 1 ) ;
dNIdy = dNdx ( i , 2 ) ;
nmass ( i d ) += N( i ) ∗Mp( p i d ) ;

28
7 EXAMPLES 7.3 Impact of two elastic bodies

nmomentum ( i d , : ) += N( i ) ∗Mp( p i d ) ∗ vp ( p i d , : ) ;
n i f o r c e ( id , 1 ) −= Vp ( p i d ) ∗ ( s t r e s s ( 1 ) ∗ dNIdx + s t r e s s ( 3 ) ∗ dNIdy ) ;
n i f o r c e ( id , 2 ) −= Vp ( p i d ) ∗ ( s t r e s s ( 3 ) ∗ dNIdx + s t r e s s ( 2 ) ∗ dNIdy ) ;
end
end
end

Listing 9: Two dimensional MPM: nodes to particles.

f o r e = 1 : elem Co u n t % loop over elements
e s c t r = element ( e , : ) ;
en o d e = n o d e ( e s c t r , : ) ;
m p ts = m p o i n t s { e } ;
% loop over p a r t i c l e s
f o r p = 1 : l e n g t h ( m p ts )
p i d = m p ts ( p ) ;
p t ( 1 ) = ( 2 ∗ xp ( p i d , 1 ) − ( en o d e ( 1 , 1 ) + en o d e ( 2 , 1 ) ) ) / d e l t a X ;
p t ( 2 ) = ( 2 ∗ xp ( p i d , 2 ) − ( en o d e ( 2 , 2 ) + en o d e ( 3 , 2 ) ) ) / d e l t a Y ;
[N, dNdxi ] = l a g r a n g e _ b a s i s ( ’Q4’ , p t ) ;
J0 = enode ’ ∗ dNdxi ;
invJ0 = inv ( J0 ) ;
dNdx = dNdxi ∗ i n v J 0 ;
Lp = z e r o s ( 2 , 2 ) ;
for i =1: length ( e s c t r )
id = e sc tr ( i ) ; % n o d e ID
vI = [0 0 ] ;
i f nmass ( i d ) > t o l
vp ( p i d , : ) += d t i m e ∗ N( i ) ∗ n i f o r c e ( i d , : ) / nmass ( i d ) ;
xp ( p i d , : ) += d t i m e ∗ N( i ) ∗ nmomentum ( i d , : ) / nmass ( i d ) ;
vI = nmomentum ( i d , : ) / nmass ( i d ) ;% n o d a l v e l o c i t y
end
Lp = Lp + v I ’ ∗ dNdx ( i , : ) ; % p a r t i c l e g r a d i e n t v e l o c i t y
end

F = ( [ 1 0 ; 0 1 ] + Lp∗ d t i m e ) ∗ r e s h a p e ( Fp ( p i d , : ) , 2 , 2 ) ;
Fp ( p i d , : ) = r e s h a p e ( F , 1 , 4 ) ;
Vp ( p i d ) = d e t ( F ) ∗ Vp0 ( p i d ) ;
dEps = d t i m e ∗ 0 . 5 ∗ ( Lp+Lp ’ ) ;
d sig m a = C ∗ [ dEps ( 1 , 1 ) ; dEps ( 2 , 2 ) ; 2 ∗ dEps ( 1 , 2 ) ] ;
s ( p i d , : ) = s ( p i d , : ) + dsigma ’ ;
e p s ( p i d , : ) = e p s ( p i d , : ) + [ dEps ( 1 , 1 ) dEps ( 2 , 2 ) 2∗ dEps ( 1 , 2 ) ] ;
end
end

The movement of two disks is given in Fig. 18. The collision occurs in a physically realistic
fashion, although no contact law has been specified. Fig. 19 plots the evolution of the kinetic,
strain and total energy.

29
7.3 Impact of two elastic bodies 7 EXAMPLES

time: 0.930000 time: 1.205000

time: 1.930000 time: 2.480000

Figure 18: Impact of two elastic bodies: time snapshots of the bodies. Top figures: twp bodies
move towards each other and collide. Bottom figures: they bounce back and move far from
each other. These figures were created in Matlab using the scatter command.

30
7 EXAMPLES 7.4 Impact of two elastic rings

2.5

2
Energy

1.5
kinetic
strain
1
total

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5
Time

Figure 19: Impact of two elastic bodies: evolution of strain, kinetic and total energies in time.

Figure 20: Impact of two elastic bodies: problem description. Dimensions are in millimeters.

7.4 Impact of two elastic rings

This example problem involves two hollow elastic cylinders, under the assumption of plane
strain, impacting each other Fig. 20. The material is a compressible Neo-Hookean with bulk
modulus K = 121.7 MPa, shear modulus µ = 26.1 MPa, the density is ρ = 1010 × 10−12
kg/mm3. The Cauchy stresses are computed as follows

1
σ= [µ0 (b − I) + λ0 ln JI] (46)
J
where J = det F and b = FFT is the left Cauchy strain tensor; λ0 = K − 2/3µ0 is the first
Lamé constant.
A grid of 100 × 60 elements is adopted
p and each cylinder is discretized by 5488 particles.
The time step was chosen as ∆t = 0.2 E/ρ = 1.5 × 10−6 s. For this problem, it is very

31
 8 GENERALIZED INTERPOLATION MATERIAL POINT METHOD-GIMP

hard to find a cutoff value in Equation (22). We therefore used Equation (24) to compute the
nodal velocities which are subsequently used to compute the particle velocity gradient Lp . For
visualization, the particle data (positions and stresses) are written to a VTP file (PolyData) and
processed in Paraview using the filter Glyph.

8 Generalized Interpolation Material Point Method-GIMP

The original MPM with C 0 shape functions suffers from the so-called cell crossing instability.
In order to illustrate this phenomenon, consider a 1D MPM discretisation shown in Fig. 22 in
which all particles have the same stress (i.e., uniform stress state), and the same volume and a
uniform element size. Each element has two particles, Fig. 22a. The internal force at node 2 is
identically zero in the absence of body force. When a particle have just moved to a new cell,
Fig. 22b, the internal force at node 2 is non-zero resulting in non-equilibrium. One solution is
not to reset the grid at the end of the time step. However doing so can result in mesh distortion
which makes MPM resemble to FEM. The cell-crossing error can be mitigated using high order
B-splines basis functions or the generalized interpolation material point (GIMP) method.
These functions, which can be interpreted as averages of the basis function and its gradient
over the particle domain, are used for mapping and interpolating data between particles and grid
nodes as well as calculating internal and external forces.

1
Z
φIp ≡ φI (xp ) = S̄Ip = χ(x − xp )NI (x)dΩ (47)
Vp Ωχ
1
Z
∇φIp ≡ ∇φI (xp ) = ∇S̄Ip = χ(x − xp )∇NI (x)dΩ (48)
Vp Ωχ

where χ(x) is the particle characteristic function centered at the particle position and Ωχ denotes
the current support of this function; NI are the standard linear shape functions. Typically piece-
wise constant particle characteristic function is used
(
1 if x ∈ Ωp
χ(x) = (49)
0 otherwise
More precisely in 1D the following particle characteristic function is usually employed
(
1 if |x| ≤ lp /2
χ(x)p = (50)
0 otherwise
Here lp is the current particle size. The initial particle size is determined by dividing the cell
spacing L by the number of particles per cell. In order to have a closed form expression for
the weighting functions given in Equation (48) so as to avoid numerical integration of these
functions, it is assumed that Ωp is always a rectangle which is aligned with the grid. In this
case, one can explicitly compute the weighting function φ as

32
8 GENERALIZED INTERPOLATION MATERIAL POINT METHOD-GIMP

Figure 21: Impact of two elastic bodies: simulation snapshots.

33
 9 IMPLICIT TIME INTEGRATION

node

(a)

(equilibriuum)
particle

(b)

Figure 22: Cell crossing issue in MPM.

1 − (4x2 + lp2 )/(4hlp )



 if |x| < 0.5lp

1 − |x|/h if lp ≤ |x| ≤ h − 0.5lp
φ(x) = (51)


(h + lp /2 − |x|)2 /(2hlp ) if h − 0.5lp ≤ |x| < h + 0.5lp
0

otherwise
of which an example is given in Fig. 23. Note that φI = φ(x − xI ). As can be seen, if there are
many particles per element (lp is getting smaller), the GIMP functions resemble the standard
FE hat functions.
Since the GIMP functions depend on the particle size which in turn evolves in time, the
GIMP functions are particle specific and time-dependent. To simplify things, Ωp is kept fixed
and the resulting approach is referred to as uGIMP (unchanged GIMP). A more complicated
approach, known as cpGIMP (contiguous particle GIMP), is to update the particle domain using
the deformation gradient F, cf. Fig. 24. Still, the updated particle domain is a rectangle in 2D
and space cannot be tiled.
The implementation of uGIMP follows closely the MPM except two things: (i) the GIMP
functions φI replace the FE hat functions NI and (ii) a larger connectivity array of the elements
i.e., the particles not only contribute to the nodes of the element they locate but also nodes of
neighboring elements.

χ(x)p = Vp δ(x − xp ) (52)

9 Implicit time integration

The explicit code operates with a time step restricted by the usual Courant-Friedrichs-Lewy
(CFL) condition, and thus the code resolves all waves. Simulating many manufacturing and

34
9 IMPLICIT TIME INTEGRATION

1
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3 3.5 4
node

1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
0 0.5 1 1.5 2 2.5 3 3.5 4

1
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3 3.5 4
Figure 23: One dimensional GIMP basis functions.

Figure 24: Tracking particle domain in GIPM: space cannot be tiled in a general multi-
dimension domain using rectilinear Ωp .

35
 10 FLUID-STRUCTURE INTERACTION

5 5

1 2

(a) MPM (b) GIMP

Figure 25: Larger element connectivity in GIMP (b) than in MPM (a). The connectivity of an
element in GIMP can be determined by getting firstly its neighboring elements (very fast for
structured grids used in MPM/GIMP) and then the nodes of these elements.

other problems with an implicit method can result in a considerable reduction in computational
cost if one is interested only in the low-frequency, bulk motion, and not in the elastic wave
motion.

10 Fluid-structure interaction
• Action of the air on aeronautic structures (aeroelasticity, flutter)

• Effect of the wind on civil structures (bridges, suspended cables, skyscrapers)

• Effect of water movement on dam or sloshing of a fluid in a container

• Fluid-membrane interaction (parachutes, vehicle airbags, blood vessels etc.)

The numerical procedures to solve the Fluid-structure interaction (FSI) problems may be
broadly classified into two approaches

• Monolithic approach: the equations governing the fluid flow and the structure deforma-
tion are solved in a single solver;

• Partitioned approach: the fluid and solid equations are solved separately using different
solvers.

The monolithic approach requires a code developed for this particular FSI problem whereas the
partitioned approach preserves software modularity because an existing flow solver and solid
solver are coupled. On the other hand, development of stable and accurate coupling algorithms
is required in partitioned approaches. Particularly, the fluid-solid interface location is not known
a priori and usually changed in time; thus, the partitioned approach requires the tracking of the
new interface location and its related quantities, which can be cumbersome and error-prone.

36
 REFERENCES REFERENCES

The partitioned approach is also referred to as coupling method as there are two models
which are coupled. Within the coupling method, one can differentiate between one-way and
two-way coupling. The coupling is one-way if the motion of a fluid flow influences a solid
structure but the reaction of a solid upon a fluid is negligible. The other way around is also
possible. Moreover, the coupling can be loose or tight coupling.
2
σ f = 2µǫ̇ − µtr(ǫ̇)I − p̂I (53)
3

p̂ = (γ − 1)ρe (54)
e is the specific internal energy and γ is the ratio of specific heats.

References
S.G. Bardenhagen. Energy Conservation Error in the Material Point Method for Solid Mechan-
ics. Journal of Computational Physics, 180(1):383–403, 2002.

M Steffen, PC Wallstedt, and JE Guilkey. Examination and analysis of implementation choices

within the material point method (MPM). Computer Modeling in Engineering and Sciences,
31(2):107–127, 2008.

D Sulsky and HL Schreyer. Axisymmetric form of the material point method with applications
to upsetting and Taylor impact problems. Computer Methods in Applied Mechanics and
Engineering, 0457825(96), 1996.

D Sulsky, Z Chen, and HL Schreyer. A particle method for history-dependent materials. Com-
puter Methods in Applied Mechanics and Engineering, 5, 1994.

D. Sulsky, S.J. Zhou, and H. L. Schreyer. Application of a particle-in-cell method to solid

mechanics. Computer Physics Communications, 87(1-2):236–252, 1995.

D. Sulsky, H.L. Schreyer, K. Peterson, R. Kwok, and M. Coon. Using the material-point method
to model sea ice dynamics. Journal of Geophysical Research, 112(C2):C02S90, 2007.

AR York, D. Sulsky, and HL Schreyer. The material point method for simulation of thin
membranes. International Journal for Numerical Methods in Engineering, 1456(May 1998),
1999.

AR York, D. Sulsky, and HL Schreyer. FluidâĂŞmembrane interaction based on the material

point method. International Journal for Numerical Methods in Engineering, pages 901–924,
2000.

37

View publication stats