Material point method and it’s evolution over years

Ankit kumar
Under the Supervision of
Dr. Rajneesh Sharma
Associate Professor

Indian Institute of Technology Mandi

School of Civil and Environmental Engineering (SCENE), IIT Mandi
1
Introduction
Material domain is
represented by
assembly of particles
Grid particle hybrid
method
Boundary conditions
are applied on grids

Mesh tangling never

occurs
Numerical methods Control domain
Particle based
based on Continumm defined by
methods
mechanics interpolation functions
Special handling of
boundaries is required

Highly applicable to
boundary value
problems
Define the element
Mesh based methods domain using nodes
and their connectivity
Mesh tangling due to
strict boundary
definitions 2
Eulerian and Lagrangian description
• In the Eulerian description of fluid • In the Lagrangian description of fluid
motion ,we consider infinitesimal fluid motion ,we consider infinitesimal
element fixed in space with the fluid fluid element moving along a
flowing through it. streamline.

Source https://www.researchgate.net/figure/The-Eulerian-and-
Lagrangian-ways-to-describe-the-flow 3
Flow chart showing evolution of material
point method

Arbitrary particle
Improved version of Improved version of
domain interpolation
PIC METHOD MPM
• PARTICLE IN • Material point • Convected MPM
CELL METHOD • Fluid implicit method • Generalised particle domain • Improvised
particle interpolated interpolation version of
method material point MPM MPM
To study dynamic To apply FLIP in solid Improved version of
compressible fluids mechanics gimp

4
Particle in cell method ( Martha W. Evans &
Francis H. Harlow in 1957)
• A method is presented for solving hydrodynamic problems involving large
distortions and compressions of the fluid in several space dimensions.
Lagrangian
1. Idea : particles
• Particles handle advection trivially
• Grids handle interactions efficiently
2. Put these two together
• Transfer quantities to grids grid
• Solve on grids(interaction forces)
• Transfer back to particles
• Move particles(advection)
Source https://warpx.readthedocs.io/en/latest/theory/pic 5
Flow chart for pic method
Start with
particles(1)

Transfer to grid(2)

Move particles(5)

Resolve forces on grid(3)

Transfer velocity back to

particles(4)
6
Problems with particle in cell method
More numerical dissipation than pure Eulerian method.
We have to resample (average) twice.
When we average from particles to grid, simple weighted averages is only
first order.
Face challenges in handling boundary conditions.
Low accuracy in collision models.

7
Fluid implicit particle method (J.U. BRACKWILL ,D.B
KOTHE &H.M. RUPELL IN 1986)
Changes they implement
Transfer back the change of a quantity from grid to particles rather than the qualities itself.
Advantages
Each delta only averaged once
Nearly eliminated numerical dissipation
Improved accuracy

8
Why the requirement of material point method
FLIP can only handle the fluid to fluid interactions.
There is only interaction of fluid to solid at boundaries.
MPM is effective for simulating solid materials.

9
Material point method (D.SULSKY,Z.CHEN
&H.L. SCHREYER IN 1994)
1. Langragian material points carry all the state data(position, velocity, stress etc).

2. Overlying mesh defined.

3. Particle state projected to mesh.

10
Continue…
4. Conservation of momentum solved on mesh giving updated mesh velocity and
position.

5. Particle position and velocities updated from mesh solution.

6. Discard deformed mesh and define new mesh and repeat.

11
MPM codes approach
Particle
initialization
Update particle Particle data to
stress grid

Compute Boundary
particle stress conditions

Boundary Compute
conditions internal forces

Integrate nodal Compute

velocities acceleration 12
Step 1 (particle initialization)
• There are two basic method for particle locations
1. Acquire them from file e.g. image data

2. Use geometric primitives to describe geometry

and tests to determine particle placement

At time=0 , they must have

Position
Volume
Mass
Velocity
Stress
Deformation gradient
Sourcehttps://www.google.com/imgres?q=material%po 13
size int%method images
Step 2 (particle data to grid)
 Moment conservation equation
= + (1)
 To calculate the quantities at the nodes, from each material point and each grid node, the algorithm defines the
weighting function Svp, as:

(2)
 Having defined the weighting function, the nodal quantities, such as nodal momentum pv and nodal masses are
computed as the product of weighting function and material point momentum Pp:
(3)
(4)

14
Step 3 (boundary conditions)
• For simplicity, assume that the computational domain is a rigid box.
• If the velocity of the rigid walls is zero, then set the velocity on those
computational nodes to be zero.
• Also need to set the velocity on the "extra" nodes to be zero as well.
Domain boundary Extra nodes

15
Step 4 (compute internal forces)
• The internal forces for given grid node due to stress in the material points are
(5)

• If the body forces b is taken into account

(6)

Step 5 (compute acceleration)

• This is just basically inverting Newton’s Second law to get acceleration at each grid node :

(7) 16
Step 6 (integrate nodal velocities)
Using basic forward Euler integration, advance the velocities at the grid nodes

(8)

Step 7 (boundary conditions)

Any Dirichlet boundary conditions (e.g. specified velocity and in particular velocity
equal to zero), it is applied here by overwriting the computed grid momentum.

Step 8 (compute particle stress)

The nodal velocities gradients can be used to compute the rate of change in the field of strain at
any point x of the domain :
= (9)
17
By using consecutive relationship we can compute stresses
Step 9 (update particle state)

Update particle velocity :

(10)
Update particle location :
(11)
Update the time
(12)
Step 10 (return to step 2)
Repeat step 2 to step 9 until desired result were obtained.

18
Problems with material point method
When a particle crosses the grid boundary, the
force at the node due to particles stress changes
sign, as the shape function gradient changes
sign.
Cell crossing error

Linear shape function (left) and its

gradient discontinuous at the node.
Representation of cell
crossing errors
19
Generalised interpolated material point method
(Bardenhagen &Kober in 2004)
Also known as Contiguous Particle GIMP
method

Changes
Shape function constructed by
MPM & GIMP Interpolation function
Integrating linear FEM shape function.

Results :
Less sensitive to mesh distortion
Improved accuracy
MPM &GIMP shape function’s gradient

Source : Bardenhagen and Kober ,2004 20

Problems with GIMP
• curved geometry is described by a ragged stair-step edge.
• the effect of shear deformation on the particle domain is not accounted.
• at large deformation, a numerical fracture is still possible as a gap between
particles develops in tensile loading.

21
Convected particle domain interpolation method(A.
Sadeghirad,R.M. Bannon & J. Burghardt in 2011)
In CPDI method, a novel alternative to
standard grid basis functions is proposed
By using these alternative grid basis
functions, dividing particle domains along
cell boundaries is avoided.
results
CPDI method is capable of more accurately
tracking particle domains with a
computational cost comparable with original
GIMP.

Source : Sadeghirad et. al 2020 22

Problems in CPDI
• particle domains are tracked as parallelograms in 2 ‐D (or
parallelepipeds in 3‐D).
• In the large deformation regime, the CPDI still may not retain
the expected convergence rate.
• It’s solution are not converging at second order.

Source 23
A Convected particle least square interpolation material
point method (Quoc-Anh Tran & Wojciech Solowski in
2020)
 It,s solution converges at second order.

 B-spline MPM can reduce grid-crossing

errors.

 Higher converging rate.

 Retains the advantages of CPDI method

QUADRATIC B-splines
Graph from:Moving least squares reconstruction for
B-spline Material Point Method

Source : (Tran and Solowski, 2020) 24

Arbitrary particle domain interpolation Material point
method (Takatoshi Kiriyama &Yosuke Higo in 2020)

Changes they implement in CPDI

APDI makes any type of control domain
deformation possible.
Uses numerical integration for the domain Quadrilateral & Hexahedron control domain used in
integration for any arbitrary shapes CPDI(2)

Effects
Capable of simulating the shape of a target
domain
Providing a precise trace of deformation.

Arbitrary Triangle & Tetrahedron control domain in

25
APDI
 Material
GIMP CPDI APDI
point method
Represent control uniform/unchanged CPDI(1) Control domain
domain as points GIMP • Control domain deforms as deforms any arbitrary
• Numerical oscillation • Control domain regain its parallelogram in 2D shape
initial shape • Numerical integration for
domain integration

Contiguous particle CPDI(2)

GIMP • Control domain deform as
• Control domain is updated quadrilateral
as it deform • Use direct integration for
deformed domains

26
Numerical example result :

Analysis of Vibration of bar without damping with the help of MPM compared with
analytical result

27
Future scope:

1.Challenges remain in handling accurate contact laws and complex

boundary conditions.
2.Lack of user-friendly interfaces for MPM software and fragmentation
within the MPM community.

28
THANK YOU

29