What is Peridynamics?
Pablo Castillo Usyaopin
Computational Engineering
Overview
1 Introduction
2 Basics of Peridynamics Theory
3 A glance to the coded Peridynamic equation of motion
Pablo Usyaopin (Computational Engineering) Peridynamics++ 2 / 29
Motivation
Classical continuum mechanics have
successfully been applied to simulate
a variety of material types.
Theory based on partial differential
equations to compute stresses and
strains.
Among most commonly numerical
approaches: finite element and finite
difference methods.
Pablo Usyaopin (Computational Engineering) Peridynamics++ 3 / 29
Motivation II
But ...
An important assumption of a
continuous body is made, which is not
always necessarily true!.
Solving its equations rely on a
differentiation process (local), well
suited for analytical mathematics, but
in computational mathematics an
integral approach might be prefered
(non-local).
Pablo Usyaopin (Computational Engineering) Peridynamics++ 4 / 29
Motivation III
Objective:
Focuse on Peridynamics and its
applications:
Simulating cracks propagation.
Branching phenomena.
Damage.
Introduce Peridynamics C++ code
(Peridigm).
Hands on implementation of
material/damage models.
Pablo Usyaopin (Computational Engineering) Peridynamics++ 5 / 29
Basics of Peridynamics I
Peridynamics (PD) is a reformulation of classical elasticity where stress and strain
partial derivatives are substituted by force and displacement integral equations.
PD Equation ofZ motion:
ρ(x)ü(x, t) = T x, t q − x − T q, t x − q dVq + b(x, t)
Hx
· → notation indicating the bond in which the force state T is acting on.
q − x → is the bond ξ formed by particles q and x.
Pablo Usyaopin (Computational Engineering) Peridynamics++ 6 / 29
Basics of Peridynamics II
Z
ρ(x)ü(x, t) = T x, t q − x − T q, t x − q dVq + b(x, t)
Hx
Hx → spherical neighbourhood of radius δx centered on particle x.
ρ(x) → Density of particle x.
ü(x, t) → Acceleration of particle x in time t.
dVq → Volume of particle x.
b(x, t) → External body force on particle x in time t.
Pablo Usyaopin (Computational Engineering) Peridynamics++ 7 / 29
Peridynamics and its relation with Classical Continuum Equations
Relation Peridynamic Theory Classical Theory
∂y
Kinematics Y q − x = y(q) − y(x) F= (x)
Z ∂x
Linear Momentum ρü(x) = T x, t q − x − T q, t x − q dVq + b(x) ρü(x) = O · σ(x) + b(x)
Hx
Constitutive Model T = T̂(Y) σ = σ̂(F)
Z
Angular
Pablo Momentum
Usyaopin (Computational Engineering) Y q − x Peridynamics++
× T q − x dVq = 0 σ = σT 8 / 29
How is damaged introduced?
The simplest way is through bond breakage:
|Y| − |X| e
s= = .
|X| |X|
Y
T = (1 − d) t M; M=
|Y|
Particle Damage:
Z
d(ξ)dV q
Hx
φ(x) = Z .
dV q
Hx
Pablo Usyaopin (Computational Engineering) Peridynamics++ 9 / 29
A glance to the coded Peridynamic equation of motion
Z
ρ(x)ü(x, t) = f (q, x)dVq + b(x),
Hx
in which:
f (q, x) = T x, t q − x − T q, t x − q .
Pablo Usyaopin (Computational Engineering) Peridynamics++ 10 / 29
Bond Based Peridynamics vs State Based Peridynamics I
Bond Based:
Z
ρ(x)ü(x, t) = f (q, x)dVq + b(x),
Hx
Y 18k |Y| − |X| y−x e
f (q, x) = c s ; c= ; s= = =
|Y| πδx4 |X| x x
State based:
Z
ρ(x)ü(x, t) = T x, t q − x − T q, t x − q dVq + b(x, t)
Hx
Z Z
3
θ= (ω x) · e dVq ; m = (ω x) · x dVq
m Hx Hx
3kθ 15µ d Y θx
T= ωx+ ωe ; ed = e −
m m |Y| 3
Pablo Usyaopin (Computational Engineering) Peridynamics++ 11 / 29
Diff. between Bond Based and State Based Peridynamics II
Shear modulus in term of Bulk and Poisson ratio:
" #
3k(1 − 2v) 3k
µ= =
2(1 + v) 5
v=1/4
Substituting into T:
" !#
3kθ 9k θx Y
T= ωx+ ω e−
m m 3 |Y|
" #
9k Y
T= ωe
m |Y|
1
Setting ω =and evaluating m in spherical coordinates → m = πδx4 .
x
" ! #
9k e Y Y
T= = cs
πδx4 x |Y| |Y|
Pablo Usyaopin (Computational Engineering) Peridynamics++ 12 / 29
The Linear Peridynamic Solid Material (LPS) in Peridigm
State based:
" #
3kθ 15µ d Y θx
T= ωx+ ωe ; where ed = e −
m m |Y| 3
" !#
3kθ 15µ θx Y
T= ωx+ ω e− ;
m m 3 |Y|
" ! ! #
15µ
3k 15µ Y
T = ωθ − m x+ω e ;
m 3 m |Y|
" #
Y
T = ω θ c1 x + ω α e ;
|Y|
Pablo Usyaopin (Computational Engineering) Peridynamics++ 13 / 29
Example: Tensile Test
Material Model: LPS
(State-Based)
Bulk Modulus (k):
1.5 × 1012 P a
Shear Modulus (µ):
6.923 × 1011 P a
Poisson ratio (v): 0.3
No. of Particles: 12, 625
Pablo Usyaopin (Computational Engineering) Peridynamics++ 14 / 29
Example: Tensile Test II
Time (min)
Material Model
1 Proc 4 Proc
LPS (State Based PD) 33.3 12.2
Bond Based PD 30.0 11.0
Pablo Usyaopin (Computational Engineering) Peridynamics++ 15 / 29
Bond Based Peridynamics vs State Based Peridynamics?
Bond based presented model limits the range of materials to Poisson ratio v = 1/4
(3D) or Plane Strain and v = 1/3 (Plane Stress) due do the independent
deformation computed for every bond. The model also employs one material
parameter (k).
State Based approach (Linear Peridynamic Solid) considers the collective
deformation of all bonds, therefore allows to model a bigger range of Poisson ratios.
In addition the model employs two material parameters (k and µ).
Pablo Usyaopin (Computational Engineering) Peridynamics++ 16 / 29
Classical Kinematics
Initial Configuration:
dx = x(q) − x(p)
Deformed Configuration:
dy = y(q) − y(p)
dy = φ((x(p) + dx, t)) − φ(x(p), t)
where: y(p) = φ(x(p));
dy = Fdx
y(q) = φ(x(q))
In many textbooks: y = y(x, t)
∂φ
∂y Defining: F = = ∇0 φ
Thus: F = ∂x
∂x
Pablo Usyaopin (Computational Engineering) Peridynamics++ 17 / 29
Peridynamic Kinematics
Initial Bond Configuration:
ξ = x(q) − x(p)
Deformed Bond Configuration:
Y ξ = y(q) − y(p)
Y ξ = x(q) − x(p) + u(q) − u(p) δp → Radius of interaction of
particle p with other particles.
Y ξ =ξ+η
y(p) = x(p) + u(p);
No continuous function assumption is done between
y(q) = x(q) + u(q)
initial and deformed configuration.
Pablo Usyaopin (Computational Engineering) Peridynamics++ 18 / 29
Classical Continuum Strains 1
Do you remember the classical strain
l −l
learned in school ε = f l0 0 ??
The true is that there are maaaany
strains!!!, but they all share a common
root: All can be derived from the
deformation gradient F.
Example 1: Right Cauchy Strain Example 2: Left Cauchy Strain
dx1 · dx2 = F−1 dy1 · F−1 dy2
dy1 · dy2 = F dx1 · F dx2
dx1 · dx2 = dy1 · (F−1 )T F−1 dy2
T
dy1 · dy2 = dx1 · F F dx2 −1
dx1 · dx2 = dy1 · (FFT ) dy2
dy1 · dy2 = dx1 · C dx2 −1
dx1 · dx2 = dy1 · b dy2
Pablo Usyaopin (Computational Engineering) Peridynamics++ 19 / 29
Classical Continuum Strains 2
Example 3: Green Lagrange Strain
1 1
dy1 · dy2 − dx1 · dx2 = dx1 · C dx2 − dx1 · dx2
2 2
1 1 h i
dy1 · dy2 − dx1 · dx2 = dx1 · C dx2 − dx2
2 2
1 h1 i
dy1 · dy2 − dx1 · dx2 = dx1 · C − I dx2
2 2
1
dy1 · dy2 − dx1 · dx2 = dx1 · E dx2
2
Pablo Usyaopin (Computational Engineering) Peridynamics++ 20 / 29
What about Peridynamic Strains?
There are no tensorial strains!,
however, it is still possible to
compute a deformation
gradient F using integrals!...
Z
F= ω ξ Y ξ ⊗ X ξ dVq · K−1 ,
Hx
Z
K= ω ξ X ξ ⊗ X ξ dVq .
Hx
Pablo Usyaopin (Computational Engineering) Peridynamics++ 21 / 29
Damage in Peridynamics
bond stretch
|Y| − |X| e
s= = .
|X| |X|
bond damage
(
1 if s > s0 ,
d(ξ) =
0 otherwise.
particle damage
Z
d(ξ)dV q
Hx bond force
φ(x) = Z .
dV q Y
Hx T = (1 − d) t M; M=
|Y|
Pablo Usyaopin (Computational Engineering) Peridynamics++ 22 / 29
Frechet Derivative
Definition:
Ψ(A + ∆A) = Ψ(A) + ∇A Ψ • ∆A + O(|∆A|)
Function Increment → ∆Ψ = ∇A Ψ • ∆A + O(|∆A|)
Example 1: Lets say we have: → Ψ(Y) = Y • Y
Ψ(Y + ∆Y) = (Y + ∆Y) • (Y + ∆Y)
= Y • Y + 2Y • ∆Y + ∆Y • ∆Y
Function Increment → ∆Ψ = 2Y • ∆Y
∂Ψ
Frechet Derivative → ∇Y Ψ = 2Y or = 2Y
∂Y
Pablo Usyaopin (Computational Engineering) Peridynamics++ 23 / 29
Increment in Bond Extension
Example 2: Consider the bond extension: e = |Y| − |X|
e(Y + ∆Y) = |Y + ∆Y| − |X|
Function Increment → ∆e = |Y + ∆Y| − |X| − |Y| + |X|
→ ∆e = |Y + ∆Y| − |Y|
Pablo Usyaopin (Computational Engineering) Peridynamics++ 24 / 29
Increment in Dilatation
3
Example 2: Dilatation (volumetric deformation): → θ = ωx•e
m
3
θ(e + ∆e) = ω x • (e + ∆e)
m
3 3
= ω x • e + ω x • ∆e
m m
3
Function Increment → ∆θ = ω x • ∆e
m
3 ∂θ 3
Frechet Derivative → ∇e θ = ωx or = ωx
m ∂e m
Pablo Usyaopin (Computational Engineering) Peridynamics++ 25 / 29
Increment in Deviatoric Bond Extension
x
Example 2: Deviatoric Bond extension: → ed = e − θ
3
x
Function Increment → ∆ed = ∆e − ∆θ
3
x
→ ∆ed = ∆e − ∇e θ • ∆e
3
Pablo Usyaopin (Computational Engineering) Peridynamics++ 26 / 29
Additional Properties
Property regarding Nested functions:
∇A (Φ(Ψ(A))) = Φ0 (Ψ(A)) ∇A (Ψ(A))
which can also be expressed as:
!
∂ ∂Φ ∂Ψ
Φ(Ψ(A)) =
∂A ∂Ψ ∂A
Pablo Usyaopin (Computational Engineering) Peridynamics++ 27 / 29
Linear Peridynamic Solid Material
Energy Density
∆W = ∇W • ∆y
which can also be expressed as:
!
∂ ∂Φ ∂Ψ
Φ(Ψ(A)) =
∂A ∂Ψ ∂A
Pablo Usyaopin (Computational Engineering) Peridynamics++ 28 / 29
Linear Peridynamic Solid Material
Energy Density
kθ2 α
W = + (ωed ) • ed )
2 2
∆W = kθ∇e θ • ∆e + α(ωed ) • ∆ed
x
= kθ∇e θ • ∆e + α(ωed ) • ∆e − ∇e θ • ∆e
3
h
d
x i
= kθ∇e θ + α(ωe ) • 1 − ∇e θ • ∆e
3
h α i
= kθ − ωx • e ∇e θ + α(ωed ) • ∆e
d
3
h α n θx o 3 i
= kθ − ωx • e − ωx + α(ωed ) • ∆e
3 3 m
h α n θm θωx • x o 3 i
= kθ − − ωx + α(ωed ) • ∆e
3 3 3 m
h 3kθ i
= ωx + α(ωed ) • ∆e
m
Pablo Usyaopin (Computational Engineering) Peridynamics++ 29 / 29