DISCRETIZED BOND-BASED PERIDYNAMICS FOR SOLID MECHANICS
By
Wenyang Liu

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Civil Engineering
2012

ABSTRACT
DISCRETIZED BOND-BASED PERIDYNAMICS FOR SOLID MECHANICS
By
Wenyang Liu

The numerical analysis of spontaneously formed discontinuities such as cracks is a longstanding challenge in the engineering field. Approaches based on the mathematical framework of classical continuum mechanics fail to be directly applicable to describe discontinuities
since the theory is formulated in partial differential equations, and a unique spatial derivative, however, does not exist on the singularities. Peridynamics is a reformulated theory of
continuum mechanics. The partial differential equations that appear in the classical continuum mechanics are replaced with integral equations. A spatial range, which is called the
horizon δ, is associated with material points, and the interaction between two material points
within a horizon is formed in terms of the bond force. Since material points separated by a
finite distance in the reference configuration can interact with each other, the peridynamic
theory is categorized as a nonlocal method.
The primary focus in this research is the development of the discretized bond-based
peridynamics for solid mechanics. A connection between the classical elasticity and the discretized peridynamics is established in terms of peridynamic stress. Numerical micromoduli
for one- and three-dimensional models are derived. The elastic responses of one- and threedimensional peridynamic models are examined, and the boundary effect associated with the
size of the horizon is discussed. A pairwise compensation scheme is introduced in this research for simulations of an elastic body of Poisson ratio not equal to 1/4. In order to enhance
the computational efficiency, the research-purpose peridynamics code is implemented in an

NVIDIA graphics processing unit for the highly parallel computation. Numerical studies are
conducted to investigate the responses of brittle and ductile material models. Stress-strain
behaviors with different grid sizes and horizons are studied for a brittle material model. A
comparison of stresses and strains between finite element analyses and peridynamic solutions
is performed for a ductile material. To bridge material models at different scales, a multiscale
procedure is proposed.
An approach to couple the discretized peridynamics and the finite element method is developed to take advantage of the generality of peridynamics and the computational efficiency
of the finite element method. The coupling of peridynamic and finite element subregions is
achieved by means of interface elements. Two types of coupling schemes, the VL-coupling
scheme and the CT-coupling scheme respectively, are introduced. Numerical examples are
presented to validate the proposed coupling approach including one- and three-dimensional
elastic problems and the mixed mode fracture in a double-edge-notched concrete specimen.
A numerical scheme for the contact-impact procedure ensuring compatibility between a
peridynamic domain and a non-peridynamic domain is developed. A penalty method is used
to enforce displacement constraints for transient analyses by the explicit time integration.
In the numerical examples, the impact between two rigid bodies is investigated to validate
the contact algorithm. The ballistic perforation through a steel plate is investigated, and the
residual velocities of the projectile are compared with the results by an analytical model.
Peridynamics is applied to study porous brittle materials. An algorithm is developed
to generate randomly distributed cubic voids and spherical voids for a given porosity. The
material behaviors at the macroscopic level including the resultant Young’s modulus and
the strength are studied with varying amounts of porosities. The degradations of Young’s
modulus and strength are compared with empirical and analytical solutions.

To my mother, without whom I could not have completed this work

iv

ACKNOWLEDGMENT

I would like to express my deepest appreciation to all the people who have inspired and
helped me during the course of my studies. First of all, I would like to express my gratitude
to my advisor, Professor Jung-Wuk Hong. I appreciate the opportunity he gave me to study
computational mechanics at Michigan State University. I thank him for his guidance, support, and countless hours of insightful discussion. I would like to thank Professors Venkatesh
Kodur, Nizar Lajnef, and Xinran Xiao for joining my Ph.D. guidance committee, for their
valuable time, and for their advice throughout my education at Michigan State University.
A special appreciation is given to the Department of Civil and Environmental Engineering
at Michigan State University. I thank the department for all their assistance, for supporting
me to attend a conference, and for my fellowship. I am thankful to all instructors with whom I
took courses at MSU, Professors Seungik Baek, Xinran Xiao, Thomas Pence, Giles Brereton,
Venkatesh Kodur, Jung-Wuk Hong, and Parviz Soroushian. The knowledge I learned from
them has provided me a solid background to conduct research. I thank the staff at High
Performance Computing Center at MSU for helping me using their computational resources.
I thank Zeynep Altinsel for teaching me speaking skills and all the encouragement she gave
me. I thank all of my friends for their support.
The research in this dissertation has been supported by the U.S. Air Force Office of
Scientific Research (Grant FA 9550-10-1-0222) and the MSU startup fund of Professor JungWuk Hong.

v

TABLE OF CONTENTS

List of Tables . . . . . . . . . . . . . . .
List of Figures

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . xvi

List of Symbols . . . . . . . . . . . . . . .
Chapter 1
Introduction . . . . . . . . . .
1.1 Overview of peridynamics
1.2 Research Objectives . . . .
1.3 Scope . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

. . . . . . . . . . . . . . . . . . xvii

.
.
.
.

1
1
4
5

Chapter 2
Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 Numerical predictions of crack growth . . . . . . . . . . . . . . . . . . . . .
2.2 Peridynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
10

Chapter 3
Discretized Peridynamics for Linear Elastic Solids . . . . . . . . . . . . . .
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Peridynamic formulation . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Micromodulus of elastic materials . . . . . . . . . . . . . . . . . . . .
3.2.2.1 One-dimensional model . . . . . . . . . . . . . . . . . . . .
3.2.2.2 Three-dimensional model . . . . . . . . . . . . . . . . . . .
3.2.3 Pairwise compensation scheme . . . . . . . . . . . . . . . . . . . . . .
3.3 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Comparison of peridynamics and analytical solution(ν = 1 ) . . . . . .
4
3.3.2.1 One-dimensional bar . . . . . . . . . . . . . . . . . . . . . .
3.3.2.2 Three-dimensional bar . . . . . . . . . . . . . . . . . . . . .
3.3.3 Comparison of modified peridynamics and analytical solution (ν = 1 )
4
3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14
14
15
15
20
21
23
27
30
30
31
31
32
39
42

vi

.
.
.
.

ix

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

Chapter 4
Discretized Peridynamics for Brittle and Ductile Solids
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Computation . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Numerical implementation . . . . . . . . . . . . . .
4.3.2 Implementing Peridynamics in GPU . . . . . . . .
4.4 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Brittle material . . . . . . . . . . . . . . . . . . . .
4.4.2 Ductile material . . . . . . . . . . . . . . . . . . . .
4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

Chapter 5
Coupling of Discretized Peridynamics with Finite Element Method
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Finite element formulations . . . . . . . . . . . . . . . . . . . .
5.3 Coupling between the peridynamic and finite element subregions . . . .
5.3.1 Coupling schemes . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.2 Inverse isoparametric mapping . . . . . . . . . . . . . . . . . . .
5.4 Numerical applications . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1 One-dimensional bar . . . . . . . . . . . . . . . . . . . . . . . .
5.4.2 Three-dimensional bar . . . . . . . . . . . . . . . . . . . . . . .
5.4.3 Mixed mode fracture in a tension-shear specimen . . . . . . . .
5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

. 71
. 71
. 73
. 73
. 76
. 76
. 81
. 85
. 85
. 88
. 97
. 103

Chapter 6
Discretized Peridynamics for Contact-Impact
6.1 Introduction . . . . . . . . . . . . . . . . . . .
6.2 Theory . . . . . . . . . . . . . . . . . . . . . .
6.2.1 Peridynamic short-range forces . . . .
6.2.2 Inverse isoparametric mapping . . . . .
6.2.3 Contact algorithm . . . . . . . . . . .
6.2.4 Penalty method . . . . . . . . . . . . .
6.3 Ballistic limit . . . . . . . . . . . . . . . . . .
6.4 Numerical studies . . . . . . . . . . . . . . . .
6.4.1 Material behavior and modeling . . . .
6.4.2 Rigid-body impact . . . . . . . . . . .
6.4.3 Ballistic perforation . . . . . . . . . . .
6.5 Summary . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

Problems
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

43
43
44
45
45
45
53
53
61
69

105
105
110
110
111
112
114
116
118
118
123
126
134

Chapter 7
Modeling of Porous Brittle Solids using Peridynamics . . . . . . . . . . . 136
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
7.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
vii

7.3
7.4

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

142
143
143
148
156

Chapter 8
Conclusions and Recommendations .
8.1 Conclusions . . . . . . . . . . . . . .
8.2 Contributions . . . . . . . . . . . . .
8.3 Recommendations for future research

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

157
157
161
162

7.5

Numerical implementation
Case studies . . . . . . . .
7.4.1 Cubic voids . . . .
7.4.2 Spherical voids . .
Summary . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

Appendix A
Derivation of Three-dimensional Micromodulus . . . . . . . . . . . . . . . 166
Appendix B
Derivation of Modified Peridynamic Young’s Modulus . . . . . . . . . . . 169
Appendix C
Nonlinear Least Squares Method . . . . . . . . . . . . . . . . . . . . . . . . . 172
Appendix D
Newton-Raphson Method for Multi-dimensional Nonlinear Systems . . 175
Bibliography . . . . . . . . . . . . . . .

viii

. . . . . . . . . . . . . . . . . . . 178

LIST OF TABLES

Table 3.1

Micromodulus c1 in one-dimensional domain. . . . . . . . . . . . . .

23

Table 3.2

Micromodulus c3 in three-dimensional domain. . . . . . . . . . . . .

26

Table 4.1

Dependent loops (Algorithm 1). . . . . . . . . . . . . . . . . . . . .

49

Table 4.2

Independent loops (Algorithm 2). . . . . . . . . . . . . . . . . . . .

50

Table 4.3

The wall-clock time to calculate bond forces for one time-step using
a node in the CPU. The three-dimensional peridynamic model has
4,725 nodes, the horizon δ = 2∆x and ∆x = 0.5 mm. . . . . . . . .

51

The wall-clock time to run peridynamic models for 500,000 timesteps on GPU and CPU, respectively. Peridynamic models have the
horizon δ = 2∆x and ∆x = 0.5 mm. . . . . . . . . . . . . . . . . . .

51

Table 4.4

Table 6.1

Explicit contact algorithm. . . . . . . . . . . . . . . . . . . . . . . . 116

Table 6.2

Numerical results of residual velocities. . . . . . . . . . . . . . . . . 127

Table 7.1

Algorithm to generate randomly distributed spherical voids. . . . . . 143

Table 7.2

Young’s modulus and tensile strength of specimens containing cubic
voids. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

Table 7.3

Young’s modulus and tensile strength of specimens containing spherical voids within the boundary. . . . . . . . . . . . . . . . . . . . . . 149

Table 7.4

Young’s modulus and tensile strength of specimens containing spherical voids intersecting with boundaries. . . . . . . . . . . . . . . . . 150

ix

LIST OF FIGURES

Figure 1.1

Schematic of peridynamics. (For interpretation of the references to
color in this and all other figures, the reader is referred to the electronic version of this dissertation.) . . . . . . . . . . . . . . . . . . .

2

(a) Relationships among relative position vector and the relative displacement vector within a peridynamic horizon. (b) Pairwise force
vector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

Figure 3.2

Discretized domain for computation. . . . . . . . . . . . . . . . . . .

19

Figure 3.3

Volume calculation scheme for discretized peridynamics. The volume
is reduced on the boundary of a horizon. . . . . . . . . . . . . . . .

19

Volumetric ratio in a horizon. The ratio decreases to 1/2 at the border
of a horizon [135]. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

Pairwise forces acting through a cross section of a node for δ = 2∆x
in a one-dimensional domain. Each node, represented by a sphere,
has a volume of (∆x)3 . . . . . . . . . . . . . . . . . . . . . . . . . .

22

Pairwise forces acting through a cross section of a node for δ = 2∆x
in a three-dimensional domain. Each node, represented by a sphere,
has a volume of (∆x)3 . . . . . . . . . . . . . . . . . . . . . . . . . .

24

Number of bonds acting through a cross section of a node in a threedimensional domain. . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

Figure 3.8

Micromodulus for different horizons (E = 70 GPa, and ν = 1/4). . .

26

Figure 3.9

Pairwise compensation scheme. . . . . . . . . . . . . . . . . . . . . .

27

Figure 3.10

Three-dimensional bar subjected to a tension.

29

Figure 3.1

Figure 3.4

Figure 3.5

Figure 3.6

Figure 3.7

x

. . . . . . . . . . . .

Figure 3.11

Three-dimensional bar subjected to a tension. Lateral forces are
added to take account of ν = 1 . . . . . . . . . . . . . . . . . . . . .
4

29

Figure 3.12

Strain distribution of the one-dimensional bar in x-direction. . . . .

32

Figure 3.13

Three-dimensional bar discretized by 21 × 15 × 15 nodes. . . . . . .

33

Figure 3.14

Strain distribution along the x-axis in the three-dimensional bar. (a)
Numerical micromodulus c3 given in Table 3.2 is used for each horizon
and (b) analytical micromodulus c3 given in Equation (3.14) is used.
˜
The classical (local) elasticity solution of strain εx = 0.005. . . . . .

34

(a) Strain- and (b) stress-distributions on the cross section at x =
5.0 mm. The classical (local) elasticity solution of strain εx = 0.005
and stress σx = 350 MPa. The horizon δ = 2∆x, and the numerical
micromodulus c3 is utilized in the simulation. . . . . . . . . . . . . .

35

Poisson ratio along the three-dimensional peridynamic bar as shown
in Figure 3.13. The given Poisson ratio of the material is ν = 1 . The
4
horizon δ = 2∆x, and the numerical micromodulus c3 is utilized in
the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

Strain distribution at the end surface of the bar for different horizons
using the numerical micromodulus: (a) δ = 2∆x, (b) δ = 3∆x, (c)
δ = 4∆x, and (d) δ = 5∆x. The classical (local) elasticity solution
of strain εx = 0.005. . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

(a) Strain εx and (b) Poisson ratio ν along the bar by the modified
peridynamic formulation. The classical (local) elasticity solution of
strain εx = 0.005, and the Poisson ratio ν of the material is 0.3. The
horizon δ = 2∆x, and the modified numerical micromodulus c3 for
ˆ
pairwise compensation scheme is utilized in the simulation. . . . . .

41

(a) Bond force as a function of bond stretch in a brittle material
model. (b) Bond force as a function of bond stretch in a ductile
material model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

Figure 4.2

NVIDIA Tesla Block Diagram [186]. . . . . . . . . . . . . . . . . . .

47

Figure 4.3

A benchmark problem of matrix multiplication. . . . . . . . . . . . .

48

Figure 3.15

Figure 3.16

Figure 3.17

Figure 3.18

Figure 4.1

xi

Figure 4.4

Figure 4.5

Figure 4.6

Figure 4.7

Figure 4.8

Figure 4.9

Figure 4.10

Figure 4.11

Comparison between the GPU enabled version of peridynamics code
and the serial version. (a) Wall-clock time and (b) the speedup of the
GPU version over the serial version. . . . . . . . . . . . . . . . . . .

52

(a) The geometry of a three-dimensional bar and the boundary conditions; (b) Displacement of the boundary region versus time. . . . .

54

Stress on the cross section at x = 5.0 mm for (a) different horizon δ
with the grid size ∆x = 0.5 mm and (b) different grid size ∆x with
the horizon δ = 2∆x. εx is the engineering strain imposed by the
displacement at both ends of the bar. . . . . . . . . . . . . . . . . .

56

Strain εx at time t = 1.0 × 10−3 s as denoted by I in Figure 4.6(b).
The horizon δ = 2∆x. (a) ∆x = 0.50 mm; (b) ∆x = 0.25 mm; (c)
the cross section at z = 3.5 mm, ∆x = 0.50 mm; (d) the cross section
at z = 3.5 mm, ∆x = 0.25 mm. . . . . . . . . . . . . . . . . . . . . .

57

Strain εx at time t = 1.6 × 10−3 s as denoted by II in Figure 4.6(b).
The horizon δ = 2∆x. (a) ∆x = 0.50 mm; (b) ∆x = 0.25 mm; (c)
the cross section at z = 3.5 mm, ∆x = 0.50 mm; (d) the cross section
at z = 3.5 mm, ∆x = 0.25 mm. . . . . . . . . . . . . . . . . . . . . .

58

Strain εx at time t = 1.8 × 10−3 s as denoted by III in Figure 4.6(b).
The horizon δ = 2∆x, the grid spacing ∆x = 0.50 mm. (a) bird’s-eye
view of bar; (b) the cross section at z = 3.5 mm; (c) the cross section
at x = 5.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

Strain εx at time t = 1.8 × 10−3 s as denoted by III in Figure 4.6(b).
The horizon δ = 2∆x, the grid spacing ∆x = 0.25 mm. (a) bird’s-eye
view of bar; (b) the cross section at z = 3.5 mm; (c) the cross section
at x = 5.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

A comparison between peridynamic and FEA solutions. (a) The
material model used in FEA which is identical to the constitutive
for bonds in peridynamics; (b) Stress σx on the cross section at
x = 5.0 mm; (c) Stress σx on an element which has a node at the
center of the bar (x = 5.0 mm, y = 3.5 mm and z = 3.5 mm) and
the stress σx on the peridynamic node at the center of the bar; (d)
Strain εx at the center of the bar (x = 5.0 mm, y = 3.5 mm and
z = 3.5 mm). ux is the displacement of the boundary region as
shown in Figure 4.5. . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

xii

Figure 4.12

Figure 4.13

Figure 5.1

Figure 5.2

Figure 5.3

Figure 5.4

A comparison between peridynamic and FEA solutions. (a) The modified material model used in FEA to incorporate retrieved macroscopic
material properties from the peridynamic model; (b) Stress σx on the
cross section at x = 5.0 mm; (c) Stress σx on an element which has
a node at the center of the bar (x = 5.0 mm, y = 3.5 mm and
z = 3.5 mm) and the stress σx on the peridynamic node at the center of the bar; (d) Strain εx at the center of the bar (x = 5.0 mm,
y = 3.5 mm and z = 3.5 mm). ux is the displacement of the boundary
region as shown in Figure 4.5. . . . . . . . . . . . . . . . . . . . . .

66

Strain εx when the boundary displacement ux = 0.025 mm. (a) LSDYNA and (b) peridynamics. . . . . . . . . . . . . . . . . . . . . . .

68

Partition of the domain. The FE subregion and the peridynamic
(PD) subregion are bridged by interface elements. . . . . . . . . . .

77

Interface element for the coupling of FE subregion and peridynamic
subregion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

VL-coupling scheme that divides a coupling force f cp among FE nodes
comprising the interface element. . . . . . . . . . . . . . . . . . . . .

80

CT-coupling scheme that divides a coupling force f cp among FE nodes
comprising the interface segment. . . . . . . . . . . . . . . . . . . .

80

Figure 5.5

Projection of an embedded peridynamic node on an interface segment. 83

Figure 5.6

Discretization of a one-dimensional bar.

. . . . . . . . . . . . . . .

86

Figure 5.7

Axial displacement along the bar using (a) VL-coupling scheme and
(b) CT-coupling scheme. . . . . . . . . . . . . . . . . . . . . . . . .

87

Figure 5.8

Three-dimensional bar subjected to tension. . . . . . . . . . . . . . .

89

Figure 5.9

Displacements (a) ux and (b) uy at different locations using the VLcoupling scheme. Coordinates of the measuring position are given in
parenthesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

Displacements (a) ux and (b) uy at different locations using the CTcoupling scheme. Coordinates of the measuring position are given in
parenthesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

Figure 5.10

xiii

Figure 5.11

Axial displacement along the edge using the CT-coupling scheme. .

94

Figure 5.12

Strain εx distributions (a) on the bar and (b) on the interface when
the applied traction σx = 700 MPa. The CT-coupling scheme is used
in the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

Summation of coupling forces in the longitudinal direction on the
interface elements at the left end of the bar. The CT-coupling scheme
is used in the simulation. . . . . . . . . . . . . . . . . . . . . . . . .

96

Mixed mode fracture test: (a) geometry of the specimen; (b) subregions of the coupling model. . . . . . . . . . . . . . . . . . . . . . .

98

Figure 5.13

Figure 5.14

Figure 5.15

Numerical prediction of crack paths using the VL-coupling scheme.
Boundary displacement: (a) un = 0.0220 mm; (b) un = 0.0225 mm;
(c) un = 0.03 mm; (d) un = 0.09 mm. . . . . . . . . . . . . . . . . 101

Figure 5.16

Numerical prediction of crack paths using the CT-coupling scheme.
Solid curves are crack paths observed in the experiments from [130].
Boundary displacement: (a) un = 0.0220 mm; (b) un = 0.0225 mm;
(c) un = 0.03 mm; (d) un = 0.09 mm. . . . . . . . . . . . . . . . . 102

Figure 6.1

Slave node and four master segments. . . . . . . . . . . . . . . . . . 111

Figure 6.2

Projection of the vector g onto a master segment. . . . . . . . . . . 113

Figure 6.3

Slave node projected on the intersection of two master segments. . . 113

Figure 6.4

Location of a contact point. . . . . . . . . . . . . . . . . . . . . . . . 114

Figure 6.5

Perforation of a thin plate by a blunt projectile. . . . . . . . . . . . 116

Figure 6.6

Peridynamic bar subjected to tension. . . . . . . . . . . . . . . . . . 117

Figure 6.7

(a) Constitutive model defined for peridynamic bonds; (b) macroscopic material response of the peridynamic bar (strain εx in the
range of 0.00 to 0.05); (c) macroscopic material response of the peridynamic bar; (d) comparison of peridynamic material behavior at
the macroscopic level to the experimental result of Weldox 460 E
steel (strain rate in the range of 0.00074 s−1 to 1522 s−1 ) presented
in [21]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
xiv

Figure 6.8

Strain εx contour: (a) boundary displacement ux = 0.12 mm; (b)
boundary displacement ux = 0.50 mm; (c) boundary displacement
ux = 1.00 mm; (d) boundary displacement ux = 2.50 mm. . . . . . . 122

Figure 6.9

Impact between two rigid bodies. . . . . . . . . . . . . . . . . . . . . 123

Figure 6.10

(a) Displacement of the projectile and (b) velocity of the projectile.

Figure 6.11

Modeling of a peridynamic plate impacted by a cylindrical projectile. 126

Figure 6.12

(a) Projectile displacement and (b) residual velocity. . . . . . . . . . 129

Figure 6.13

Indentation formed by the projectile with an initial velocity vi =
250 m/s: (a) bird’s-eye view of the plate; (b) lateral view of the plate. 130

Figure 6.14

Perforation of the plate by the projectile an the initial velocity vi =
300 m/s: (a) t = 30 µs, bird’s-eye view of the plate; (b) t = 30 µs,
lateral view of the plate; (c) t = 70 µs, bird’s-eye view of the plate;
(d) t = 70 µs, lateral view of the plate; (e) t = 120 µs, bird’s-eye view
of the plate; (f) t = 120 µs, lateral view of the plate. . . . . . . . . 132

Figure 6.15

Perforation of the plate by the projectile an the initial velocity vi =
400 m/s: (a) t = 30 µs, bird’s-eye view of the plate; (b) t = 30 µs,
lateral view of the plate; (c) t = 70 µs, bird’s-eye view of the plate;
(d) t = 70 µs, lateral view of the plate; (e) t = 120 µs, bird’s-eye view
of the plate; (f) t = 120 µs, lateral view of the plate. . . . . . . . . 133

Figure 7.1

Finite and heterogeneous solid [169]. . . . . . . . . . . . . . . . . . . 136

Figure 7.2

(a) Relative density versus angle of coalescence and (b) theoretical
solution of Young’s modulus of porous materials. . . . . . . . . . . 141

Figure 7.3

Spherical pores with annular flaws [99]. . . . . . . . . . . . . . . . . 142

Figure 7.4

Cubic voids in a specimen of dimensions of 25 mm×25 mm×25 mm.
Porosity P = 5% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

Figure 7.5

Cross-sectional stress versus strain of specimens containing cubic voids.146

Figure 7.6

(a) Normalized Young’s modulus and (b) normalized tensile strength
of specimens containing cubic pores. . . . . . . . . . . . . . . . . . 147
xv

125

Figure 7.7

Spherical pore distributions in the specimen and distributions at the
cross section along the longitudinal direction. Porosities (a) P = 5%,
(b) P = 10% , (c) P = 15%, and (d) P = 20%. Pores are generated
within the specimen boundaries. . . . . . . . . . . . . . . . . . . . . 151

Figure 7.8

Cross-sectional stress versus strain of specimens containing spherical
pores within the boundary. . . . . . . . . . . . . . . . . . . . . . . . 152

Figure 7.9

(a) Normalized Young’s modulus and (b) tensile strength of specimens
containing spherical pores within the boundary. . . . . . . . . . . . 153

Figure 7.10

Spherical pore distributions in the specimen. Porosities (a) P = 5%,
(b) P = 10% , (c) P = 15%, and (d) P = 20%. Pores intersect with
specimen boundaries. . . . . . . . . . . . . . . . . . . . . . . . . . . 154

Figure 7.11

Cross-sectional stress versus strain of specimens containing spherical
pores intersecting with boundaries. . . . . . . . . . . . . . . . . . . . 154

Figure 7.12

(a) Normalized Young’s modulus and (b) tensile strength of specimens
containing spherical pores intersecting with boundaries. . . . . . . . 155

Figure 8.1

State-based peridynamics [156]. . . . . . . . . . . . . . . . . . . . . 163

xvi

LIST OF SYMBOLS

x

Material point in the reference configuration

x′

Neighboring material point in the reference configuration

ρ

Density

u

Displacement field

∇

Divergence

σ

Stress matrix

b

Body force density vector

Hx

Neighborhood of the material point at x

Vx ′

Volume of neighboring material point

ξ

Relative position vector

η

Relative displacement vector

f

Bond force

δ

Horizon

∆x

Grid spacing

c

Micromodulus

s

Bond stretch

rj

Volume reduction distance

w

Micropotential
xvii

c1
˜

One-dimensional analytical micromodulus

c3
˜

Three-dimensional analytical micromodulus

c1

One-dimensional numerical micromodulus

c3

Three-dimensional numerical micromodulus

E

Young’s modulus

ν

Poisson ratio

A

Nodal cross-sectional area

NL

Set of bonds

L
fV

Summation of bond forces

Epd

Peridynamic Young’s modulus

ˆ
f

Compensation force

ˆ
Epd

Modified peridynamic Young’s modulus

ˆ
σ

Lateral stress

c3
ˆ

Modified three-dimensional numerical micromodulus

NH

I

Number of neighboring nodes

µ

Scalar function to determine bond failure

sy

Yield stretch

s0

Critical stretch

Gf

Energy release rate

k

Bulk modulus

C

Elastic constitutive matrix

ε

Strain

H

Displacement interpolation matrix

xviii

B

Strain-displacement matrix

M

Mass matrix

K

Stiffness matrix

F

Force vector

Fint

Internal force

Fext

External force

Fcp

Coupling force

φ

Shape function

Ueb

Displacement of an embedded node

ξ, η, ψ

Natural coordinates of an embedded node

ξ c , ηc

Natural coordinates of a projection point

I

Converged solution of natural coordinates of an embedded node

t

Position vector of an embedded node

r

Position vector of a projection point

fs

Peridynamic short-range force

ds

Maximum distance for none-zero short-range force

ms

Master node

ns

Slave node

g

Vector from master node to slave node

c

Vectors along segment edge

m

Unit normal vector of a segment

s

Projection of vector g

l

Penetration length

xix

kcs

Contact stiffness

m1

Nodal mass of a slave node

m2

Segment mass

Ξp

Penalty scale factor

Ξs

Scale factor of slave stiffness

Ξm

Scale factor of master stiffness

Fp

Penalty force

vi

Initial velocity

vbl

Ballistic limit velocity

vr

Residual velocity

mp

Projectile mass

mt

Plug mass

P

Porosity

X

Relative density

θ

Coalescence angle

D

Pore diameter

xx

Chapter1
Introduction
1.1

Overview of peridynamics

The analysis of problems involving discontinuities such as cracks is a long-standing challenge
in the engineering field. Fundamentals of linear elastic fracture mechanics were established
around 1960 [5]. The development of fracture mechanics provides a more reliable methodology than the traditional strength-based approach for engineering designs. However, the
primary concern of the classical fracture mechanics is related to problems with pre-existing
defects within a body, rather than spontaneous formations of discontinuities in the material.
With the advent of computers, the finite element method has been developed for solving a
wide range of engineering problems, for example, structural mechanics, heat transfer, and
fluid flows [11]. Despite the effectiveness and applicability of the finite element method in
many engineering analyses, it is difficult to use the finite element method for numerical predictions of crack growth and damage since the finite element formulations are based on the
partial differential equations in classical continuum mechanics.
Peridynamics [153] is a reformulated theory of continuum mechanics. In the peridynamic
theory, the partial differential equations that appear in the classical continuum mechanics
are replaced with integral equations. Specifically, the stress divergence term in the equation
of motion for a continuum is replaced with an integral function. Consider a body occupying
1

the region R as shown in Figure 1.1. A spatial range, which is called the horizon δ, is
associated with a material point x in the reference configuration. The material point x can
interact with all material points within the horizon δ. The interaction between the material
point x and the material points x′ is formed in terms of the bond force . In the mathematical
framework, the total interacting force on the material point x is determined by integrating all
bond forces over the spatial domain Hx , which represents the neighborhood of the material
point x within the horizon δ. In dynamic analyses, the acceleration of the material point x
is determined by the total interacting force and applied external forces.

Horizon δ

X

f

X’

y

x

z
Figure 1.1: Schematic of peridynamics. (For interpretation of the references to color in this
and all other figures, the reader is referred to the electronic version of this dissertation.)

The main advantage of peridynamics is that it is directly applicable to represent discontinuities. This is in contrast to approaches based on the mathematical framework of classical
continuum mechanics. Special techniques are needed to treat discontinuities in the classical
2

approaches since the spatial derivatives in the equation of motion are undefined along discontinuities. For example, cohesive crack models [10, 44] were introduced to take account of
nonlinearity near crack tips. In order to incorporate cohesive zone models into finite element
analyses, interface elements are inserted a priori into the finite element mesh [39]. In the
extended finite element method [13], enrichment functions, which work as an additional set
of functions for the approximate displacement field, are added to the global displacement
field for simulations of crack propagation.

In local continuum models, material points can only interact with adjacent material
points by means of contact forces [113]. In contrast, materials points separated by a finite
distance in the reference configuration can interact with each other in the peridynamic theory. Therefore, the peridynamic theory is categorized as a nonlocal method [153]. Although
a substantial amount of research on nonlocal methods has been done before the development
of peridynamics, most nonlocal modelings involve spatial derivatives [153]. Prior to the development of peridynamics, research progress has been made on mesh free particle methods
[106] such as smoothed particle hydrodynamics and meshfree Galerkin methods. The differences between peridynamics and meshfree particle methods deserve to be noticed. In the
smoothed particle hydrodynamics, for example, the acceleration of a particle is governed by
partial differential equations [124]. On the contrary, the motion of a node in peridynamics is
governed by an integral equation. In numerical implementations of peridynamics, a material
region is discretized into nodes. Since peridynamics is a continuum model in essence, it
might be considered as a continuum version of molecular dynamics [115].
3

1.2

Research Objectives

Compared with the classical theory, peridynamics is a relatively new methodology. The
thrust of research on peridynamics is providing an alternative theory to the classical continuum mechanics that is directly applicable for numerical simulations of spontaneously formed
discontinuities. In this dissertation, the primary focus is the development of the discretized
bond-based peridynamics for solid mechanics. The following research objectives are identified:
• Establish a connection between the classical elasticity and the discretized peridynamics
for one- and three-dimensional models.
• Investigate material behaviors of peridynamic models and the influence of the horizon.
• Develop a numerical scheme for peridynamic simulations of materials with Poisson’s
ratio not equal to 1/4.
• Develop an algorithm enabling optimal parallelization and implement the peridynamics
code in graphics processing units for high performance computing.
• Compare the macroscopic material behaviors of peridynamic models and finite element
analyses and propose a multiscale approach to bridge material models defined at the
different scales.
• Develop a coupling approach of discretized peridynamics and finite element method
to take advantage of the generality of peridynamics in the presence of discontinuities
and the computational efficiency of finite element method. For the validation of the
proposed coupling approach, elastic and fracture problems will be examined.
4

• Propose numerical scheme for the modeling of contact between a peridynamic domain
and a non-peridynamic domain such as conventional finite elements and rigid bodies.
Rigid-body impact and perforation of thin plate will be investigated to validate the
numerical scheme for contact-impact problems.
• Develop an algorithm to generate voids for a given porosity in peridynamic models to
study the degradations of Young’s modulus and strength of porous brittle materials.

1.3

Scope

The dissertation is organized as follows: In Chapter 1, the overview of peridynamics and the
differences compared with the classical theory are summarized. In Chapter 2, the literature
review regarding numerical predictions of crack growth and the development of peridynamics
is provided. In Chapter 3, a connection between the classical elasticity and the discretized
peridynamics is addressed in terms of peridynamic stress. Micromoduli are determined numerically for one- and three-dimensional simulations. A pairwise compensation scheme is
proposed for the simulations of materials with Poisson ratio other than 1/4. In Chapter 4,
peridynamics is employed to study brittle and ductile materials. In order to enhance the
computational efficiency, a peridynamics code is parallelized using an NVIDIA graphics processing unit in the form of a high-level implicit programming model. A multiscale procedure
is proposed to bridge the peridynamic material model at the bond scale and the material
model at the macroscopic scale for finite element analyses. In Chapter 5, an approach to
couple discretized peridynamics with finite elements is proposed. An interface element is
introduced to couple peridynamic and finite element subregions, and two types of coupling
schemes are examined. The coupling approach is employed to study elastic deformations
5

of one- and three-dimensional models and the mixed-mode fracture in a concrete specimen.
In Chapter 6, a contact-impact procedure is developed for the modeling of contact between
a peridynamic domain and a non-peridynamic domain such as conventional finite elements
and rigid bodies. The proposed procedure is applied to study the impact between two rigid
bodies and the ballistic perforation of a thin plate. In Chapter 7, porous brittle material is
studied using peridynamics. The effects of porosity on the degradation of Young’s modulus
and strength are investigated. The last chapter is the conclusion and descriptions of future
works.

6

Chapter2
Literature review
Fracture is one of major concerns in engineering field for a long time, and researchers
have made substantial efforts in order to understand material failures and alleviate potential
dangers. Inglis, Griffith, and other researchers made contributions to the early development
of fracture analyses [5]. Irwin extended the Griffith approaches by developing the energy
release rate [5]. The fracture mechanics provides a more reliable methodology for engineering design than the traditional strength based approach. However, the classical fracture
mechanics has its limitations. For example, a pre-existing crack needs to be defined, and the
fracture process zone is required to be small compared to geometrical dimensions [86, 87].

2.1

Numerical predictions of crack growth

Numerical predictions of crack initiation and growth have been considered as a class of challenging problems. Approaches based on the mathematical framework of classical continuum
mechanics fail to be directly applicable to describe discontinuities since the theory is formulated in partial differential equations, and a unique spatial derivative, however, does not exist
on the singularities. In order to model crack growth and material damage, a considerable
amount of research has been done.
Barenblatt [10] and Dugdale [44] first introduced the cohesive crack models, which address
7

relationships between the cohesive tractions resisting the separation of cracks and the crack
opening displacement [151]. The cohesive zone models are incorporated into finite element
models using, for example, interface elements and contact surfaces. Numerous researchers
have adopted cohesive zone models in numerical simulations. For example, Foulk et al. [60]
implemented a cohesive zone model in nonlinear finite element formulations. Ruiz et al. [151]
studied dynamic mixed-mode fracture using three-dimensional cohesive models. Overviews
of cohesive crack models have been presented by Elices et al. [48] and Planas et al. [141].
The partition of unity finite element method (PUFEM) was presented by Melenk and
Babuska [120]. Belytschko and collaborators [13, 123, 167] investigated the partition of
unity principle for numerical simulations of fracture problems and introduced the extended
finite element method (XFEM) to alleviate shortcomings of the conventional finite element
method. XFEM allows the discontinuity not constrained to element boundaries, and it
can model the discontinuity without remeshing [35]. Mo¨s and Belytschko [122] used the
e
extended finite element method to model the growth of arbitrary cohesive cracks. Sukumar
et al. [168, 83] implemented XFEM with Dynaflow and conducted crack growth simulations
of channel-cracking in thin films. Mariani and Perego [117] presented a method for the
simulation of quasi-static cohesive crack propagation using XFEM for quasi-brittle materials.
Cox [35] proposed enrichment functions to represent the discontinuity using an analytical
investigation of the cohesive crack problem. Considering a global level of minimizing the
total energy of the system, Meschke and Dumstorff [121] proposed a variational form of
XFEM for the propagation of cohesive and cohesionless cracks in quasi-brittle solids.
A large number of different techniques have been carried out by researchers. For example,
Ibrahimbegovic and Delaplace [85] used microscale and mesoscale discrete models to study
dynamic fracture. Kozicki et al. [98] showed a lattice based discrete approach to model
8

fracture in brittle materials. Moslemi and Khoei [125] employed an adaptive finite element
analysis to model curved crack growth. Ooi and Yang [132] presented the scaled boundary
finite element method to study crack propagation problems. Citarella and Buchholz [33]
investigated crack growth by the dual boundary element method.

Different from FEA in which elements are connected by a topological mesh, meshfree
particle methods employ a finite number of discrete particles to describe the state of a system [106]. The meshfree particle methods can be classified into microscopic, mesoscopic
and macroscopic meshfree particle methods according to the length scales [109]. The difficulties for simulating many engineering problems such as penetration and fragmentation
include remeshing and mapping the state variables from the old mesh to the new mesh [106].
Unlike finite element methods, meshfree particle methods eliminate mesh constraints, and
demonstrate advantages in many applications. Chen et al. [30] presented large deformation
analysis of nonlinear elastic and inelastic structures based on Reproducing Kernel Particle
Method (RKPM). The application of RKPM includes, for example, elastic-plastic deformation and hyperelasticity [111]. Meshfree particle methods [106] demonstrate capability for
numerical simulations of material failures. Molecular dynamics, for example, is capable of
investigating nonlinearities in the vicinity of cracks, the bond breaking between atoms, and
the formation of extended defects [150, 63, 25]. Holian et al. [78] simulated the openingmode fracture under the tensile loading using molecular dynamics. Research on coupling
finite element methods and meshless methods [108] has also been conducted. De et al. [37]
presented development of the method of finite spheres. Hong et al. [79, 80] used analytical
transformations before the numerical integration to improve the method of finite spheres and
proposed a technique to couple finite element and finite sphere discretizations.
9

2.2

Peridynamics

Compared with the classical continuum mechanics, peridynamics is a relatively new development. Researchers have studied the peridynamic theory, its convergence to the classical
theory, and numerical applications. In the following, the development of peridynamics and
some crucial articles are going to be reviewed.
The first paper on the topic of peridynamics was published by Silling [153]. The development of peridynamic theory was motivated to propose an alternative theory of continuum
mechanics that is useful for solving problems involving spontaneous formations of discontinuities without any special treatment on discontinuities. The main difference between
peridynamics and the classical theory and other nonlocal methods is that spatial derivatives
are eliminated in the peridynamic theory. The theory in [153] is specified as the bond-based
peridynamics since a pairwise force function is used to describe the interaction between two
material points.
The peridynamic formulation was used to study the deformation of an infinite elastic
bar by Silling et al. [161]. In this work, it is found that the peridynamic solution converges
to the classical solution as the horizon goes to zero. In addition, the solution shows some
special features of peridynamics that are not presented in the classical theory. For example,
the oscillation of the displacement field decays from the loading region and spreads out
to infinity. The same smoothness between the displacement field and the body force field
is found in the peridynamic theory. On the other hand, the displacement field is twoorder smoother in derivatives than the body force field in the classical theory. On the onedimensional bar problem, Bobaru et al. [20] discussed three types of numerical convergence
in peridynamics, the adaptive refinement, and scaling. Other theoretical works include
10

the dynamic responses of a peridynamic bar [183] and the well-posedness and structural
properties in the peridynamic equation [50].
A generalized formulation of bond-based peridynamics was described by Silling et al.
[156], which is referred to as the state-based peridynamics. A force state similar to the stress
tensor in the classical continuum mechanics is introduced. Different from the bond-based
peridynamics, the interaction between two material points might not be along the direction
of the deformed bond in the state-based peridynamics. The convergence of peridynamic
state to the classical elasticity was studied by Silling and Lehoucq [157]. It is shown that the
peridynamic stress tensor convergences to a Piola-Kirchhoff stress tensor as the length scale
goes to zero. Using the state-based peridynamic method, elastic deformation and fracture
of a bar were studied by Warren et al. [181], and viscoplasticity was studied by Foster et al.
[58, 59]. Silling and Lehoucq [158] presented the development of peridynamic theory of solid
mechanics.
With the general applicability of peridynamics, many applications of the theory have
been studied. The first publication on the numerical simulation using the peridynamic
model was authored by Silling and Askari [155]. In their work, the bond-based peridynamics is employed to study the numerical convergence in an opening-mode fracture problem
and the impact of a sphere on a brittle target. Dayal and Bhattacharya [36] studied the
kinetics of phase transformations using peridynamics without any additional kinetic relation
or the nucleation criterion. By adding pairwise peridynamic moments, Gerstle et al. [65]
proposed a micropolar peridynamic model. Demmie and Silling [40] reviewed the development in peridynamics and conducted simulations of extreme loading on concrete structures
using peridynamics. Kilic and Madenci [93] studied structural stability and failure analysis,
and employed peridynamics to predict crack paths in a quenched glass plate [92]. Using
11

the bond-based peridynamics, Ha and Bobaru [70, 69] reproduced dynamic fracture phenomenon observed in experiments and found the source of asymmetry in the crack path in
a perfectly symmetric computational model. Silling et al. [160] proposed a material stability condition for crack nucleation in an elastic peridynamic body. By comparing with
experimental studies, Agwai et al. [3] conducted a comparative study of the extended finite
element method, cohesive zone model, and peridynamics. Peridynamics can be implemented
within the framework of molecular dynamics [134].
Peridynamics has also been applied to study composites. By introducing fiber bonds,
matrix bonds, and interply bonds, Xu et al. [8, 188] applied peridynamics to model damage
and failure in composites. Bonds break irreversibly when the bond stretch exceeds a critical
value. Kilic et al. [91] used a random number generator to determine fiber locations in a
lamina, and studied the damage in center-cracked laminates with different fiber orientations
using peridynamics. Hu et al. [81, 82] derived the critical stretches for fiber bonds and matrix
bonds considering intralamina fracture energy for longitudinal and transverse loadings, and
proposed a homogenized peridynamic description of fiber-reinforced composites.
Compared with the finite element method, peridynamics is computationally expensive.
Macek and Silling [115] compared computer wall-clock run times between the EMU and FEA
implementations of peridynamics. The FEA implementation is much faster than the direct
meshless method. Askari et al. [8] applied peridynamics to model material and structural
failures using the Columbia Supercomputer at NASA Advanced Supercomputing Division.
In order to take advantage of the general applicability of peridynamics and the efficiency of
the finite element method, Macek and Silling [115] used embedded nodes and elements to
couple the peridynamic domain and conventional meshes. Similar technique was employed
to study shock and vibration reliability of leadfree electronics [103]. Kilic and Madenci [94]
12

proposed a coupling approach using overlapping regions in which both peridynamic and
finite element equations are utilized. From the perspective of energy equivalence, Lubineau
et al. [113] developed a morphing strategy for the coupling of non-local and local continuum
mechanics.

13

Chapter3
Discretized Peridynamics for Linear Elastic Solids
3.1

Introduction

Peridynamics is a new development of continuum mechanics that can simulate fractures
and other discontinuities [153]. It reformulates the mathematical description of continuum
mechanics in the forms of integral equations rather than partial differential equations [159].
The essence of peridynamics is that it computes forces on a material point using integral
equations, and it might be considered as a continuum version of molecular dynamics [155].
The interacting forces between material points over certain distance render the methodology
in the category of the nonlocal theory [46, 49, 183]. Unlike classical continuum mechanics,
the definitions of the stress and strain have not been established in the peridynamic theory
due to the nonlocal interactions. Although some efforts have been made to introduce the
peridynamic stress [104], the complexity of the definition itself inhibits implementation in
numerical simulations.
The bond-based peridynamics has limitations in representing continuum material properties. It can only simulate materials with Poisson ratio of 1/4 due to the property of Cauchy
crystal that only involves two-particle interactions [155]. In order to extend the capability
of peridynamics to model various materials that have different material properties, a gen14

eralized peridynamics was proposed, which directly incorporates a constitutive model from
conventional solid mechanics [156, 154, 181, 58]. Gerstle et al. [65] proposed the micropolar
peridynamics by adding rotational degrees of freedom in the links. However, straightforward
numerical schemes to model peridynamic materials of Poisson ratio other than 1/4 have been
under development.
In this chapter, micromoduli for one- and three-dimensional discretized peridynamic models are derived by equilibrating the peridynamic Young’s modulus and the conventional
Young’s modulus under constant strain condition. The peridynamic stress corresponds to
the stress in the classical (local) continuum theory, and the conventional constitutive law is
utilized to obtain the peridynamic Young’s modulus. Numerical studies are performed using
the numerically derived micromoduli for one- and three-dimensional simulations, respectively. Comparisons of strain distributions between the peridynamic solutions and classical
(local) elasticity solutions are conducted, and we discuss the boundary effect observed in the
simulations. In addition, a new pairwise compensation scheme for discretized peridynamics
is proposed to model materials of Poisson ratios other than 1/4, and the methodology is
verified with a numerical example.

3.2
3.2.1

Theory
Peridynamic formulation

The equation of motion in the classical continuum mechanics is derived from the principle
of linear momentum that the rate of change of linear momentum equals the force applied on
15

the body as [118]
ρ¨ = ∇ · σ + b,
u

(3.1)

where ρ is the density, u is the acceleration vector, σ is the stress matrix, and b is the body
¨
force vector. The formulation requires a unique spatial derivative which, however, does not
exist along discontinuities. In contrast, peridynamics uses integration to compute the force
on a material point, and the equation of motion of the material point at x in the reference
configuration at time t, as shown in Figure 1.1, is written as [135]

ρ¨ (x, t) =
u
Hx

f (η, ξ)dVx′ + b(x, t),

(3.2)

where f is a pairwise force vector that the material point at x′ exerts on the material point
at x, and Hx is a neighborhood of the material point at x.

u(x’, t)
x’

η+ξ

ξ

u( x , t) x
-f

f

δ

xJ

xI

(a)

(b)

Figure 3.1: (a) Relationships among relative position vector and the relative displacement
vector within a peridynamic horizon. (b) Pairwise force vector.

The relative position vector in the reference configuration shown in Figure 3.1 is expressed
16

as [155]
ξ = x′ − x,

(3.3)

and the relative displacement vector at time t is written as

η = u(x′ , t) − u(x, t).

(3.4)

For each material point, a scalar δ, called the horizon, is assumed to exist to determine
the interacting spatial range between the material point at x and the material point at x′
such that
f (η, ξ) = 0 ∀η,

if ξ > δ.

(3.5)

The pairwise force vector f has the direction of η + ξ, which connects the material point at
x and the material point at x′ in the deformed body, as

f (η, ξ) = f (η, ξ)

η+ξ
,
η+ξ

(3.6)

where f is a scalar-valued pairwise force function, and · is the Euclidean norm. The force
vector f has following properties [153]: The first property is

f (−η, −ξ) = −f (η, ξ)

∀η, ξ,

(3.7)

that describes the balance of linear momentum. The other property

(η + ξ) × f (η, ξ) = 0

17

(3.8)

arises from the balance of angular momentum. It should be noted that the balance of angular
momentum is satisfied by the sum of force couples which produces zero moment [156]. For
the material point at x′ within the horizon of the material point at x, the scalar-valued
pairwise force function is expressed as

f (η, ξ) = c × s(t, η, ξ),

(3.9)

where c is the micromodulus, and s is the bond stretch which possesses a similar concept of
the strain in elasticity. The bond stretch s is defined as

s(t, η, ξ) =

η+ξ − ξ
,
ξ

(3.10)

where ξ is the original bond length in the reference coordinate, and η + ξ is the current
bond length. If the stretch s = 0, then there is no pairwise force f between material points.
To evaluate the interactive forces among material points and to solve peridynamic equation of motion, the material domain is discretized with a number of nodes as shown in
Figure 3.2. The distance between two adjacent nodes is ∆x in x−, y−, and z−directions.
All nodes have the same volume (∆x)3 .
In order to consider the volume reduction of a node which has an intersection with the
horizon boundary as illustrated in Figure 3.3, a volume reduction scheme is introduced as
[135]

VJ ( ξ ) =




 δ− ξ + 1
 2r

2

j



 VJ






0


VJ

if (δ − rj ) ≤ ξ ≤ δ
if ξ ≤ (δ − rj )
otherwise

18

,

(3.11)

∆x
∆x
J
I
δ

Figure 3.2: Discretized domain for computation.

where (δ −rj ) is the distance from which the volume is reduced, and rj is chosen to be half of
the grid spacing ∆x in numerical implementations. Figure 3.4 illustrates the volume change
of the node J as the distance between the node I and the node J increases.

VJ
rj

.
.
.

δ −rj

δ

VI

Figure 3.3: Volume calculation scheme for discretized peridynamics. The volume is reduced
on the boundary of a horizon.

19

Volumetric ratio

1.5

1.0

0.5

0.0

δ-rj

0

δ

Distance
Figure 3.4: Volumetric ratio in a horizon. The ratio decreases to 1/2 at the border of a
horizon [135].

3.2.2

Micromodulus of elastic materials

A peridynamic material is considered to be microelastic if the bond force is defined as
[153, 189]
f (η, ξ) =

∂w
(η, ξ)
∂η

∀η, ξ,

(3.12)

where w is the micropotential. Setting the strain energies of peridynamics and classical
linear elasticity identical, the constant micromodulus c1 for one-dimensional peridynamics
˜
is obtained as [20]
c1 =
˜

2E
,
Aδ 2

(3.13)

where A is the cross-sectional area, and c3 for three-dimensional models is [49, 50]
˜

c3 =
˜

18k
,
πδ 4

20

(3.14)

where k is the material bulk modulus. In the derivation, the Poisson ratio ν is a fixed value
of 1/4.

3.2.2.1

One-dimensional model

Consider a bar subjected to a uniaxial tension in the longitudinal direction as illustrated in
Figure 3.5. The distance between nearest nodes is ∆x, and each node has a constant volume
(∆x)3 . The bonds can be considered as springs for elastic materials. Let L be the set of
bonds passing through or ending at the cross section AI of a node from the positive side.
The total force per unit volume acting through the cross section of the node I that has the
horizon δ = 2∆x is written as

L
fV (xI ) =

NL

f (η, ξ)VJ
(3.15)

J=1

1
1
= f VJ + f VJ + f VJ = 2c1 sVJ ,
2
2
where NL is the number of bonds in the set L. The force contribution is reduced to 50% if
the distance between two interacting nodes is equal to the horizon by applying the volume
reduction scheme in Equation (3.11).
The resultant force through the cross section of the node I is given by

F = 2c1 sVJ VI ,

(3.16)

and the peridynamic stress σx at the node I is calculated as

σx = 2c1 sVJ VI /AI = 2c1 s∆x4 .
21

(3.17)

δ=2 ∆x

AI

∆x

n

xI
∆x

Figure 3.5: Pairwise forces acting through a cross section of a node for δ = 2∆x in a onedimensional domain. Each node, represented by a sphere, has a volume of (∆x)3 .

By applying Hooke’s law, the corresponding peridynamic Young’s modulus Epd in onedimensional peridynamic models can be obtained in terms of the peridynamic stress and the
bond stretch as
Epd =

σx
= 2c1 ∆x4 .
s

(3.18)

Setting the peridynamic Young’s modulus Epd equal to the Young’s modulus E of the
isotropic material, the micromodulus c1 for one-dimensional peridynamic models is obtained
as
c1 =

E
.
2∆x4

(3.19)

The micromodulus c1 for different horizons in one-dimensional models are listed in Table 3.1. We designate c as numerical micromodulus since the number of bonds passing
through the cross-sectional area of a node needs to be calculated numerically for different
horizons. Compared with the one-dimensional micromodulus given in Equation (3.13) which
is obtained from the elastic energy density in the classical theory, it can be shown that the
one-dimensional analytical micromodulus yields the numerical micromoduli summarized in
22

Table 3.1 if j∆x (j = 2, 3, 4 and 5) is substituted to δ, and (∆x)2 is substituted to A in
Equation (3.13).
Horizon δ

Micromodulus c1 (N/m6 )

2∆x

E
2∆x4
2E
9∆x4
E
8∆x4
2E
25∆x4

3∆x
4∆x
5∆x

Table 3.1: Micromodulus c1 in one-dimensional domain.

3.2.2.2

Three-dimensional model

In the three-dimensional domain, an elastic body is subjected to the external force such that
the bond stretch is a constant value s in the domain. Let L be the set of bonds passing
through or ending at the cross section AI of a node from the positive side as shown in
Figure 3.6. The horizons of 2∆x, 3∆x, 4∆x, and 5∆x are selected since the horizon smaller
than 1∆x would lead to no peridynamic bonds within the horizon, and horizons larger than
5∆x will make the computation expensive. On the other hand, if δ = 1∆x, the peridynamic
model reduces to a truss-type system, in which each node is only connected to its nearest
neighbors. Figure 3.7 represents the number of bonds in the set L for different horizons.
The number of bonds is 11 for the horizon δ = 2∆x, and it increases to 631 for the horizon
δ = 5∆x.
The projection of the pairwise force f in x-direction is expressed as

fx = f

|xJ − xI |
,
ξ

(3.20)

where f = c3 s is the magnitude of the pairwise force between node I and node J, c3 is the
23

AI

xI

n

δ

Figure 3.6: Pairwise forces acting through a cross section of a node for δ = 2∆x in a threedimensional domain. Each node, represented by a sphere, has a volume of (∆x)3 .

micromodulus for the three-dimensional model, and |xJ − xI | is the distance between two
nodes in the x-coordinate. It should be noted that the bond force is reduced to 1/2 if it
passes through the edge of the cross section, and is reduced to 1/4 if it passes through the
corner of the cross section since the portion of the forces within the domain reduce by 50%
and 75%, respectively. Considering all nodes have the same volume (∆x)3 , the peridynamic
stress σx on the cross section of the nodes I is given by

1
σx =
AI

NL

(fx VJ )VI ,

(3.21)

J=1

where NL is the number of bonds in the set L, and the cross-sectional area AI = ∆x2 in the
24

700

Number of bonds

600
500
400
300
200
100
0
2

3

4

5

δ/ x
Figure 3.7: Number of bonds acting through a cross section of a node in a three-dimensional
domain.

uniformly discretized grid. Since the elastic body is subjected to the isotropic expansion,
the stresses in y− and z−directions are identical to the x−directional stress σx , expressed
as σy = σz = σx . By applying Hooke’s law of linear elasticity, the peridynamic Young’s
modulus Epd can be obtained as

Epd =

σx ν(σy + σz )
−
,
s
s

(3.22)

where Poisson ratio ν = 1/4. By substituting Equations (3.20) and (3.21) to Equation (3.22)
and setting the peridynamic Young’s modulus Epd in Equation (3.22) equal to the Young’s
modulus of the material E, the micromodulus c3 can be calculated. For example, the micromodulus for the horizon δ = 2∆x is

c3 = 0.302942
25

E
,
(∆x)4

(3.23)

which has dimensions of force per unit volume squared. The micromoduli c3 for horizons
ranging from 2∆x to 5∆x are summarized in Table 3.2. Figure 3.8 compares the micromoduli c3 and c3 for an elastic material which has Young’s modulus E = 70 GPa and Poisson
˜
ratio ν = 1/4. The numerically determined micromoduli c3 for discretized peridynamics are
larger than the analytical micromoduli c3 . For the horizon δ = 2∆x, numerically determined
˜
micromodulus c3 is 1.27 times larger than the micromodulus c3 . As the horizon to grid spac˜
ing ratio increases, the numerical micromodulus c3 converges to the analytical micromodulus
c3 .
˜
Horizon δ
2∆x

Micromodulus c3 (N/m6 )
0.302942 E 4
∆x

0.052385 E 4
∆x
0.017290 E 4

3∆x
4∆x

∆x

0.006819 E 4
∆x

5∆x

Table 3.2: Micromodulus c3 in three-dimensional domain.

c
~3
c3

6

Micromodulus (N/m )

4e+23

3e+23

2e+23

1e+23

0e+00
2

3

4

5

δ/ x
Figure 3.8: Micromodulus for different horizons (E = 70 GPa, and ν = 1/4).

26

3.2.3

Pairwise compensation scheme

For an elastic body which has Poisson ratio not equal to 1/4, we introduce a compensation
force ˆ for each node subjected to the pairwise force f due to the stretch of peridynamic bonds.
f
A set of compensation forces ˆ is superposed in the transverse directions perpendicular to
f
the direction of the force f , resulting in additional transverse deformations which simulate
the effect of Poisson ratios other than 1/4. It should be noted that two pairs of compensation
forces in y− and z−directions are added to nodes in three-dimensional domains.
If a node is enclosed completely in an elastic body, the added compensation forces ˆ are
f
balanced, and the resultant nodal force becomes zero. On the other hand, there exists a
force added to the nodes located at the boundary. Therefore, the direction of the resultant
force is towards inside of the body if Poisson ratio is larger than 1/4, as shown in Figure 3.9.

x

f
−f

^
f

x’
^
f
x

f

−f

x’

^
−f

Figure 3.9: Pairwise compensation scheme.

Consider a three-dimensional rectangular bar, which has Young’s modulus E and Poisson
ratio ν. The bar is subjected to a uniaxial tension as shown in Figure 3.10. The compensation
27

forces acting on the lateral surfaces are superposed on the nodes at boundaries as shown in
ˆ
Figure 3.11. The peridynamic Young’s modulus Epd needs to be determined such that
the strains in Figure 3.11 are identical to the strains in the configuration of Figure 3.10.
ˆ
For Poisson ratio larger than 1/4, the modified peridynamic Young’s modulus Epd should
be larger than Epd to keep the strain identical. The longitudinal strain in x−direction in
Figure 3.11 equals the strain in Figure 3.10, which is written as

σx
1
σx
−
=
ˆ
Epd
Epd 4

σy
ˆ
σz
ˆ
+
ˆ
ˆ
Epd Epd

,

(3.24)

where σy and σz are superposed stresses from the compensation forces for ν = 1/4. The
ˆ
ˆ
strains in y− and z−directions also need to be identical to the strains in the original peridynamic model shown in Figure 3.10. Therefore, the equilibrium of the strain εy is expressed
as
−ν

σy
ˆ
σx
1
−
=
ˆ
Epd
Epd 4

σz
ˆ
σx
+
ˆ
ˆ
Epd Epd

.

(3.25)

Since the compensation forces are identical in all lateral directions perpendicular to the
pairwise force f , the lateral stresses are equal as

σy = σz .
ˆ
ˆ

(3.26)

Solving Equations (3.24) to (3.26), the transverse lateral stresses in y− and z−directions are
obtained as
σy =
ˆ

1 − 4ν
σx ,
3 − 2ν

σz =
ˆ

1 − 4ν
σx .
3 − 2ν

(3.27)

Substituting Equation (3.27) to Equation (3.24), the modified peridynamic Young’s modulus

28

can be derived as
ˆ
Epd =

5
E .
6 − 4ν pd

(3.28)

ˆ
Hence, the micromodulus c3 corresponding to Epd should be modified as
ˆ

c3 =
ˆ

5
c ,
6 − 4ν 3

(3.29)

where c3 is the numerically determined micromodulus summarized in Table 3.2.

y
x

E

σx

pd

, ν

σx

z
Figure 3.10: Three-dimensional bar subjected to a tension.

^
σy

y
x

^

E

σx

pd

, ν=1/4

σx

z

^
σz
Figure 3.11: Three-dimensional bar subjected to a tension. Lateral forces are added to take
1
account of ν = 4 .

29

3.3
3.3.1

Case studies
Numerical implementation

The material region is discretized isotropically in all the directions in space such that the
distance between two adjacent nodes is ∆x in x−, y−, and z−directions as shown in Figure 3.2. All nodes have the same volume (∆x)3 . The equation of motion in Equation (3.2)
is rewritten as
ρ¨ t =
uI

NH

I

f (η t , ξ)VJ + bt ,
I

(3.30)

J=1

where ut is the acceleration of the node I at time t, f (η t , ξ) is the pairwise force, NH is
¨I
I

the total number of nodes within the horizon of the node I, and bt is the body force at
I
time t. The velocity-Verlet scheme [51] is used for the time integration to update positions,
velocities and accelerations as

1
x(t + ∆t) = x(t) + v(t)∆t + a(t)∆t2 ,
2
1
v(t + ∆t/2) = v(t) + a(t)∆t,
2
1
a(t) = fV (t),
ρ
1
v(t + ∆t) = v(t + ∆t/2) + a(t + ∆t)∆t,
2

(3.31)
(3.32)
(3.33)
(3.34)

where x is the position vector, v is the velocity vector, a is the acceleration vector, and fV
is the scaled force vector which is defined as

fV =

NH

I

J=1

30

f (t)VJ .

(3.35)

The acceleration is calculated by dividing the scaled pairwise force fV , which has the unit
of force per volume, by the density ρ.
For the numerical simulation, peridynamics is implemented in the framework of LAMMPS
[135]. The research code is compiled and run in the parallel mode in a workstation which
has dual quad-core 64 bit CPU’s at 2 GHz, 16 GB memory and a 2 TB hard disk drive.

3.3.2

1
Comparison of peridynamics and analytical solution(ν = 4 )

3.3.2.1

One-dimensional bar

Consider a one-dimensional bar subjected to a uniaxial tension. The length of the bar is
10 mm, and the size of the horizon is set to 1.0 mm. The grid spacing ∆x is 0.5 mm, which
corresponds to δ/∆x = 2. Young’s modulus E is 70 GPa, and the density ρ = 2700 kg/m3 .
The magnitude of nodal force applied to the nodes at both ends is 87.5 N, and the force is
gradually increased for 1000 time steps, setting the time step dt = 5×10−8 sec. The classical
(local) stress σ = 87.5 N/(0.0005 m)2 = 350 MPa, and the classical (local) strain ε = 0.005.
Substituting Young’s modulus and the grid spacing to Equation (3.19), the micromodulus is
calculated as c1 = 5.6 × 1023 N/m6 .
Figure 3.12 shows the strain distribution from the middle (x=5.0 mm) to the right end of
the bar (x=10.0 mm). The strain in the middle of the bar (x=5.0 mm) is 0.00498, and the
strain increases to 0.00516 at the end of the bar. The strain distribution is very close to the
classical (local) elasticity solution. However, the boundary effect exists near the end of the
one-dimensional bar. The total force per unit volume acting through the cross section of the
L
bar fV = 6.98 × 1011 N/m3 . The resultant force obtained by multiplying the volume VI is

87.2 N. The peridynamic stress on the cross section is calculated as σ = 348.8 MPa. Dividing
31

the stress by the measured strain, the resultant Young’s modulus through back-calculation
is 70.043 GPa. The peridynamic results show good agreement in the strain and the resultant
Young’s modulus.

0.010

εx
analytical solution

0.008

εx

0.006
0.004
0.002
0.000
5

6

7

8

9

10

x (mm)
Figure 3.12: Strain distribution of the one-dimensional bar in x-direction.

3.3.2.2

Three-dimensional bar

We examine material responses of a three-dimensional rectangular bar of dimensions 10 mm
×7 mm ×7 mm with horizons of 1.0 mm, 1.5 mm, 2.0 mm and 2.5 mm. The bar is discretized
with nodes distributed uniformly in a fixed 21 × 15 × 15 grid as shown in Figure 3.13. The
grid spacing ∆x is 0.5 mm, which leads the horizon to grid spacing ratios δ/∆x to be 2, 3, 4,
and 5. Young’s modulus E of the elastic material is 70 GPa, Poisson ratio ν = 0.25, and the
density ρ = 2700 kg/m3 . Nodes at both ends of the bar are subjected to a tensile loading
gradually increasing up to 87.5 N. The corresponding classical (local) stress σ = 350 MPa,
and the analytical (local) strain εx = 0.005.
Figures 3.14(a) and (b) show the comparison of longitudinal strain εx along the center
32

Figure 3.13: Three-dimensional bar discretized by 21 × 15 × 15 nodes.
line of the bar using c3 and c3 in the numerical studies, respectively. With c3 , as shown in
˜
Figure 3.14(a), the longitudinal strain for δ = 2∆x is close to the classical (local) solution
εx = 0.005. The strain εx = 0.00486 in the middle of the bar (x = 5.0 mm), and the
strain εx = 0.00492 at the right end (x = 10.0 mm). As shown in Figure 3.14(b), the
strain obtained using the analytical micromodulus c3 in the simulation is larger than the
˜
classical elasticity solution, and the increase of the horizon lowers the strain in all the domain
except the boundary. Therefore, the corresponding resultant Young’s modulus through backcalculation is smaller than the Young’s modulus E of the material. The strain increases in
the three layers of the nodes close to the ends.

33

0.010
0.008

εx

0.006
0.004
23

6

δ=2∆x, c3=3.3929×1022 (N/m6)
δ=3∆x, c3=5.8671×1022 (N/m6)
δ=4∆x, c3=1.9365×1021 (N/m6)
δ=5∆x, c3=7.6374×10 (N/m )

0.002
0.000
5

6

7

8

9

10

x (mm)
(a)
0.010
0.008

εx

0.006
0.004
23

6

δ=2∆x, ~3=2.6738×1022 (N/m6)
c
δ=3∆x, ~3=5.2816×10 (N/m )
c
22
6
δ=4∆x, ~3=1.6711×1021 (N/m6)
c
~ =6.8449×10 (N/m )
δ=5∆x, c3

0.002
0.000
5

6

7

8

9

10

x (mm)
(b)
Figure 3.14: Strain distribution along the x-axis in the three-dimensional bar. (a) Numerical
micromodulus c3 given in Table 3.2 is used for each horizon and (b) analytical micromodulus
c3 given in Equation (3.14) is used. The classical (local) elasticity solution of strain εx =
˜
0.005.

34

The strain and stress distributions on the cross section at x = 5.0 mm are shown in
Figures 3.15(a) and (b), respectively for the horizon δ = 2∆x. The numerical micromodulus
is utilized in the simulation. The strain is 0.00486 at the center of the cross section at x =
5.0 mm, and the strain is 0.00494 at the corner. The average of the strain over the cross
section at x = 5.0 mm is 0.00491. The total force per volume acting through the cross
L
section of the bar fV = 1.53 × 1014 N/m3 , and the peridynamic stress σx = 339.9 MPa.

The resultant Young’s modulus calculated by dividing the peridynamic stress by the strain
is 69.27 GPa in the middle of the bar (x = 5.0 mm). For the horizon δ = 2∆x, the average
strain at the end of the bar (x = 10.0 mm) is 0.00503, and the resultant Young’s modulus
through back-calculation is 69.58 GPa. In y−direction, the strain εy in the middle of the
bar is -0.00121, and the resultant Poisson ratio is ν =0.249. Figure 3.16 shows Poisson ratio
obtained from the longitudinal and transverse strains from the middle to the right end of

16
12
8
4
7
56
0
34
0 1 2
12 z (mm)
3 4 5
6 7 0
y (mm)

2

-3

εx (×10 )

σx (×10 MPa)

the bar.

(a)

4
3
2
1
7
56
0
34
0 1 2
12 z (mm)
3 4 5
6 7 0
y (mm)
(b)

Figure 3.15: (a) Strain- and (b) stress-distributions on the cross section at x = 5.0 mm.
The classical (local) elasticity solution of strain εx = 0.005 and stress σx = 350 MPa. The
horizon δ = 2∆x, and the numerical micromodulus c3 is utilized in the simulation.

Figure 3.17 shows the strain distributions at the end of the bar (x = 10.0 mm) for
35

0.500
0.400

ν

0.300
0.200
0.100
0.000
5

6

7

8

9

10

x (mm)
Figure 3.16: Poisson ratio along the three-dimensional peridynamic bar as shown in Figure 3.13. The given Poisson ratio of the material is ν = 1 . The horizon δ = 2∆x, and the
4
numerical micromodulus c3 is utilized in the simulation.

different horizons using the numerical micromodulus summarized in Table 3.2. As the horizon
increases, the strain also increases at the end of the bar. For the horizon δ = 2∆x, the strain
at the center of the cross section at the end is 0.00492, the strain is 0.00540 at the corner,
and the average strain over the cross section is 0.00503 as shown in Figure 3.17(a). For the
horizon δ = 5∆x, the strain at the center of the cross section at x = 10.0 mm is 0.00732,
and the strain increases significantly to 0.01580 at corners as shown in Figure 3.17(d). The
average strain over the cross section at the end of the bar for the horizon δ = 5∆x is 0.00893,
which is much larger than the classical (local) solution εx = 0.005.
Compared with the classical (local) elasticity solution εx = 0.005, the boundary effect
increases as the horizon increases since the larger horizon causes larger loss of bond contribution near the ends of the bar where tractions are applied and near the lateral surfaces.
Consequently, the bond stretch s between those nodes close to the boundaries gets larger
due to the smaller fraction of the contributing bonds among the total number of bonds.
36

-3

εx (×10 )

-3

εx (×10 )

16
12
8
4
67
0
45
3
0 1 2
12 z (mm)
3 4 5
6 7 0
y (mm)

16
12
8
4
67
0
45
3
0 1 2
12 z (mm)
3 4 5
6 7 0
y (mm)

-3

-3

εx (×10 )

(b)

εx (×10 )

(a)

16
12
8
4
67
0
45
0 1 2
23 z (mm)
3 4 5
01
6 7
y (mm)
(c)

16
12
8
4
67
0
45
0 1 2
23 z (mm)
3 4 5
01
6 7
y (mm)
(d)

Figure 3.17: Strain distribution at the end surface of the bar for different horizons using the
numerical micromodulus: (a) δ = 2∆x, (b) δ = 3∆x, (c) δ = 4∆x, and (d) δ = 5∆x. The
classical (local) elasticity solution of strain εx = 0.005.

It has been observed that the peridynamic solution converges to the classical (local)
elasticity solution as the horizon decreases in one-dimensional models [161, 20]. Similarly,
for the selected range of the horizons (2∆ to 5∆) of three-dimensional models, material
responses, using the horizon δ = 2∆x and the corresponding numerical micromodulus c3 ,
are very close to classical (local) solutions. However, the use of the micromodulus c3 = 18k
˜
4
πδ

yields larger responses in strains since the analytical micromodulus c3 is smaller than c3 ,
˜
summarized in Table 3.2.
37

It should be noted that a different horizon may be chosen for different types of problems
such as large-strain plastic response and fracture phenomena. For example, in crack branching problems, the horizon δ = 4∆x is proven to be effective [70, 69] since adequate number
of nodes inside the horizon are required for the development of crack path. On the other
hand, larger error from the use of larger horizons at the vicinity of the boundaries might
be alleviated by normalizing forces or adding correction factors as proposed by researchers
[115, 90].

38

3.3.3

Comparison of modified peridynamics and analytical solu1
tion (ν = 4 )

We introduce the concept of compensation force ˆ for each pairwise force f to extend the
f
1
discretized bond-based peridynamics to simulate materials that have Poisson ratios ν = 4 .

Consider a rectangular bar of dimensions 10 mm×7 mm×7 mm and choose the horizon
δ = 1.0 mm. The bar is discretized with the grid spacing ∆x = 0.5 mm so that the
peridynamic solutions might be close to the classical (local) elasticity solutions. Each node
on the ends of the bar is subjected to a tension of 87.5 N in the longitudinal direction. The
density ρ = 2700 kg/m3 , Young’s modulus E = 70 GPa, and Poisson ratio ν = 0.3. The
applied traction on the ends is σ = 350 MPa, and the corresponding classical, local, elastic
strain εx is 0.005.
For each pairwise bond force on a node, two pairs of compensation forces are applied in
the perpendicular directions of the force vector f as shown in Figure 3.9. The compensation
forces are self-balanced for the nodes inside the body. Therefore, the resultant compensation
forces vanish. However, for the nodes on the boundary, the forces are superposed to the
surface normal direction.
Additional tractions are applied on all the lateral surfaces to simulate Poisson effect
as shown in Figure 3.11. The applied additional traction on lateral surfaces is given by
Equation (3.27). Substituting the target Poisson ratio ν = 0.3 and traction σx = 350 MPa
to Equation (3.27), we obtain σy = -29.2 MPa. After multiplying the cross-sectional area
ˆ
(∆x)2 , the forces applied on each node on the lateral boundary surfaces are −7.289 N in y−
and z−directions.
Setting δ = 2∆x and substituting c3 = 0.302942 E/∆x4 given in Table 3.2 to Equa39

tion (3.29), the micromodulus is modified to c3 = 3.5343 × 1023 N/m6 . Figure 3.18 shows
ˆ
the strain and the resultant Poisson ratio calculated from the measured longitudinal and
transverse strains along the center line of the bar. The strain εx at at the center of the cross
section at x = 5.0 mm is 0.00502, and the average strain over the cross section is 0.00509.
The sum of forces acting through the cross section of the bar is 19842 N, and the corresponding peridynamic stress σx on the cross section of the bar is calculated by dividing the
sum of forces by the cross-sectional area of the bar. The corresponding peridynamic stress
is 352.7 MPa, and the resultant Young’s modulus is 70.26 GPa through back-calculation.
Figure 3.18(b) shows the back-calculated Poisson ratio along the bar. Compared with
the Poisson ratio of the material ν = 0.3, the numerical results show good agreement. The
boundary effect is also observed in the distribution of Poisson ratio. Poisson ratio ν decreases
close to the end of the bar since the longitudinal strain increases at locations close to the
end of the bar.

40

0.010
0.008

εx

0.006
0.004
0.002
0.000
5

6

7

8

9

10

8

9

10

x (mm)
(a)
0.500
0.400

ν

0.300
0.200
0.100
0.000
5

6

7

x (mm)
(b)
Figure 3.18: (a) Strain εx and (b) Poisson ratio ν along the bar by the modified peridynamic
formulation. The classical (local) elasticity solution of strain εx = 0.005, and the Poisson ratio ν of the material is 0.3. The horizon δ = 2∆x, and the modified numerical micromodulus
c3 for pairwise compensation scheme is utilized in the simulation.
ˆ

41

3.4

Summary

The one-dimensional micromodulus c1 and the three-dimensional micromodulus c3 are derived for discretized peridynamic models. The micromoduli are calculated by introducing the
peridynamic stress corresponding to the classical (local) stress, calculating the corresponding
peridynamic Young’s modulus, and equilibrating the peridynamic Young’s modulus with the
conventional Young’s modulus. For three-dimensional models, the numerically determined
micromodulus c3 is larger than the analytically obtained micromodulus c3 . The micromoduli
˜
c3 and c3 converge as the horizon to grid spacing ratio δ/∆x increases.
˜
Material responses of peridynamic bars with numerical micromoduli and conventional
analytical micromoduli are investigated. Compared with classical (local) elasticity solutions,
the peridynamic results of strain and the resultant Young’s modulus by back-calculation
show good agreement for the horizon δ = 2∆x using the numerical micromodulus. The
corresponding peridynamic solutions of strain using the analytical micromodulus are larger
than the classical (local) elasticity solutions. The errors due to the boundary effect are
almost negligible for the horizon δ = 2∆x, but increase in the results of strain as the horizon
δ increases. Therefore, the horizon δ = 2∆x can be recommended in the linear elastic regime
although larger horizons might be chosen for fracture problems [70, 69]. In the numerical
implementation, a horizon might be determined first, and then the grid spacing can be
selected considering the range of material behaviors.
To model materials that have Poisson ratios other than 1/4, a force compensation scheme
adding additional forces on each pairwise force between nodes is proposed for discretized peridynamics. The numerical results match classical (local) solutions in the material responses
including the stress, the strain, back-calculated Young’s modulus, and Poisson ratio.
42

Chapter4
Discretized Peridynamics for Brittle and
Ductile Solids
4.1

Introduction

Compared with macroscopic modeling such as FEA, three-dimensional peridynamic analyses
of large simulations are computationally intensive due to the nonlocal property of peridynamics. Macek et al. [115] compared computer wall-clock run times between the EMU
and FEA implementations of peridynamics. The FEA implementation is much faster than
the direct meshless method. Askari et al. [8] applied peridynamics to model material and
structural failures using the Columbia Supercomputer at NASA Advanced Supercomputing
Division.
In this chapter, a peridynamics code for three-dimensional simulations is implemented
using an NVIDIA Tesla C1060 GPU for parallel computing. A significant speedup over the
serial calculation in a CPU is achieved, varying with the number of discretization nodes in
peridynamic models. The effects of different horizons and grid sizes are investigated for a
brittle material. The stresses and strains of FEA and peridynamic solutions, respectively,
are compared for a ductile material. A multiscale procedure is proposed to bridge plastic
simulations at different scales. The results after applying the proposed procedure show good
agreements between FEA and peridynamics.
43

4.2

Theory

The peridynamic formulations are introduced in Section 3.2, where elastic peridynamic models are examined. In this chapter, brittle and ductile materials are going to be investigated.
Figure 4.1(a) shows a brittle material model defined for peridynamics. In the elastic regime,
the bond force is a scalar-function of the bond stretch s. In order to simulate damage in
peridynamic models, a critical stretch for bond failure is introduced as s0 . Once a bond
fails, it can not sustain force any more, which makes the peridynamic model hysteretic.
Figure 4.1(b) defines a ductile material model for peridynamics. A bond yields at the yield

Bond force

Bond force

limit sy , and shows perfect plasticity upto s0 .

S0

Bond stretch

S0 Bond stretch

Sy

(a)

(b)

Figure 4.1: (a) Bond force as a function of bond stretch in a brittle material model. (b)
Bond force as a function of bond stretch in a ductile material model.

For microbrittle materials, the bond force, which is a scalar function of the bond stretch
s, is expressed as
f (η, ξ) = c × s(t, η, ξ) × µ(t, η, ξ),

(4.1)

where c is the micromodulus summarized in Tables 3.1 and 3.2 for one- and three-dimensional
models respectively, and µ is the scalar function to determine the bond failure. The scalar
44

function µ is related with the critical bond stretch as

µ(t, η, ξ) =



 1 if s(t′ , η, ξ) < s0

for all 0

t′

t,

(4.2)


 0 otherwise.

The critical bond stretch s0 for microbrittle material is obtained by setting the work required
to break all the bonds per unit fracture area identical to the energy release rate Gf [155]:
5Gf
.
9kδ

s0 =

(4.3)

By considering broken bonds, damage dependencies can be introduced into the critical bond
stretch [155, 69].

4.3
4.3.1

Computation
Numerical implementation

To calculate the bond forces on peridynamic nodes after discretization, a neighbor list is
created for each node in the reference configuration. In the subsequent calculations, the
neighbor list is referred to determine the total bond forces acting on a node. The velocityVerlet scheme [51] is used for the time integration to update displacements and velocities.
The procedures of the velocity-Verlet algorithm are explained in Section 3.3.1.

4.3.2

Implementing Peridynamics in GPU

GPUs, originally developed for the acceleration for graphics rendering, have a large number of
transistors devoted to data processing rather than data caching and flow control [131]. This
45

advantage has been taken for numerical computations. GPUs have evolved for many dataparallel and computing-intensive programs, and a variety of compilers have been developed.
Researchers have used GPUs for highly parallel computation. Harris [73] described a method
to simulate stable fluid using a GPU. Kruger et al. [101] introduced a framework for solving
sets of algebraic equations in GPUs. Manavski et al. [116] used GPUs as an accelerator
for Smith-Waterman sequence alignment. Phillips et al. [138] implemented a multi-block
turbulent flow solver in GPU processors.
We use a GPU, NVIDIA Tesla C1060, that has 4 GB memory and 30 multiprocessors,
and each multiprocessor has 8 cores. Therefore, 240 cores are available in total. Each core
can execute a sequential thread, and the clock rate is 1,296 MHz. Figure 4.2 shows the
block diagram of a GPU accelerator [186]. The GPU has its own memory which is called the
device memory, and the accelerator communicates with CPU using IO commands and DMA
memory transfer [186]. The data bandwidth is 512-bits wide. To port the peridynamics code
to a GPU platform, PGI Fortran Accelerator Programming Model [143, 142], a high-level
implicit programming model for the general purpose computation on GPUs, is utilized. A
set of compiler directives designed for GPUs is available to specify regions of the code to be
offloaded from the CPU to the GPU. The code implemented with accelerator directives is
portable even without a GPU unit since the directives are considered as comments by other
Fortran compilers.
A benchmark problem for multiplication of a matrix An×n and a vector bn×1 is written
with directives as:
!$acc region
do k=1,1000,1
do

i =1,n
46

Thread Execution Control Unit

1

2

n

3

X86 Host

...

Thread Processors
Special Function Unit
Local Memory

Host Memory
DMA
Device Memory

Figure 4.2: NVIDIA Tesla Block Diagram [186].

r(i)=0
do j = 1,n
r(i)=r(i)+A(i,j)*b(j)
enddo
enddo
enddo
!$acc end region

In the fixed format, PGI Accelerator directives are specified with the sentinel !$acc. The
directive !$acc region defines the region to be compiled for execution in the GPU [142].
For benchmarking, the multiplication is repeated one thousand times, and the summation
is stored in a vector r. Figure 4.3 compares the wall-clock times by serial calculation in
the CPU and by the parallelized calculation in the GPU. For the matrix of size n smaller
than 5,000, the GPU calculation is merely a few times faster than the serial calculation.
47

However, as the size of the matrix increases, the speedup of the GPU calculation over the
serial calculation increases significantly.

3

Wall-clock time (× 10 sec)

12

CPU
GPU

8

4

0
0

5000

10000
Matrix Size (n)

15000

20000

Figure 4.3: A benchmark problem of matrix multiplication.

For the effective parallelization in GPUs, all the necessary information in the main memory is required to be transferred to the GPU memory for the calculation. We use the data
clauses to control the transfer of selected data between the CPU and the GPU, which reduces
the system overhead in copying data. The peridynamic model is initialized in the CPU, and
then variables and arrays are copied to the GPU for the calculation within the data region
using directives as
!$acc data region copyin(list), local(list)
!$acc updateout(list)
!$acc end data region
The clause copyin(list) is used to copy the listed variables and arrays from the CPU to the
GPU memory, and the clause local(list) declares local variables and arrays to be allocated
48

only in the GPU memory [142]. Calculation results in the GPU are copied to the main
memory once per multiple calculation steps using the directive !$acc updateout(list).
It should be noted that copying the whole array is faster than copying the noncontiguous
subarrays [186].
In the subroutine to calculate bond forces on each node in a discretized peridynamic
model, the conventional algorithm includes two dependent loops as shown in Table 4.1,
which is referred to as Algorithm 1. However, the loop dependency disables optimal parallelization on GPUs. Thus the nested loops are separated into two independent loops as shown
in Table 4.2, which is referred to as Algorithm 2. Peridynamics codes using Algorithm 1 and
Algorithm 2, respectively, are compiled with PGFORTRAN v10.3 for the serial calculation
in the CPU. We compared the wall-clock time to calculate bond forces in a three-dimensional
peridynamic model, which has 4,725 nodes with the horizon δ = 2∆x. As summarized in
Table 4.3, it takes 13,113 milliseconds to complete 1 time-step using Algorithm 1 for bond
force calculations, and Algorithm 2 spends 21,080 milliseconds. By explicitly setting the
memory considering the maximum number of neighbors in the peridynamic model, Algorithm 2 consumes more memory and time than Algorithm 1 if no GPU is utilized. Overall,
Algorithm 1 is 1.61 times faster than Algorithm 2 in the CPU.
do i = 1, all nodes
M = number of neighbors of the node i
{Calculate bond forces on the node i}
do j = 1, M
f(i) = f(i) + c*s
end do
end do
Table 4.1: Dependent loops (Algorithm 1).
For the comparison of the performance by the GPU enabled calculation using Algorithm
49

N = all nodes
M = maximum number of neighbors of all nodes
{Calculate each bond force}
do i = 1, N*M
f_bond(i)=c*s
end do
{Sort bond forces to each node}
do j = 1, all nodes
f(j) = summation of f_bond sorted on the node j
end do
Table 4.2: Independent loops (Algorithm 2).
2 and the serial calculation executed in the CPU using Algorithm 1, three-dimensional peridynamic simulations are conducted, and the wall-clock times are summarized in Table 4.4.
With the grid size ∆x = 0.5 mm, the first peridynamic model consists of 11 × 8 × 8 nodes,
the second model consists of 21 × 15 × 15 nodes, and the third model consists of 41 × 29 × 29
nodes. Young’s modulus E is 70 GPa, Poisson ratio is 0.25, and the critical bond stretch is
s0 = 0.004. The time-step is set to dt = 1×10−8 s. Each simulation is carried out for 500,000
time-steps. Figure 4.4(a) shows the wall-clock time to run the peridynamic models. For the
first peridynamic model which has 704 nodes, the GPU calculation takes 3.5 minutes while
the serial calculation spends 19.5 minutes in the CPU. The GPU calculation is 5.6 times
faster than the serial calculation of the peridynamics code. For the second peridynamic
model which has 4,725 nodes, the GPU calculation spends 10.5 minutes while the serial calculation spends 144.0 minutes. For the third peridynamic model which has 34,481 nodes, the
GPU calculation spends 65.7 minutes while the serial calculation uses 1294.1 minutes. The
GPU calculation is 19.7 times faster than the serial calculation. In general, as the number
of nodes in a peridynamic model increases, the speedup of the GPU enabled version of the
peridynamics code over the serial version in the CPU increases as shown in Figure 4.4(b).
50

Algorithm 4.1

Algorithm 4.2

13,113

21,080

Wall-clock time (millisecond)

Table 4.3: The wall-clock time to calculate bond forces for one time-step using a node in the
CPU. The three-dimensional peridynamic model has 4,725 nodes, the horizon δ = 2∆x and
∆x = 0.5 mm.

Number of nodes Time (minute) by the GPU Time (minute) by the CPU
11×8×8
21×15×15
41×29×29

3.5
10.5
65.7

19.5
144.0
1294.1

Table 4.4: The wall-clock time to run peridynamic models for 500,000 time-steps on GPU
and CPU, respectively. Peridynamic models have the horizon δ = 2∆x and ∆x = 0.5 mm.

51

Wall-clock time (min)

1500

CPU
GPU

1000

500

0
0

10000

20000
30000
Number of nodes

40000

(a)
25

Speedup

20
15
10
5
0
0

10000

20000
30000
Number of nodes

40000

(b)
Figure 4.4: Comparison between the GPU enabled version of peridynamics code and the
serial version. (a) Wall-clock time and (b) the speedup of the GPU version over the serial
version.

52

4.4
4.4.1

Case studies
Brittle material

In the peridynamic framework, a three-dimensional rectangular bar is modeled as shown in
Figure 4.5(a). The dimensions of the peridynamic bar are 10 mm × 7 mm × 7 mm. Boundary
conditions are set by imposing a constant velocity to a region of nodes on both ends of the
bar such that vx = −10 mm/s at the left-end boundary region and vx = 10 mm/s at the
right-end boundary region. For the material, Young’s modulus E is 70 GPa, and Poisson
ratio ν is 0.25. Time-step dt is set to dt = 5 × 10−8 s. To simulate brittle failure of materials,
the critical stretch s0 for bond failure is set as 0.01.
To study the response of the peridynamic bar for different horizons, the rectangular bar
is uniformly discretized with nodes distributed in a 21 × 15 × 15 grid. The grid size ∆x
is 0.5 mm, and the horizons δ of 2∆x, 3∆x and 4∆x are used in the model. In order to
calculate the cross-sectional stress, the summation of all interaction forces between peridynamic nodes passing through the cross section of the bar is projected on the x-axis. Then
the cross-sectional stress is obtained by dividing the resultant force by the sectional area.
The macroscopic engineering strain is measured by dividing the variation of the length by
the original length. Figure 4.6(a) shows the stress σx versus the strain εx on the cross section at x = 5.0 mm for different horizons. By calculating the tangent of the stress-strain
curves for the horizon δ = 2∆x in Figure 4.6(a), the back-calculated Young’s modulus E is
70.46 GPa. Compared with the exact value E = 70.00 GPa of the material, the error by
the back-calculated Young’s modulus is 0.66%. For the horizon δ = 3∆x and δ = 4∆x, the
back-calculated Young’s moduli are smaller than the exact value as shown in Figure 4.6(a).
53

boundary region

boundary region
7 mm

y

x

m

m

vx

vx

7

z

3 mm

4 mm

3 mm

ux (mm)

0.05

(a)

Time (s)

0.005

(b)
Figure 4.5: (a) The geometry of a three-dimensional bar and the boundary conditions; (b)
Displacement of the boundary region versus time.

To investigate the effect of grid size, the peridynamic bar is discretized in a 21 × 15 × 15
grid with a grid size ∆x = 0.50 mm and in a 41×29×29 grid with a grid size ∆x = 0.25 mm,
respectively. Figure 4.6(b) compares the cross-sectional stress σx at x = 5.0 mm for the grid
size ∆x = 0.25 mm and ∆x = 0.50 mm, where the horizon δ is set as 2∆x. As evident
in Figure 4.6(b), the back-calculated Young’s modulus E for the grid size ∆x = 0.25 mm
is very close to the back-calculated Young’s modulus E for the grid size ∆x = 0.50 mm.
Since stresses and strains are not intrinsic variables in the peridynamic theory, the average
bond stretch within the grid width ∆x is calculated for each node as nodal strain. Strain
distributions at times t = 1.0 × 10−3 s, t = 1.6 × 10−3 s and t = 1.8 × 10−3 s denoted by
54

I (no bond failure), II (a few bond failure) and III (massive bond failure) in Figure 4.6(b),
respectively, are plotted. As the strain contours demonstrate, strain εx is about 0.005 at
time t = 1.0 × 10−3 s as shown in Figure 4.7. At time t = 1.6 × 10−3 s, the cross-sectional
stress σx reaches its maximum value, and the strain εx near the boundary regions is larger
than 0.01 as shown in Figure 4.8(b). On the cross section at z = 3.5 mm, those nodes
having the nodal strain εx in the range of 0.007 to 0.010 form diagonal patterns as shown
in Figure 4.8(d). At time t = 1.8 × 10−3 s when the cross-sectional stress σx significantly
decreases due to massive bond failure, nodes having the nodal strain εx over 0.010 form an
elliptical distribution on the cross section at z = 3.5 mm as shown in Figure 4.9(b) and
Figure 4.10(b). On the cross section at x = 5.0 mm, interior nodes have the nodal strain
εx over 0.01, and a failure core is formed on the cross section as shown in Figure 4.9(c) and
Figure 4.10(c).

55

700

δ=2∆x
δ=3∆x
δ=4∆x

600

σx (MPa)

500
400
300
200
100
0
-100
0.000

0.005

0.010

0.015

0.020

0.025

εx

(a)
700

∆x=0.25 mm
∆x=0.50 mm

600

σx (MPa)

500
400
300
200
100

I

II III

0
-100
0.000

0.005

0.010

0.015

0.020

0.025

εx

(b)
Figure 4.6: Stress on the cross section at x = 5.0 mm for (a) different horizon δ with the
grid size ∆x = 0.5 mm and (b) different grid size ∆x with the horizon δ = 2∆x. εx is the
engineering strain imposed by the displacement at both ends of the bar.

56

(a)

(b)

(c)

(d)

Figure 4.7: Strain εx at time t = 1.0 × 10−3 s as denoted by I in Figure 4.6(b). The horizon
δ = 2∆x. (a) ∆x = 0.50 mm; (b) ∆x = 0.25 mm; (c) the cross section at z = 3.5 mm,
∆x = 0.50 mm; (d) the cross section at z = 3.5 mm, ∆x = 0.25 mm.

57

(a)

(b)

(c)

(d)

Figure 4.8: Strain εx at time t = 1.6 × 10−3 s as denoted by II in Figure 4.6(b). The horizon
δ = 2∆x. (a) ∆x = 0.50 mm; (b) ∆x = 0.25 mm; (c) the cross section at z = 3.5 mm,
∆x = 0.50 mm; (d) the cross section at z = 3.5 mm, ∆x = 0.25 mm.

58

(a)

(b)

(c)
Figure 4.9: Strain εx at time t = 1.8×10−3 s as denoted by III in Figure 4.6(b). The horizon
δ = 2∆x, the grid spacing ∆x = 0.50 mm. (a) bird’s-eye view of bar; (b) the cross section
at z = 3.5 mm; (c) the cross section at x = 5.0 mm.

59

(a)

(b)

(c)
Figure 4.10: Strain εx at time t = 1.8 × 10−3 s as denoted by III in Figure 4.6(b). The
horizon δ = 2∆x, the grid spacing ∆x = 0.25 mm. (a) bird’s-eye view of bar; (b) the cross
section at z = 3.5 mm; (c) the cross section at x = 5.0 mm.

60

4.4.2

Ductile material

We consider a rectangular bar subjected to a stretch by imposing a constant velocity on
nodes located at the boundary regions such that vx = −10 mm/s at the left-end boundary
region and vx = 10 mm/s at the right-end boundary region as shown in Figure 4.5(a). The
dimensions of the bar are 10 mm in length, 7 mm in width, and 7 mm in thickness. The bar
is uniformly discretized with nodes in a 21 × 15 × 15 grid, and the horizon δ = 2∆x. Young’s
modulus E is 70 GPa, Poisson ratio ν is 0.25, and the density ρ is 2700 kg/m3 . Peridynamic
bonds are assumed to yield at the bond stretch sy = 0.004, and the critical stretch for bonds
to break is set as s0 = 0.02. The calculation time-step is dt = 5 × 10−8 s.
The peridynamics results are compared with the results by the finite element analysis
(FEA) using a commercial FEA software, LS-DYNA [112]. In the FEA, the identical boundary conditions to the peridynamic model are imposed, and the time-step dt is automatically
calculated. The material model used in the FEA is identical to the material model in peridynamics: Young’s modulus E is 70 GPa, and Poisson ratio ν is 0.25. The yield stress
σy is 280 MPa obtained by multiplying the yield bond stretch sy = 0.004 and Young’s
modulus. The limit strain is set to be 0.020, and the corresponding plastic strain is 0.016.
Figure 4.11(a) shows the material properties implemented in LS-DYNA.
Figure 4.11(b) compares the average stress on the cross section at x = 5.0 mm by
peridynamics and the finite element analysis. After the yielding, the cross-sectional stress
continues to increase as the boundary displacement ux increases in the peridynamic model
while the FE model shows almost perfect plasticity. The cross-sectional stress demonstrates
a trilinear response in the peridynamics. As shown in Figure 4.11(b), the end of the second
piecewise linear line approximately corresponds to the displacement ux = 0.017 mm, and
61

the stress is σx = 420 MPa. For the stress of a node in the peridynamic model, the total
bond forces acting through the cross section of a node are calculated. Figure 4.11(c) shows
the comparison of the stress σx of the element which has a node at the center of the bar
(x = 5.0 mm, y = 3.5 mm and z = 3.5 mm) and the stress σx of the peridynamic node located
at the center of the peridynamic bar. Similar to the average stress on the cross section, the
stress of the peridynamic node at the center of the bar shows a trilinear response while the
stress in the element in the FEA has a bilinear response. It is noticed that the outer bonds of
a node yield earlier than the inner bonds in the peridynamic model, and the stiffness of the
peridynamic model gets closer to the value of the FE model after all the bonds yield, as shown
in Figure 4.11(c). The nodal strains εx at the center of the bar (x = 5.0 mm, y = 3.5 mm and
z = 3.5 mm) obtained by peridynamics and FEA are compared in Figure 4.11(d). The cause
of the differences might be the scales at which the material model is defined. Therefore, to
bridge the different scales, it is proposed to modify the scheme to define the material in the
finite element analysis rather than using the identical material model to peridynamics. The
material response obtained in Figure 4.11(b) is interpolated to a trilinear line as shown in
Figure 4.12(a). The point B is determined from the macroscopic response of the peridynamic
model.
A piecewise linear plastic model available in LS-DYNA [112] is used in the FE simulation
p

in which the effective stress versus effective plastic strain εef f curve is defined, where the
p

t
effective plastic strain is given as εef f = 0

2 εp εp dt [72]. The plastic strain starts to
3 ˙ij ˙ij

increase after the strain reaches the yield strain ε = 0.004. The effective plastic strain can
be approximated by integrating the strain rate after yielding, imposed by displacements at
p

both ends. By evaluating the integration, the effective plastic strain εef f corresponding to
p

σx = 420 MPa is obtained as εef f = 0.0037 . The average stress on the cross section at
62

x = 5.0 mm and the stress at the center point of the section are plotted in Figures 4.12(b)
and (c). The nodal strains at the center of the bar are compared in Figure 4.12(d), and the
results show good agreement. As shown in Figures 4.12(b) to (d), the differences between
the solutions by peridynamics and FEA are reduced substantially. The developed scheme
can bridge the simulations at different scales at which the material model is defined. First, a
mesoscale modeling is performed using peridynamics, and material properties are retrieved
from the macroscopic material response. Then the properties are described as a modified
material model at the macroscale for the FEA. Strains εx by peridynamics and FEA when
the displacement at the boundary ux = 0.025 mm (t = 2.5 × 10−3 s) are compared in
Figure 4.13, and the distribution of strains shows almost identical results.

63

σ (MPa)

B

A

280

0.004

0.020

ε
(a)

700

Peridynamics
LS-DYNA

600
σx (MPa)

500
400
300
200
100
0
-100
0.00

0.01

0.02
0.03
ux (mm)

0.04

0.05

(b)
Figure 4.11: A comparison between peridynamic and FEA solutions. (a) The material model
used in FEA which is identical to the constitutive for bonds in peridynamics; (b) Stress σx on
the cross section at x = 5.0 mm; (c) Stress σx on an element which has a node at the center
of the bar (x = 5.0 mm, y = 3.5 mm and z = 3.5 mm) and the stress σx on the peridynamic
node at the center of the bar; (d) Strain εx at the center of the bar (x = 5.0 mm, y = 3.5 mm
and z = 3.5 mm). ux is the displacement of the boundary region as shown in Figure 4.5.

64

Figure 4.11 (cont’d)

700

Peridynamics
LS-DYNA

600
σx (MPa)

500
400
300
200
100
0
-100
0.00

0.01

0.02
0.03
ux (mm)

0.04

0.05

(c)

0.020

εx

0.015
0.010
0.005
0.000
0.00

Peridynamics
LS-DYNA
0.01

0.02
0.03
ux (mm)
(d)

65

0.04

0.05

B

σ (MPa)

420

C

A

280

p

εeff =0.0037

0.004

0.020

ε
(a)

700

Peridynamics
LS-DYNA

600
σx (MPa)

500
400
300
200
100
0
-100
0.00

0.01

0.02
0.03
ux (mm)

0.04

0.05

(b)
Figure 4.12: A comparison between peridynamic and FEA solutions. (a) The modified
material model used in FEA to incorporate retrieved macroscopic material properties from
the peridynamic model; (b) Stress σx on the cross section at x = 5.0 mm; (c) Stress σx on an
element which has a node at the center of the bar (x = 5.0 mm, y = 3.5 mm and z = 3.5 mm)
and the stress σx on the peridynamic node at the center of the bar; (d) Strain εx at the
center of the bar (x = 5.0 mm, y = 3.5 mm and z = 3.5 mm). ux is the displacement of the
boundary region as shown in Figure 4.5.

66

Figure 4.12 (cont’d)

700
600

Peridynamics
LS-DYNA

σx (MPa)

500
400
300
200
100
0
-100
0.00

0.01

0.02
0.03
ux (mm)

0.04

0.05

(c)

0.020

εx

0.015
0.010
0.005
0.000
0.00

Peridynamics
LS-DYNA
0.01

0.02
0.03
ux (mm)
(d)

67

0.04

0.05

(a)

(b)
Figure 4.13: Strain εx when the boundary displacement ux = 0.025 mm. (a) LS-DYNA and
(b) peridynamics.

68

4.5

Summary

A peridynamics code is developed for three-dimensional simulations using numerical micromoduli for discretized peridynamic models. To enhance the efficiency in the calculation,
we use the PGI Accelerator Programming model on an NVIDIA Tesla C1060 GPU for the
parallelization. A benchmark problem of matrix multiplication is conducted to compare the
wall-clock time between the parallelized algebraic computation on the GPU and the serial
computation on the CPU. In order to eliminate the loop dependency in the peridynamics
code for efficient parallelization, an algorithm that separates nested loops into independent
loops is adopted in the GPU implementation. We compared the wall-clock times of the
parallelized calculation in GPU and the serial calculation in CPU for different peridynamic
models. As the number of nodes in the peridynamic model increases, a significant speedup
over the serial calculation is achieved in the GPU.
The responses of the peridynamic model for the horizon δ of 2∆x, 3∆x and 4∆x are
investigated for a brittle material. For the horizon δ = 2∆x, the back-calculated Young’s
modulus matches the exact value of the material. On the other hand, for the horizon larger
than 2∆x, the back-calculated Young’s moduli are smaller than the exact value. With the
horizon δ = 2∆x, the responses of the peridynamic bar discretized with different grid sizes
(∆x = 0.50 mm and ∆x = 0.25 mm) are compared, and numerical results of the backcalculated Young’s modulus show good agreements. The strain distributions at different
times in the simulations are demonstrated. After a large number of bonds reach the critical
stretch, a failure core is formed inside the peridynamic model.
For a ductile material, the peridynamic stresses and strains are compared with the results
by FEA using LS-DYNA. When the material model used in FEA is identical to that of the
69

peridynamic model, the stresses of FEA show bilinear responses while the peridynamic model
demonstrates trilinear responses with the horizon δ = 2∆x. The nodal strains by FEA and
peridynamics at the center of the bar are not identical. Therefore, to bridge material models
at different scales, the material properties are retrieved from the macroscopic response of
the peridynamic model. These retrieved properties are used to modify the material model in
FEA. By applying the proposed procedure, the differences between FEA and peridynamic
solutions reduce substantially. This proposed scheme will be useful for coupling peridynamics
and FEM. Instead of directly applying the material model defined for peridynamic bonds
to FEA in a coupled model, which may not lead to an accurate result, the plastic material
model in FEA might be obtained from the macroscopic responses of peridynamic models.

70

Chapter5
Coupling of Discretized Peridynamics
with Finite Element Method
5.1

Introduction

As a reformulated theory of continuum mechanics, peridynamics eliminates the spatial
derivatives, and it is valid regardless of discontinuities [153, 159]. Therefore, peridynamics is useful to solve problems involving spontaneously emerged discontinuities. Compared
with the finite element method (FEM), peridynamics is computationally expensive. Macek
and Silling [115] implemented peridynamics in a commercial finite element analysis code,
ABAQUS, using truss elements. The conventional FE mesh is coupled with the peridynamic
truss mesh using the embedded element feature available in the finite element analysis code.
Lall et al. [103] used the peridynamics based finite element model to study shock and
vibration reliability of electronics. Kilic and Madenci [94] presented a coupling approach
using overlapping regions in which both peridynamic and FE equations are utilized. Agwai
et al. [2] and Oterkus [133] employed the submodeling approach to couple the FEM with
peridynamics. In their approach, the global analysis by means of finite element analysis
is performed first, and then peridynamics is used for submodeling. A morphing strategy
based on the energy equivalence was proposed by Lubineau et al. [113]. In this chapter,
we introduce a coupling approach of discretized peridynamics with finite elements. Differ71

ent from the approach in [115, 103] implementing peridynamic model in the framework of
the conventional FEM and the submodeling approach [2, 133], the peridynamic subregion
is directly coupled to the finite element subregion in the present approach. An interface
element is introduced to calculate coupling forces instead of using overlapping regions [94] or
the morphing strategy [113] to couple peridynamic and FE subregions. Depending on how
coupling forces are divided to FE nodes of an interface element, we further discuss two types
of coupling schemes.
A great deal of research effort has been made to study many fracture problems. One
common benchmark problem characterized by the mixed mode fracture is the test of a
double-edge-notched concrete specimen conducted by Nooru-Mohamed et al. [130]. The
test of Nooru-Mohamed was adopted by De Borst [38] in the discussion of computational
modeling of concrete fracture. For the analyses, the finite element smeared-crack approach
with the gradient Rankine plasticity model, Cruch-Crack model, and Ottosen’s model has
been used for numerical studies of Nooru-Mohamed’s experiment by Di Prisco et al. [41].
A comparative study of three-dimensional constitutive models for the double-edge-notched
test was performed by Pivonka et al. [140]. Gasser and Holzapfel [64] employed the cohesive
crack model with the PUFEM for the numerical modeling of the test. The XFEM was utilized by Cox [35], Meschke and Dumstorff [121], and Unger et al. [174] for the simulations.
An adaptive mesh refinement technique applied to a nonlocal version of anisotropic damage
model was employed by Patz´k and Jir´sek [136]. R´thor´ et al. [146] used a hybrid ana
a
e
e
alytical and XFEM to study the propagation of curved cracks in the double-edge-notched
concrete specimen.
The remainder of this chapter is organized as follows. In Section 5.2, the finite element
formulations are summarized. In Section 5.3, we present two types of coupling schemes de72

pending on the distributing schemes of coupling forces to FE nodes of interface elements, and
the inverse isoparametric mapping that is essential in the coupling approach is summarized.
In Section 5.4, numerical examples under quasi-static conditions are provided to validate the
proposed coupling approach including one- and three-dimensional elastic problems and the
mixed mode fracture in a double-edge-notched concrete specimen. The results by applying
two types of coupling schemes for numerical simulations are discussed. In the last section,
concluding remarks are summarized.

5.2
5.2.1

Theory
Finite element formulations

The equation of motion in the conventional continuum mechanics is derived from the principle
of linear momentum. The temporal change rate of linear momentum is equal to the force
applied on the body as [118]

ρ¨i = σij,j + fiB
u

in Ω,

(5.1)

where ρ is the density, ui is the acceleration field, σij is the Cauchy stress tensor, and fiB
¨
is the body force density field. The essential boundary condition Γu and natural boundary
condition Γf are defined, respectively, as [11]

ui = Ui

on Γu ,

(5.2)

σij nj = Fis

on Γf ,

(5.3)

73

where the surface of the body Γ = Γu ∪ Γf , Γu ∩ Γf = 0, and nj means the components of
the unit outer normal vector on Γ. The constitutive equation for continuum is stated as

σij = Cijkl εkl ,

(5.4)

where Cijkl is the elastic constitutive coefficient, and the components of strains are defined
as
εij =

1
u + uj,i .
2 i,j

(5.5)

Applying the principle of virtual work to Equation (5.1), we have

Ω

(σij,j + fiB − ρ¨i )δui dΩ = 0,
u

(5.6)

where δui is the virtual displacement. After integrating by parts and applying the divergence
theorem, the weak formulation is obtained as [110]

−

Ω

σij δui,j dΩ +

Γ

σij nj δui dS +

Ω

fiB δui dΩ −

Ω

ρ¨i δui dΩ = 0.
u

(5.7)

Considering the symmetry of the stress tensor (σij = σji ) and applying the boundary conditions, we have

Ω

ρ¨i δui dΩ +
u

Ω

σij δεij dΩ =

74

Γf

Fis δui dS +

Ω

fiB δui dΩ.

(5.8)

Substituting the constitutive law in Equation (5.4) into Equation (5.8), we obtain the finite
element formulation

Ω

ρ¨i δui dΩ +
u

Ω

Cijkl εkl δεij dΩ =

Γf

Fis δui dS +

Ω

fiB δui dΩ.

(5.9)

In the finite element analysis, the displacements within each element are interpolated by
means of shape functions as [11]
u(e) = H(e) U(e) ,

(5.10)

where H(e) is the displacement interpolation matrix, the superscript e represents the element
e, and the nodal displacement vector U(e) is expressed as U(e)T = {u1 v1 w1 · · · un vn wn }
for an element of n nodes. The strain vector is evaluated by

ε(e) = B(e) U(e) ,

(5.11)

where B(e) is the strain-displacement matrix which is written as










(e) = 
B










∂h1
∂x
∂h1
∂y

∂hn
∂y
∂h1
∂z

∂h1
∂y

···

∂h1
∂x
∂h1
∂z

∂h1
∂z



∂hn
∂x

∂hn
∂z
∂hn
∂y

∂h1
∂y

∂hn
∂x
∂hn
∂z

∂h1
∂x

∂hn
∂z

75

∂hn
∂y
∂hn
∂x










.









(5.12)

Substituting Equations (5.10) and (5.11) into Equation (5.9), we have the weak formulation
in matrix form as
¨
MU + KU = F,

(5.13)

where M is the mass matrix, K is the stiffness matrix, and F is the force vector. The
assembled matrices following the convention of direct stiffness method [11] are summarized
as

M(e) ,

M=

M(e) =

e

K(e) ,

K=

K(e) =

e
(e)

e
(e)
Fs =

5.3

e

Γ

(e)
f

Ω(e)

ρ(e) H(e)T H(e) dΩ(e) ,

B(e)T C(e) B(e) dΩ(e) ,

(5.14)
(5.15)

(e)

Fs +

F=

Ω(e)

FB ,
(e)
FB =

H(e)T FS(e) dS (e) ,

(5.16)
Ω(e)

H(e)T f B(e) dΩ(e) .

Coupling between the peridynamic and finite element subregions

5.3.1

Coupling schemes

To gain the efficiency from finite element analyses and exploit the generality of peridynamics in the presence of discontinuities, a domain is partitioned into a conventional FE
subregion and a peridynamic subregion as shown in Figure 5.1. With the coupling approach,
the subregion where failure is expected can be modeled using peridynamics. The overall
computational cost is reduced by using FEM for the remainder of the domain.
76

FE subregion

PD subregion

conventional element
Interface element
PD node
Figure 5.1: Partition of the domain. The FE subregion and the peridynamic (PD) subregion
are bridged by interface elements.
In the conventional FE subregion, the lumped mass matrix is formed by distributing the
total element mass to nodes of the element [11]. The internal forces on FE nodes can be
calculated as
Fint =

Fint(e) =
e

K(e) U(e) ,

(5.17)

e

where K(e) is the element stiffness matrix. The equation of motion of an FE node is obtained
as
¨
MI UI = Fext − Fint ,
I
I

(5.18)

¨
where MI is the lumped mass of node I, UI is the acceleration vector field, Fext is the
I
external force applied on the node I by evaluating the corresponding components of F in
Equation (5.16), and Fint is the internal force vector of the node I. Because the mass matrix
I
is diagonal, Equation (5.18) can be solved without factorizing a global stiffness matrix.
To bridge the FE subregion and the peridynamic subregion, we introduce an interface
element. A three-dimensional interface element consisting of eight FE nodes is illustrated
in Figure 5.2. In an interface element, a number of peridynamic nodes are embedded for
the calculation of coupling forces. The interacting forces between embedded peridynamic
77

1(5)

4 (8)
∆x
δ
Interface element
FE node
PD node
embedded PD node

2(6)

3 (7)

Figure 5.2: Interface element for the coupling of FE subregion and peridynamic subregion.

nodes and peridynamic nodes out of interface elements are called coupling forces. It is
worth noting that interactions between embedded peridynamic nodes are not considered as
coupling forces. The number of embedded peridynamic nodes is determined by the size of
the horizon, and there should be sufficient embedded nodes within the horizon of nodes near
the interface of the peridynamic subregion and the FE subregion as shown in Figure 5.2. To
evaluate coupling forces, each embedded peridynamic node represents a material volume of
(∆x)3 inside an interface element. However, embedded peridynamic nodes are not involved
in the global equation. In other words, the displacements of embedded peridynamic nodes
are not calculated by solving the equation of motion. Therefore, the mass of an interface
element is equally distributed to FE nodes of the interface element. Consider an embedded
peridynamic node subjected to the coupling force f cp as shown in Figure 5.3, the coupling
force is then divided to FE nodes of the interface element by means of shape functions as

cp

fi = φi (ξ, η, ψ)f cp ,
78

i = 1, · · · , 8,

(5.19)

where φi is the shape function of the node i belonging to the interface element, (ξ, η, ψ)
are the natural coordinates of the embedded node in the interface element, which should
be determined by the inverse isoparametric mapping. We designate this coupling scheme
as the VL-coupling scheme since the whole volume of the interface element is subjected
to coupling forces. On the other hand, different from the VL-coupling scheme, we might
divide coupling forces only to the FE nodes on the interface segment as shown in Figure 5.4.
Therefore, FE nodes not on the interface segment are subjected to internal forces arising from
the element stiffness only. Since the interface between the peridynamic and FE subregions
is similar to a contact surface, the scheme demonstrated in Figure 5.4 is designated as
the CT-coupling scheme. To implement the CT-coupling scheme, interfaces between the
peridynamic subregion and the FE subregion have to be defined prior to analyses. Coupling
forces on embedded nodes are divided to those FE nodes on the interface segment as shown
in Figure 5.4 by
cp

fi = φi (ξc , ηc )f cp ,

i = 3, 4, 7, 8,

(5.20)

where φi is the shape function on the interface segment, and (ξc , ηc ) are the natural coordinates of the projection of an embedded node onto the interface segment.
In general, the equation of motion for FE nodes of an interface element is written as

ˆ
¨
MI UI = Fext − Fint ,
I
I

(5.21)

where Fext is the external force by evaluating Equation (5.16), and the internal force is given
I
as
cp
ˆ
Fint = Fint + fI =
I
I

cp

K(e) U(e)
e

79

+ fI ,
I

(5.22)

1(5)

σ

4 (8)

σ

f cp

Interface element
FE node
PD node
embedded PD node

2(6)

3 (7)

Figure 5.3: VL-coupling scheme that divides a coupling force f cp among FE nodes comprising
the interface element.

1(5)

4 (8) interface segment

f cp

σ

σ
Interface element
FE node
PD node
embedded PD node

2(6)

3 (7)

Figure 5.4: CT-coupling scheme that divides a coupling force f cp among FE nodes comprising
the interface segment.

in which [·]I denotes the corresponding components of a vector associated with the node I
cp

of the interface element, and fI is the summation of coupling forces on the node I. The
explicit algorithm is employed for the transient dynamic analyses. Nodal accelerations are
calculated first, and nodal velocities and displacements are updated subsequently. After the
displacements of FE nodes of the interface elements are calculated, the displacements of
80

embedded peridynamic nodes are then determined by

Ueb = φi (ξ, η, ψ)Ui ,

i = 1, · · · , 8,

(5.23)

in which (ξ, η, ψ) are the natural coordinates of an embedded peridynamic node in the
interface element, and Ui is the nodal displacement of an interface element. For peridynamic
nodes out of the interface elements, Equation (3.30) is used to update nodal accelerations.

5.3.2

Inverse isoparametric mapping

To couple peridynamic subregions with finite element subregions, a certain number of peridynamic nodes are embedded in the interface elements. If the Cartesian coordinates of an embedded peridynamic node are known, the natural coordinates of the embedded peridynamic
node in an interface element should be determined by the inverse isoparametric mapping,
which is essential especially for random discretizations. However, the inverse isoparametric
mapping from the Cartesian coordinates to the natural coordinates is nontrivial since equations to be solved are nonlinear. Murti and Valliappan [127] presented a numerical technique
by bisecting a line passing a point and a node of known natural coordinates, and this method
was extended to the three-dimensional space by Murti et al. [128]. However, the bisection
method has its limitations [31]. A more generalized approach for the inverse isoparametric
mapping is presented by Chinnaswamy et al. [31]. For the inverse mapping of an embedded
peridynamic node with known Cartesian coordinates (ˆ, y , z ), the equation to be solved is
x ˆ ˆ
81

written as


 n


 


φi (ξ, η, ψ)xi − x  
ˆ
  
  
 


  0
f1 (ξ, η, ψ)  i=1
  
 

  
  n

  
 


 

φi (ξ, η, ψ)yi − y = 0 ,
ˆ
f = f2 (ξ, η, ψ) =
  
 

  
  i=1

  
  n

  
 

  
 

  0
f (ξ, η, ψ) 


3

φi (ξ, η, ψ)zi − z 
ˆ




(5.24)

i=1

where (ξ, η, ψ) are the natural coordinates of an embedded peridynamic node in an interface
element to be determined. By expanding the vector f in Taylor’s series and omitting the
second and higher order terms, it can be shown that [31]


 ∂f1
 ∂ξ




I = I0 −  ∂f2
 ∂ξ



 ∂f3
∂ξ

∂f1
∂η
∂f2
∂η
∂f3
∂η

−1

∂f1 
∂ψ 





∂f2 
∂ψ 



∂f3 
∂ψ

f0 ,

(5.25)

0

where I0 = [ξ0 , η0 , ψ0 ]T is an approximate solution, and f0 is the vector f evaluated at the
approximate solution I0 . Equation (5.25) can be rewritten as








I = I0 − 






∂φi
∂ξ xi

∂φi
∂η xi

∂φi
∂ξ yi

∂φi
∂η yi

∂φi
∂ξ zi

−1 



∂φi 

xi  


∂ψ 

 

 


∂φi 

∂ψ yi 
 

 

 


 
∂φi


zi

∂ψ

∂φi
∂η zi

0

and it can be simplified as

I = I0 + ∆I,

82






φ i xi − x
ˆ







φ i yi − y  ,
ˆ









φi zi − z 
ˆ

(5.26)

0

(5.27)

where

 
 
 ∆ξ 
 
 
 
 
∆I = ∆η = −J−1 f0 .
0
 
 
 
 
 
∆ψ 

(5.28)

The updated solution of I is used as the value of I0 in Equation (5.27) for the next iteration.
A few iterations are performed till the solution of I converges.

η
ns

t
z

ns’
r

y

∂r
∂η

interface segment

∂r
∂ξ

ξ

x
embeded PD node
FE node
Figure 5.5: Projection of an embedded peridynamic node on an interface segment.

If the CT-coupling scheme is employed, coupling forces on embedded peridynamic nodes
are only distributed to FE nodes on the interface segment as shown in Figure 5.4. Hence, the
natural coordinates of the projection of an embedded peridynamic node onto the interface
segment have to be determined. Let t be the position vector of an embedded node ns , and
the projection of ns onto the interface segment is denoted by ns′ as shown in Figure 5.5.
The position vector r of the point ns′ on the interface segment can be expressed as

r = f1 (ξc , ηc )e1 + f2 (ξc , ηc )e2 + f3 (ξc , ηc )e3 ,
83

(5.29)

where (ξc , ηc ) are the natural coordinates of the point ns′ on the interface segment, and
4

fi (ξ, η) =

j

φ j xi ,

(5.30)

j=1

in which φj is the shape function of the node j on the interface segment. The natural
coordinates (ξc , ηc ) of the point ns′ on the interface segment must satisfy [72]
∂r
· (t − r) = 0,
∂ξ
∂r
· (t − r) = 0.
∂η

(5.31)
(5.32)

However, there is no analytical solution to Equations (5.31) and (5.32). To solve numerically,
a few iterations of the least-squares projection are used to generate an initial guess as

ξ0 = 0,

η0 = 0,

(5.33)



 

 

r,ξ · (r − t)
∆ξ 
r,ξ · r,ξ r,ξ · r,η 
,
=−



 

r · (r − t)
∆η 
r,η · r,ξ r,η · r,η
,η


ξi+1 = ξi + ∆ξ,

ηi+1 = ηi + ∆η.

(5.34)

(5.35)

With an initial guess, Newton-Raphson method is then utilized to find the solution of Equations (5.31) and (5.32) as [187]

H

 
 
∆ξ 
 
∆η 

=−





r,ξ · (r − t)


r · (r − t)
,η

84

,

(5.36)









r,ξ · r,ξ r,ξ · r,η  (r − t) · r,ξξ (r − t) · r,ξη 
H=
+
,
r,η · r,ξ r,η · r,η
(r − t) · r,ξη (r − t) · r,ηη
ξi+1 = ξi + ∆ξ,

ηi+1 = ηi + ∆η.

(5.37)

(5.38)

The solutions of ξi+1 and ηi+1 are used to update the value of the position vector r, and
then Equation (5.36) is evaluated again for the next iteration. The converged solutions of
ξi+1 and ηi+1 are the natural coordinates of the point ns′ on the interface segment.

5.4
5.4.1

Numerical applications
One-dimensional bar

For benchmarking, the present coupling approach of peridynamics with FEM is employed
to study the axial deformation of a one-dimensional bar. The length of the bar is 9.5 mm,
and dimensions of the cross section are 0.5 mm by 0.5 mm. Young’s modulus E of the bar is
70 GPa, and the density is 2700 kg/m3 . Figure 5.6 shows the multiscale discretization of the
bar. The finite element mesh size is 1.5 mm, and the conventional bar elements are utilized.
The stiffness of a bar is k = EA , where A is the cross-sectional area and L is the length of
L
the bar. The peridynamic grid spacing is 0.5 mm, which is equal to the width of the bar.
By setting the horizon to 2∆x, the one-dimensional micromodulus c1 is 5.6 × 1023 N/m6 .
Two interface elements are used to couple peridynamic and finite element subregions, and
each interface element has two embedded peridynamic nodes for the calculation of coupling
forces as shown in Figure 5.6.
85

9.5 mm

F

F
δ = 2∆x

x
Interface element
FE node
PD node
embedded PD node

Figure 5.6: Discretization of a one-dimensional bar.
A tensile loading with the magnitude of F = 175 N is applied on both ends of the bar under the quasi-static condition. The force is gradually increased during 50, 000 steps with the
calculation time step dt = 5 × 10−8 s, which is less than the critical time step for the explicit
time integration. Figure 5.7(a) shows the displacement along the bar using the VL-coupling
scheme. Nodal displacements in the peridynamic subregion show good agreement with the
quasi-static solution. The displacements of the FE nodes, however, show small discrepancies.
The reason is that coupling forces in an interface element are divided among all FE nodes of
the element using the VL-coupling scheme so that FE nodes at the interfaces of subregions
(x = ±1.75 mm) only receive partial coupling forces. FE nodes at the other end of the
interface elements (x = ±3.25 mm) receive the rest of coupling forces, and are also subjected
to internal forces contributed by the element stiffness. Consequently, the displacements of
FE nodes at the interfaces are slightly overestimated, and displacements of other FE nodes
are underestimated as indicated in Figure 5.7(a). In contrast, the solution using the CTcoupling scheme, which distributes coupling forces only to FE nodes at interfaces, is almost
identical to the quasi-static solution as shown in Figure 5.7(b). Hence, for the calculation of
axial displacement, the CT-coupling scheme is more effective than the VL-coupling scheme
to achieve the coupling between the peridynamic subregion and the finite element subregion.

86

0.05
0.04
Displacement (mm)

0.03
0.02
0.01
0.00

interface

interface

-0.01
-0.02
-0.03

quasi-static solution
FE node (coupling model)
PD node (coupling model)

-0.04
-0.05
-5

-4

-3

-2

-1

0
1
x (mm)

2

3

4

5

4

5

(a)
0.05
0.04
Displacement (mm)

0.03
0.02
0.01
0.00

interface

interface

-0.01
-0.02
-0.03

quasi-static solution
FE node (coupling model)
PD node (coupling model)

-0.04
-0.05
-5

-4

-3

-2

-1

0
1
x (mm)

2

3

(b)
Figure 5.7: Axial displacement along the bar using (a) VL-coupling scheme and (b) CTcoupling scheme.

87

5.4.2

Three-dimensional bar

A three-dimensional bar subjected to tension is examined to compare the solutions of the
present coupling approach and the classical (local) elasticity solutions. The dimensions of
the bar are taken to be 10 mm in length, 7 mm in width, and 7 mm in thickness as shown
in Figure 5.8. The three-dimensional model is partitioned into two FE subregions and
one peridynamic subregion. Each FE subregion consists of four eight-node solid interface
elements, and the mesh size of the interface element is 3.5 mm. The peridynamic subregion
is uniformly discretized with the grid spacing ∆x = 0.5 mm, and the size of the horizon is
set to δ = 1.0 mm. In the interface element, two additional layers of peridynamic nodes
are embedded along the longitudinal direction to ensure sufficient nodes in the horizon of
peridynamic nodes near the interfaces. Tractions at both end surfaces are gradually applied
up to σx = 700 MPa during 100, 000 steps as quasi-static loading, and the calculation time
step is set to dt = 5 × 10−8 s. Young’s modulus of the bar is 70 GPa, Poisson ratio is 0.25,
and the density is 2700 kg/m3 .
We first examine the numerical solutions using the VL-coupling scheme. The longitudinal displacements ux at three measuring positions are compared to the quasi-static results
as shown in Figure 5.9(a). The displacement ux at the position p1 , which is the at the
end surface, is underestimated compared with the quasi-static solution. On the other hand,
the displacement at p2 , which is at the interface of the peridynamic subregion and the FE
subregion, and the displacement at p3 , which is inside the peridynamic subregion, show
good agreement with the quasi-static solutions. To look into Poisson effect, the transverse
displacements measured at different positions are plotted in Figure 5.9(b). The transverse
displacement q1 at the end surface is smaller than the quasi-static value since the longitudinal
88

7m

m

Interface element PD subregion Interface element

y

σx

x

σx

7 mm

z

3.5 mm

3 mm

3.5 mm

Figure 5.8: Three-dimensional bar subjected to tension.

displacement at the end surface is underestimated. Nevertheless, the transverse displacement
q2 at the interface demonstrates close agreement with the quasi-static solution, which indicates that Poisson effect is preserved at the interface. The transverse displacement at q3 ,
which is inside the peridynamic subregion, also matches the quasi-static solution.
For the comparison of two types of coupling schemes, the CT-coupling scheme, which
achieves coupling by considering the interfaces of subregions similar to contact surfaces, is
then employed to study this quasi-static problem. As shown in Figure 5.10(a), the longitudinal displacement p1 at the end surface is close to the quasi-static solution with the error less
than 2%. The longitudinal displacements at interface and inside the peridynamic subregion,
denoted by p2 and p3 respectively in Figure 5.10(a), are almost identical to the quasi-static
solutions. The transverse displacements measured at different positions are plotted in Figure 5.10(b). With the improved solution of longitudinal displacement at the end surface,
89

the transverse displacement q1 at the end surface turns to be very close to the quasi-static
solution. On the other hand, the transverse displacement q2 at the interface is overestimated. This phenomenon occurs due to the reason that decompositions of coupling forces
in the transverse direction are divided only among FE nodes at the interface. By comparing
Figures 5.9 and 5.10, it is observed that the CT-coupling scheme is effective in resolving displacements normal to the interface of peridynamic and FE subregions. On the other hand,
the VL-coupling scheme, which divides coupling forces among all FE nodes comprising the
interface element, is capable to preserve Poisson effect at the interface.
Figure 5.11 shows the longitudinal displacement along the edge of the bar solved by
the CT-coupling scheme. It is noted that a smooth curve can be obtained if the nodal
displacements are connected, which is different from the results in [94] where a jump in
displacement is observed at the interface. The strain εx distributions on the bar and at the
interface are plotted in Figure 5.12. Nodal strains in the FE subregions are obtained by
evaluating Equation (5.5), and nodal strains in the peridynamic subregion are calculated as
the average bond stretch within the grid width ∆x. Considering the quasi-static condition,
the analytical value of strain εx is 0.01. The strain of the coupling model is in the range of
0.0099 to 0.0105 as shown in Figure 5.12, which agrees with the quasi-static solution.
The essence of coupling forces is interactions between nodes in the peridynamic subregion
and material volumes of an adjacent continuous body (i.e. interface elements) represented
by embedded peridynamic nodes. Therefore, the ratio of the grid spacing of peridynamic
nodes to the mesh size of interface elements should not have affect the results. To illustrate
this point, Figure 5.13 shows the summation of coupling forces in the longitudinal direction
on the interface elements at the left end of the bar as the applied traction increases. As
indicated by the comparison in Figure 5.13, differences in the results using the grid spacing
90

∆x = 1.0 mm and ∆x = 0.5 mm are insignificant.

91

0.05

quasi-static solution
coupling model

Displacement (mm)

0.04
0.03

p1(5,0,0)

0.02
0.01

p2(1.5,0,0)
p3(0.25,0.25,0.25)

0.00
0

100

200
300
400
500
Applied traction (MPa)

600

700

(a)
0.000

Displacement (mm)

q3(0.25,0.25,0.25)

-0.003
q1(5,3.5,0)
q2(1.5,3.5,0)

-0.006

quasi-static solution
coupling model

-0.009
0

100

200
300
400
500
Applied traction (MPa)

600

700

(b)
Figure 5.9: Displacements (a) ux and (b) uy at different locations using the VL-coupling
scheme. Coordinates of the measuring position are given in parenthesis.

92

0.05

quasi-static solution
coupling model

Displacement (mm)

0.04
0.03

p1(5,0,0)

0.02
0.01

p2(1.5,0,0)
p3(0.25,0.25,0.25)

0.00
0

100

200
300
400
500
Applied traction (MPa)

600

700

(a)
0.000

Displacement (mm)

q3(0.25,0.25,0.25)

-0.003

q2(1.5,3.5,0)

-0.006

q1(5,3.5,0)

quasi-static solution
coupling model

-0.009
0

100

200
300
400
500
Applied traction (MPa))

600

700

(b)
Figure 5.10: Displacements (a) ux and (b) uy at different locations using the CT-coupling
scheme. Coordinates of the measuring position are given in parenthesis.

93

0.06

Displacement (mm)

0.04
0.02
interface

0

interface

-0.02
-0.04

quasi-static solution
coupling model (FE nodes)
coupling model (PD nodes)

-0.06
-6

-4

-2

0
x (mm)

2

4

6

Figure 5.11: Axial displacement along the edge using the CT-coupling scheme.

94

(a)

(b)
Figure 5.12: Strain εx distributions (a) on the bar and (b) on the interface when the applied
traction σx = 700 MPa. The CT-coupling scheme is used in the simulation.

95

35

Coupling force (kN)

30
25
20
15
10
5

∆x=0.5 mm
∆x=1.0 mm

0
0

100

200
300
400
500
Applied traction (MPa)

600

700

Figure 5.13: Summation of coupling forces in the longitudinal direction on the interface
elements at the left end of the bar. The CT-coupling scheme is used in the simulation.

96

5.4.3

Mixed mode fracture in a tension-shear specimen

A benchmark problem of mixed mode crack propagation in a concrete specimen has been
investigated experimentally by Nooru-Mohamed et al. [130]. The double-edge-notched specimen is illustrated in Figure 5.14(a). The dimensions of the specimen are taken to be
200 mm in both length and height, 50 mm in width, and two notches at edges are 25 mm
in length, 5 mm in height, and 50 mm in width. For the numerical study, the specimen
is partitioned into two finite element subregions and one peridynamic subregion as shown
in Figure 5.14(b). The FE subregion is discretized with two mesh sizes, 10 mm× 10 mm×
16.25 mm and 10 mm× 10 mm× 13 mm, respectively. The peridynamic subregion is discretized with the grid spacing ∆x = 5.0 mm, which is 1/10 of the specimen thickness, and
the horizon is set to δ = 1.0 mm. The notches in the peridynamic subregion are introduced
by deleting nodes along two notches and removing all bonds across the notches. Since the
ratio of the horizon to the grid spacing is equal to two, two layers of peridynamic nodes adjacent to interfaces are embedded in interface elements for the calculation of coupling forces.
Young’s modulus, Poisson ratio, and fracture energy were not measured in the experiments.
Therefore, we adopt the material properties E = 30 GPa and Gf = 110 J/m2 as in [35]. For
Poisson ratio, it is estimated as ν = 0.2 in the numerical studies in [121, 35] and ν = 0.3
is used in [125]. In the present study, Poisson ratio ν = 0.25 is assumed. The density is
calculated from the concrete composition given in [130] as ρ = 2265 Kg/m3 . By applying
Equation (4.3), the critical bond stretch is obtained as s0 = 5.5277 × 10−4 . The brittle
material model illustrated in Figure 4.1(a) is utilized for the peridynamic subregion, and the
linear elastic model is employed for the remaining FE subregions, where material parameters
E = 30 GPa and ν = 0.25 are used.

97

un

5 mm

5 mm

200 mm

F

25 mm

25 mm

F

50

un

mm

200 mm

32.5 mm

(a)

135 mm

FE subregion

32.5 mm

PD subregion

FE subregion
10 mm

10 mm
200 mm

(b)
Figure 5.14: Mixed mode fracture test: (a) geometry of the specimen; (b) subregions of the
coupling model.

98

For the comparison, the load-path 1b (specimen 46-05) in [130] is considered. A horizontal shear force of 10 kN is applied first, and then the vertical displacement un is applied
on the top and bottom of the specimen as shown in Figure 5.14(a). For the numerical calculation, the time step is set to dt = 1 × 10−7 s, which is less than the critical time step
for the explicit time integration. The vertical displacement un is applied by imposing a
constant velocity of 10 mm/s as quasi-static loading. We first examine the numerical predictions of crack paths using the VL-coupling scheme. The cracks initiate near the notches
as shown in Figure 5.15(a), and propagate along the horizontal direction for about 50 mm.
As the boundary displacement increases, the direction of propagation changes as shown in
Figures 5.15(b) and 5.15(c). At the boundary displacement un = 0.09 mm, two cracks are
connected as shown in Figure 5.15(d). The numerical prediction of crack paths using the
VL-coupling scheme shows differences with the experimental observations in [130].
We then apply the CT-coupling scheme for the numerical simulation. Interfaces are
predefined in the reference configuration, and there are 280 interface segments with the given
mesh sizes. Natural coordinates of the projections of embedded nodes onto interface segments
are saved to a list. In the subsequent calculations, the natural coordinates of projected points
are referred, and the coupling forces on embedded peridynamic nodes are divided among FE
nodes at interface segments. The damage evolution is shown in Figure 5.16. Crack initiation
occurs at the left and right notches as shown in Figure 5.16(a), and materials ahead of crack
tips are damaged for the length of around 15 to 20 mm. Due to the angle change of the
principle stress, cracks propagate with an angle as shown in Figure 5.16(b). As the boundary
displacement increases, two curvilinear crack paths are clearly observed in Figure 5.16(c). An
enclose area is gradually formed between two curvilinear cracks as indicated in Figure 5.16(d).
The numerical prediction of cracks shows agreement, especially for the lower crack, with
99

the experimental observation presented in [130], which is illustrated using the solid line in
Figure 5.16(d). The small discrepancies appeared in Figure 5.16(d) might be caused by the
perfectly brittle material model adopted in the numerical simulation. More investigation
of material models for concrete in the framework of peridynamics is required in the future.
Note that the numerically predicted crack paths are in perfect symmetry. On the contrary,
the cracks in the experiments do not show perfect symmetry, and we might speculate on
symmetry breaking caused by the slight differences in the loading condition and material
heterogeneity.
Compared with the numerical results using the VL-coupling scheme, the results using the
CT-coupling scheme demonstrate better agreement with the experimental observation. The
reason for the inferiority of the VL-coupling scheme in this mixed mode fracture problem is
that the horizontal displacement at edges caused by the constant shear loading is underestimated after the applied vertical displacement at boundaries reaches the critical value for the
crack initiation at two notches. Consequently, the opening-mode fracture dominates, and
the large area ahead of crack tips is damaged in the plane of notches. After the intact region
in the middle of the specimen becomes relatively small, rotations of crack paths then take
place as shown in Figure 5.15.

100

(a)

(b)

(c)

(d)

Figure 5.15: Numerical prediction of crack paths using the VL-coupling scheme. Boundary
displacement: (a) un = 0.0220 mm; (b) un = 0.0225 mm; (c) un = 0.03 mm; (d) un =
0.09 mm.

101

(a)

(b)

(c)

(d)

Figure 5.16: Numerical prediction of crack paths using the CT-coupling scheme. Solid
curves are crack paths observed in the experiments from [130]. Boundary displacement: (a)
un = 0.0220 mm; (b) un = 0.0225 mm; (c) un = 0.03 mm; (d) un = 0.09 mm.

102

5.5

Summary

A coupling approach of discretized peridynamics with FEM is presented in this chapter.
To bridge conventional FE subregions and peridynamic subregions, an interface element is
introduced. The proposed coupling approach is different from other methods in the sense of
direct coupling via interface elements. Depending on the size of the horizon, a number of
peridynamic nodes are embedded in an interface element. The embedded peridynamic nodes
are not involved in the global equation, but essential in the calculation of coupling forces.
The coupling forces describe interactions between embedded peridynamic nodes in interface
elements and peridynamic nodes in peridynamic subregions. Two types of coupling schemes
are introduced. In the VL-coupling scheme, coupling forces on embedded peridynamic nodes
are divided to FE nodes of interface elements. On the other hand, coupling forces are divided
only to FE nodes on interface segments in the CT-coupling scheme. The inverse isoparametric
mapping techniques to determined the natural coordinates of embedded peridynamic nodes
in the interface elements and the natural coordinates of projected points on the interface
segments are summarized.
Numerical simulations are conducted to compare the computational results using the
coupling approach to the classical elasticity solutions. The axial deformation of a onedimensional bar under quasi-static loading is studied. It is found that the displacements at
interfaces of subregions are slightly overestimated, and the displacements in the FE subregions are underestimated using the VL-coupling scheme. On the other hand, the numerical
solution using the CT-coupling scheme is almost identical to the quasi-static solution in all
subregions and interfaces. For three-dimensional simulations, a bar subjected to quasi-static
tension is partitioned into a peridynamic subregion and two FE subregions consisting of
103

eight-node interface elements at both ends of the bar. By measuring displacements at different positions, the CT-coupling scheme is found to be effective in resolving displacements
normal to the interface of peridynamic and FE subregions, and the VL-coupling scheme is
capable to preserve Poisson effect at the interfaces. Longitudinal strain distributions in the
bar and at the interface using the CT-coupling scheme demonstrate good agreement with
the quasi-static solution.
The last numerical example is the mixed mode fracture in a concrete specimen subjected
to quasi-static loading. The region where failure is expected is modeled using peridynamics,
and the remaining region is modeled using conventional FEM to reduce the computational
cost. Numerical predictions of crack paths using the VL-coupling scheme and the CTcoupling scheme are studied. Two independent curvilinear crack paths are observed in the
results using the CT-couping scheme, and the numerical predictions of crack patterns are
close to the experimental observations presented in [130].

104

Chapter6
Discretized Peridynamics for ContactImpact Problems
6.1

Introduction

Impact, which involves the collision of two or more objects, encompasses a wide range of phenomena such as automobile accidents, drop of electric devices, and even molecular collisions
[66]. When bodies collide, the contact force is developed to prevent them from overlapping
in material regions. A compatible contact surface enveloping the initial contact points is
formed by the reaction force [165]. The stress waves developed by the collision travel away
from the region of contact [66]. As the initial velocity increases, penetration of targets occurs.
Penetration involves perforation, embedment, and ricochet [9].
Although the dynamic behavior of materials in the process of impact is complicated,
only a few highly influential mechanisms play the governing roles [145]. Therefore, it is
possible to use relatively simple analytical models to describe the impact process. Bishop
et al. [17] first studied the analytical methods for penetration mechanics. Recht and Ipson
[145] developed analytical equations of ballistic perforation dynamics, and confirmed the
proposed analytical models by comparing with experimental data. Tate [170, 171] modified
the hydrodynamic theory to incorporate strength effects for the prediction of the deceleration
of a long rod penetrating into a target, and studied the deformation of a soft rod striking
105

a rigid target. Jones et al. [88] modified the one-dimensional penetration theory proposed
by Tate by taking account of mass transfer and the mushroom strain into the analysis.
Cinnamon et al. [32] presented a one-dimensional analysis to predict penetration depths and
profile hole diameters for rod penetrators in semi-infinite targets. Grace [68] proposed a onedimensional model to predict the penetration of long rods into targets. Warren and Forrestal
[179] presented penetration models using the spherical cavity-expansion approximation for
rigid spherical-nosed rods that penetrate aluminum targets. Analytical models that describe
forces and penetration depths for rigid long rod and aluminum targets were developed by
Forrestal et al. [55, 54, 57]. In addition, Forrestal et al. [52, 53] proposed a dimensionally
consistent empirical equation for the depth of penetration into concrete targets based on
experimental data, and studied the deceleration of penetration in concrete targets. An
empirical relationship describing the penetration hole diameter in thin plate was proposed
by Hill [76]. A review of simplified analytical models for ballistic penetration into different
media was presented by Ben-Dor et al. [16].
Most research progress on penetration has been made in experimental investigations, and
a large number of studies can be found. Stock and Thompson [164] studied the penetration
of aluminum alloys by projectiles and showed that the formation of bands of intense shear
reduces the ability of the material to withstand further penetration. Calder and Goldsmith
[27] studied the plastic deformation and perforation of thin plates subjected to projectile
impact using a high speed framing camera. Doyle [42] performed experimental measurements of the contact force during the transverse impact of an aluminum plate using dynamic
strain gage, and comparisons were made with the finite element analysis. Bless et al. [19]
studied hypervelocity penetration of ceramics and penetration phenomena, which are dependent on ceramic properties. The penetration of hard layers by rod projectiles was analyzed
106

by Bless and Anderson [18]. The penetration of long rod into aluminum was investigated
under different impact velocities by Hohler et al. [77]. Trucano and Grady [173] presented
experimental, analytical, and computational techniques to study impact shock and penetration in low density media. Subramanian et al. [166] conducted reverse impact experiments
for aluminum penetrations to study the penetration rate, consumption velocity, and total
penetration. Rosenberg et al. [148] carried out ballistics experiments of ceramic tiles and
studied a simplified Johnson-Holmquist failure mode for brittle materials. Murr et al. [126]
studied the low velocity to hypervelocity penetration transition in metal targets. Frew et
al. [62, 61] conducted depth-of-penetration experiments of steel projectile into concrete, and
conducted experiments in which penetration depths of limestone targets were measured and
compared to an analytical penetration equation. Piekutowski et al. [139] performed depthof-penetration experiments using steel projectiles and aluminum targets. Similarly, Forrestal
and Piekutowski [56] conducted penetration experiments with 6061-T6511 aluminum targets
and spherical-nose steel projectiles. Lundberg et al. [114] studied the interaction of a metallic projectiles and a ceramic target. Gomez and Shukla [67] performed an experimental
study of projectile penetration into concrete. Warren and Poormon [180] presented experimental, analytical, and numerical results of penetration of aluminum targets by ogive-nosed
projectiles at oblique angles. Børvik and colleagues [24, 21, 23, 22] examined the ballistic
penetration of steel plates by cylindrical projectiles by experimental, analytical, and numerical investigations, and studied a constitutive model for non-linear finite element analyses.
Behner et al. [12] experimentally studied the penetration and failure of glass by rod impact.
Anderson et al. [7] investigated the failure and penetration response of glass impacted by
short rods so that there is no steady driving stress.
Due to the complexity and the experimental cost in laboratory tests, researchers resort
107

to effective numerical methods for the analyses of contact-impact problems. Wilkins [185]
investigated influence of material properties on perforation of targets using a finite difference
program. Anderson et al. [6] used Eulerian wave propagation program to examine the longrod penetration as a function of impact velocity. Taylor et al. [172] used the hydrocode
modeling to study hypervelocity impact on brittle materials. Warren and Tabbara [182]
performed simulations of penetration without the need for discretizing the target and the
need for a contact algorithm. Heinstein et al. [75] proposed a method for the contact-impact
modeling of large deformation in explicit transient dynamics. Rosenberg and Dekel [147]
presented two-dimensional simulations to investigate penetration of rigid short projectiles
into semi-infinite targets and perforation of thin metal plates.
As a well-established method for analyzing a wide range of engineering problems, the
finite element method has been employed to study contact and impact problems. Hughes et
al. [84] presented the theory and numerical implementation of the finite element method for
contact-impact problems using a node-to-node scheme. A numerical method for the analysis
of contact between three-dimensional bodies was presented by Chaudhary and Bathe[29].
Belytschko and Lin [14] presented an algorithm for three-dimensional impact simulations,
and Belytschko and Neal [15] developed a pinball algorithm as a simplified contact-impact
algorithm that can be readily vectorized. Carpenter et al. [28] developed an approach
for the fulfillment of surface contact conditions in the transient non-linear finite element
analysis. Anderheggen et al. [4] introduced a numerical contact algorithm with an explicit
time integration scheme in the finite element method. Hahn and Wriggers [71] proposed an
explicit multi-body contact algorithm. Elabbasi et al. [47] addressed reliability issues that
are necessary for analyzing the accuracy and effectiveness of engineering designs involving
contact conditions.
108

The primary focus in this chapter is the numerical scheme for the modeling of contact
between a peridynamic domain and a non-peridynamic domain such as conventional finite
elements and rigid bodies. Macek and Silling [115] presented a plate perforation example,
in which the peridynamic model is implemented using truss elements and a kinetic contact
algorithm available in the commercial FEA package is employed. Littlewood [107] implemented contact modeling between the peridynamic and finite element portions by an iterative
approach. The contact algorithm utilized in [107] has to work with planar facets, and consequently each peridynamic node is represented by an icosahedron. Different from approaches
in [115, 107], the contact between a discretized peridynamic domain and a non-peridynamic
domain is directly modeled as a node-to-surface type in the present research, and a penalty
method without iterations is used for transient analyses by the explicit time integration.
The remainder of this chapter is organized as follows. First, a short-range force model
is introduced into the peridynamic theory. The inverse isoparametric mapping, which is
useful to determine whether a peridynamic node penetrates into an element, is summarized.
A contact algorithm for the calculation of the contact point and the penetration depth
is then explained. For the transient analyses by the explicit time integration, a penalty
method is employed to enforce displacement constraints on the contact surfaces. In the
numerical examples, the material behavior of a peridynamic model is examined. Then the
impact between two rigid bodies is presented to validate the contact algorithm. The last
numerical example is a steel plate, which is modeled using peridynamics, perforated by a
blunt-nosed cylindrical rigid projectile, and the residual velocities are compared with the
analytical values.

109

6.2

6.2.1

Theory

Peridynamic short-range forces

In the peridynamic theory introduced in Section 3.2, material points interact through predefined bonds. However, in applications involving damage, two material points, which are
beyond the horizon of each other in the reference configuration, might come into contact. In
order to consider contact forces between material points, a short-range force model [115] is
incorporated into the peridynamic theory as

fs =

η+ξ
min 0, cs
η+ξ

η+ξ
−1
ds

,

(6.1)

where cs is the spring constant for contact forces, and ds is a constant such that the contact
force is nonzero if the distance between two material points is smaller than ds . In numerical
implementations, cs is set to be 15c as in [115], and ds is chosen to be ds = 0.9∆x, where
∆x is the grid spacing for the discretization.
The peridynamic equation of motion after discretization is written as

ρ¨ t =
uI

NH

Ns

I

J=1

f VJ +

fs VK + bt ,
I

(6.2)

K=1

where ut is the acceleration of the node I at time t, NH is the total number of nodes within
¨I
I

the horizon of the node I in the reference configuration, Ns is the number of nodes within
the contact force range ds of the node I, V is the nodal volume, and bt is the body force at
I
time t.
110

6.2.2

Inverse isoparametric mapping

The isoparametric finite element compared to the generalized coordinate finite element is
more effective in evaluating element matrices [11]. The physical element is mapped into a
reference element that is a square or a cube, and there exists a one-to-one mapping between
the Cartesian coordinate system (x, y, z) and the natural coordinate system (ξ, η, ψ) [34].If
the natural coordinates of a point in an element are known, the Cartesian coordinates can
easily be calculated as
n

x=

Ni (ξ, η, ψ)xi ,

(6.3)

i=1

where Ni are shape functions, and xi represent the Cartesian coordinates of the node i.
However, the inverse isoparametric mapping from the Cartesian coordinates to the natural
coordinates is nontrivial, and it is often required in the dynamic analyses due to, for example,
the remeshing in the crack propagation [127]. In the contact-impact problems, the inverse
isoparametric mapping is useful to determine whether a node penetrates into an element.
The details of three-dimensional inverse isotropic mapping are introduced in Section 5.3.2

ns
s2
s3

ms

s1
s4

z
y
x
Figure 6.1: Slave node and four master segments.

111

6.2.3

Contact algorithm

Consider a slave node ns lies above a master surface as shown in Figure 6.1, and the master
surface consists of four segments s1 , s2 , s3 , and s4 . In order to find the contact point and
the corresponding master segment, the first step is to search for the master node ms that is
nearest to the slave node ns . The vector beginning at ms and ending at ns is denoted by g
as shown in Figure 6.2. Vectors c1 and c2 are along edges of a segment, and point outward
from the master node ms . The unit normal vector of the segment can be calculated by

m=

c1 × c2
,
|c1 × c2 |

(6.4)

and the projection of the vector g onto the segment is expressed as

s = g − (g · m)m.

(6.5)

A slave node ns lies on the segment if the following criteria are satisfied [72]:

(c1 × s) · (c1 × c2 ) > 0,

(6.6)

(c1 × s) · (s × c2 ) > 0.

(6.7)

However, if the slave node is above the intersection of two master segments as shown in
Figure 6.3, the criteria given in Equations (6.6) and (6.7) would fail to be satisfied. Instead,
the segment intersection on which the slave node lies can be determined by finding the
112

11
00g n
11
00 n
c
11
00 s
1111
0000
m 00
11
c

s

2

s’

s

1

z
y
x
Figure 6.2: Projection of the vector g onto a master segment.
maximum value of the following expression [72]:

g · ci
,
|ci |

i = 1, 2, · · · , n,

(6.8)

where n is the number of segment intersections at ms . In the case shown in Figure 6.3, the
value of n is 4.

111
000 n
c
111
000
111
000n c
m
11111
00000
s

2

s

1

s’

z
y
x
Figure 6.3: Slave node projected on the intersection of two master segments.

Assume that a master segment has been detected for the slave node ns as shown in
Figure 6.4. Then the contact point ns′ , which is the nearest point on the master segment to
the slave node, needs to be determined. The method to find the contact point is the same
as the technique introduced in Section 5.3.2
113

η
ns

t
z

∂r
∂η
ns’

r

y

∂r
∂ξ

ξ

x
Figure 6.4: Location of a contact point.
Penetration of the slave node ns through the master segment occurs if

l = n · (t − r) < 0,

(6.9)

where l is the penetration depth, and n is the unit normal vector of the master segment
containing the contact point.

6.2.4

Penalty method

For transient analyses by the explicit integration, the penalty method is employed to enforce the displacement constraints at contact surfaces. The contact between a peridynamic
domain and a non-peridynamic domain can be conveniently treated as a node-to-surface
contact. The nodes in the discretized peridynamic domain are assumed to be slave nodes,
and the surfaces of the non-peridynamic domain, which might come into contact with the
peridynamic domain, are considered as master surfaces. The standard penalty formulation
[72] is usually utilized to calculate the contact force, which is dependent on material constants. However, the standard penalty formulation is mainly suitable for similar materials
coming into contact. For contact between dissimilar materials, the segment-based penalty
114

formulation [72] is recommended, and the contact stiffness is given as


 Ξs




kcs = 0.5 × Ξp ×
 or




 Ξ
m







 m m
1
1 2
,
 m1 + m2 dt2






(6.10)

where Ξp is the penalty scale factor which is 0.1 by default, Ξs and Ξm , which are 1.0 by
default, are scale factors of slave and master penalty stiffnesses respectively, m1 is the nodal
mass of the slave node, and m2 is the segment mass which can be taken to be half of the
element mass [72].
If the slave node ns penetrates the master segment, i.e. l < 0 in Equation (6.9), a
penalty force is then applied to the slave node and the master segment. The force, which is
proportional to the penetration depth, is given by

Fp = (kcs × l)n,

(6.11)

and the penalty force applied to each node of the master segment can be determined using
shape functions as
Fi = φi (ξc , ηc )Fp ,
p

(6.12)

where i indicates nodes that comprise the master segment.
Nodal accelerations are calculated by considering the internal forces, external forces,
and contact forces. Subsequently, nodal velocities and displacements are updated using the
Velocity-Verlet algorithm. The flowchart of the explicit algorithm for transient dynamic
analyses is summarized in Table 6.1, which is similar to the scheme presented in [15].

115

1. Initialize.
2. Calculate the external forces.
3. Calculate the internal forces.
(a) In the peridynamic domain, compute the internal forces from bond forces.
(b) In the non-peridynamic domain, e.g. finite elements, compute the internal
forces arising from the element stress.
4. Calculate the contact forces.
5. Assemble external forces, internal forces, and contact forces into a force array.
6. Compute nodal accelerations.
7. Integrate the accelerations to calculate the nodal velocities and displacements.
8. Go to step 2 for the next time step.
Table 6.1: Explicit contact algorithm.

projectile

target

d

shear zone

Figure 6.5: Perforation of a thin plate by a blunt projectile.

6.3

Ballistic limit

The perforation of a plate depends on the geometry and the strength of material [145].
Ballistic limit velocity is the minimum perforation velocity [145]. During perforation of a
plate, damage develops in the shear zone surrounding the projectile as shown in Figure 6.5,
and a plug in the target plate is formed by indentation. The residual velocity of a blunt
116

projectile after perforating a thin plate is expressed as [145]

p

p

vr = a(vi − vbl )1/2 ,

(6.13)

where vi is the initial velocity, vbl is the ballistic limit velocity, the constant p is equal to 2,
and the mass ratio is
a=

mp
,
mp + mt

(6.14)

where mp is the mass of the projectile, and mt is the mass of the plug formed by perforation.
The value of mt can be approximated by [24]

mt ≈

π 2
d ρt h t ,
4

(6.15)

where d is the diameter of the cylindrical projectile, ρt is the density of the target plate, and
ht is the thickness of the plate. On the other hand, constants a and p in Equation (6.13)
can also be determined by fitting the experimental data [24].

boundary region

22.5 mm

boundary region

y
x

vx
m

z

.5

22
9 mm

13.5 mm

9 mm

Figure 6.6: Peridynamic bar subjected to tension.

117

m

vx

6.4
6.4.1

Numerical studies
Material behavior and modeling

Consider a peridynamic bar subjected to tension as shown in Figure 6.6. The dimensions of
the bar are taken to be 31.5 mm in length, 22.5 mm in width and thickness. The boundary
regions shown in Figure 6.6 are subjected to a constant velocity vx in the longitudinal direction. Magnitudes of the velocity vx are 10 mm/s for both boundary regions, but directions
are opposite. The horizon is set to δ = 3.0 mm, and the bar is discretized with the grid
spacing ∆x = 1.5 mm. The material model defined for peridynamic bonds is shown in Figure 6.7(a), which demonstrates perfect plasticity at the bond level. The density of the bar
is 8000 kg/m3 . Setting Young’s modulus E = 200 GPa, we can calculate the micromodulus
c = 1.1968 × 1022 N/m6 . The computational step is set to 5 × 10−8 s, which is less than the
critical time step for the explicit time integration.
The macroscopic material behavior of the peridynamic bar is shown in Figure 6.7(b), in
which σx is the cross-sectional stress in the middle of the bar and εx is the strain imposed
by the boundary displacements. As indicated in Figure 6.7(b), the peridynamic bar yields
at the value very close to the yield stress 490 MPa, which is calculated by multiplying
Young’s modulus by the yield stretch. Unlike the material model at the bond level that
shows perfect plasticity after yielding, the cross-sectional stress σx at the macroscopic level
continues increasing till the value of around 930 MPa as shown in Figure 6.7(b). The reason
is that peridynamic bonds do not yield at the same time. After reaching the value at around
σx = 930 MPa, most bonds yield and the peridynamic bar demonstrates a small amount of
strain hardening as shown in Figure 6.7(c), and the cross-sectional stress starts to decrease
118

at the strain εx = 0.3.
Figure 6.7(d) shows the comparison of peridynamic material behavior at the macroscopic
level to the experimental result of Weldox 460 E steel presented in [21]. The Weldox 460
E steel has Young’s modulus E = 200 GPa and the yield stress of 490 MPa, which are
the same as the values that we use to define the material model for peridynamic bonds.
The peridynamic model demonstrates close agreement with the experimental results before
yielding, but the flow stress of the peridynamic model is larger than the experimental results
due to the scale of the constitutive model defined at the bond level for peridynamics. Strain
εx distributions by calculating the average bond stretches within the grid width ∆x as nodal
strains are plotted in Figure 6.8. As the displacements at boundary regions increase, bonds
in the peridynamic model first yield and then break as the bond stretch reaches the critical
value s0 = 0.5. Note that if the state-based peridynamics [156, 181, 58, 59] is used with a
proper constitutive model in the classical theory, better agreement with the experimental
results in the macroscopic material properties might be obtained.

119

Bond force (N/m6)

0.00245
Bond stretch

0.50000

σx (MPa)

(a)
1000
900
800
700
600
500
400
300
200
100
Peridynamics
0
0.00
0.01
0.02

0.03

0.04

0.05

εx

(b)
Figure 6.7: (a) Constitutive model defined for peridynamic bonds; (b) macroscopic material
response of the peridynamic bar (strain εx in the range of 0.00 to 0.05); (c) macroscopic
material response of the peridynamic bar; (d) comparison of peridynamic material behavior
at the macroscopic level to the experimental result of Weldox 460 E steel (strain rate in the
range of 0.00074 s−1 to 1522 s−1 ) presented in [21].

120

σx (MPa)

Figure 6.7 (cont’d)
1000
900
800
700
600
500
400
300
200
100
0
0.0

Peridynamics
0.2

0.4

0.6

0.8

1.0

εx

σ (MPa)
x

(c)
1000
900
800
700
600
500
400
300
200
100
0
0.0

Experiment
Peridynamics
0.1

0.2

0.3
εx

(d)

121

0.4

0.5

(a)

(b)

(c)

(d)

Figure 6.8: Strain εx contour: (a) boundary displacement ux = 0.12 mm; (b) boundary
displacement ux = 0.50 mm; (c) boundary displacement ux = 1.00 mm; (d) boundary
displacement ux = 2.50 mm.

122

projectile

rigid wall

y
x

Figure 6.9: Impact between two rigid bodies.

6.4.2

Rigid-body impact

In order to verify the contact algorithm, a rigid-body impact problem shown in Figure 6.9
is examined. The length of the cylindrical projectile is 80 mm, and the diameter is 20 mm.
The mass of the rigid projectile is 0.197 kg, and the initial velocity is vi = −120 m/s.
For the sake of detecting penetration and calculating the penetration length, the cylindrical
projectile is modeled with eight-node solid elements, and the bottom surface of the projectile
is comprised of 48 quadrilateral segments. The rigid wall is discretized into nodes with a
uniform Cartesian grid of 1.5 mm in width. The nodes in the rigid wall after discretization,
which serve as slave nodes, are constrained in all degrees of freedom during the calculation,
and the calculation time step dt is set to 5 × 10−8 s.
The displacement of the projectile is plotted in Figure 6.10(a). The maximum penetration
depth into the rigid wall is around 0.20 mm, and then the projectile is gradually bounced
up by the contact forces. As indicated in Figure 6.10(b), the final velocity of the projectile
after impact is v = 120.01 m/s, which is very close to the analytical value v = 120 m/s by
considering the conservation of kinetic energy. Hence, the proposed contact-impact procedure
123

is validated for this rigid-body impact problem.

124

0.15
Displacement (mm)

0.10
0.05
0.00
-0.05
-0.10
-0.15
-0.20
0

1

2

3

4

5

6

4

5

6

Time (µs)

(a)
150

Velocity (m/s)

100
50
0
-50
-100
-150
0

1

2

3
Time (µs)

(b)
Figure 6.10: (a) Displacement of the projectile and (b) velocity of the projectile.

125

6.4.3

Ballistic perforation

The ballistic penetration of a plate impacted by a blunt-nosed cylindrical projectile is examined for benchmarking. The blunt-nosed cylindrical projectile has the mass of 197 g,
diameter of 20 mm, and length of 80 mm. The plate dimensions are taken to be 60 mm
in length, 60 mm in width, and 12 mm in thickness. The density of the steel plate is
8000 kg/m3 , Young’s modulus is 200 GPa, and the yield stress is 490 MPa, which are the
same as material properties of Weldox 460 E steel. The ballistic penetration of Weldox 460
E steel is studied in [24], where a 12 mm thick plate is clamped in a rigid frame with the
inner clamp diameter of 500 mm. Due to the high computational cost, we cannot afford
modeling of the plate with the same in-plane dimensions as in the experiment in [24].

Figure 6.11: Modeling of a peridynamic plate impacted by a cylindrical projectile.

The target steel plate is modeled using peridynamics as shown in Figure 6.11. The
126

constitutive model at the bond level is defined as perfect plasticity such that the yield
stretch sy is 0.00245, and the critical stretch s0 is set to 0.5. The horizon is chosen to be
δ = 3.0 mm, and the grid spacing is ∆x = 1.5 mm. A finer grid spacing might be utilized, but
the corresponding computational cost will increase significantly. The projectile is modeled as
a rigid body with 48 eight-node solid elements as shown in Figure 6.11. In the simulations,
four lateral surfaces of the plate are constrained in x, y, and z directions. The computational
step is set to dt = 1 × 10−7 s.
Initial velocity vi
265
280
300
330
370
400

Residual velocity vr

m/s
m/s
m/s
m/s
m/s
m/s

33.12 m/s
68.67 m/s
114.97 m/s
161.47 m/s
232.75 m/s
254.84 m/s

Table 6.2: Numerical results of residual velocities.

Displacements of the projectile with different initial velocities are shown in Figure 6.12(a).
With the initial velocity vi = 250 m/s, the projectile indents the plate without perforation.
The maximum indentation depth is 6.2946 mm at the time t = 67 µs, and then the projectile
bounces back . The indentation formed by impact is shown in Figure 6.13, and a large part
of the plate is deformed by the indentation.
By performing a set of numerical experiments, the ballistic limit velocity is found to
be 261.4 m/s, which is larger than the value of 184.5 m/s determined by the experiments
presented in [24]. A few reasons can be considered herein. First, it has been found that the
ballistic limit velocity is sensitive to the mesh size in the finite element analysis [21], and
the numerically determined value of the ballistic limit velocity converges to the experimental
value as the number of elements over the plate thickness increases [21]. Therefore, we might
127

anticipate that the numerically determined ballistic limit velocity would be closer to the
experimental result if a smaller grid spacing ∆x is used, which, however, would incur excessive
computational cost. Second, as indicated in Figure 6.7, the bond-based peridynamic model
demonstrates a larger increase in the flow stress beyond the yield point than the experimental
result. To yield better agreement with the experimental result in the ballistic limit velocity,
further investigation of the material model using the state-based peridynamics is required
in the future. Third, the in-plane dimensions of the plate in the numerical simulations are
different from those in the experiments by [24].
For the perforations by blunt projectiles, the analytical model [145] has been proven to
be effective to estimate residual velocities. To compare numerical results with analytical
solutions of the projectile residual velocity, the approximated mass ratio a given in Equation (6.14) is calculated to be 0.87. By substituting a = 0.87 and vbl = 261.4 m/s into
Equation (6.13), analytical values of the residual velocity can be obtained. The numerically
determined residual velocities for the projectile with different initial velocities are summarized in Table 6.2, and the comparison of the numerical results to the analytical solutions
of the residual velocity is plotted Figure 6.12(b). Since the analytical model is derived by
considering the conservation of energy and momentum, the numerical results of residual velocities for the projectile with the initial velocities beyond the ballistic limit still demonstrate
good agreement with the analytical model, regardless of the higher ballistic limit velocity in
the numerical simulations.

128

Displacement (mm)

0
-5
-10
-15
vi=250 m/s
vi=265 m/s
vi=280 m/s
vi=300 m/s
vi=330 m/s
vi=370 m/s
vi=400 m/s

-20
-25
-30
0

20

40

60

80

100

Time (µs)

(a)

Residual velocity (m/s)

300
250
200
150
100
50

analytical
simulation

0

260 280 300 320 340 360 380 400 420
Initial velocity (m/s)

(b)
Figure 6.12: (a) Projectile displacement and (b) residual velocity.

129

(a)

(b)

Figure 6.13: Indentation formed by the projectile with an initial velocity vi = 250 m/s: (a)
bird’s-eye view of the plate; (b) lateral view of the plate.

130

The deformed target plates during perforation are plotted in Figures 6.14 and 6.15 for the
initial velocities vi = 300 m/s and vi = 400 m/s, respectively. The overall physical behavior
of the perforation process is captured by the numerical simulations. In the first phase, the
projectile indents the target, and the bending of the plate dominates. Then, the mass in
front of the projectile is accelerated by the contact forces. Subsequently, damages occur
near the shear zone surrounding the projectile after the peridynamic bond stretches reach
the critical value. In the last phase, a plug is formed after severe indentation. The shape
of the plug in numerical simulations is similar to the typical experimental observations. By
comparing Figures 6.14 and 6.15, it is observed that as the initial velocity increases, the
bulge on the rear side of the plate becomes localized, and the radius of deformation outside
of the indented area is reduced. This observation agrees with the investigation presented in
[24].

131

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6.14: Perforation of the plate by the projectile an the initial velocity vi = 300 m/s:
(a) t = 30 µs, bird’s-eye view of the plate; (b) t = 30 µs, lateral view of the plate; (c)
t = 70 µs, bird’s-eye view of the plate; (d) t = 70 µs, lateral view of the plate; (e) t = 120 µs,
bird’s-eye view of the plate; (f) t = 120 µs, lateral view of the plate.

132

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6.15: Perforation of the plate by the projectile an the initial velocity vi = 400 m/s:
(a) t = 30 µs, bird’s-eye view of the plate; (b) t = 30 µs, lateral view of the plate; (c)
t = 70 µs, bird’s-eye view of the plate; (d) t = 70 µs, lateral view of the plate; (e) t = 120 µs,
bird’s-eye view of the plate; (f) t = 120 µs, lateral view of the plate.

133

6.5

Summary

In this chapter, a numerical scheme is proposed for the modeling of contact of a peridynamic
domain and a non-peridynamic domain (e.g. conventional finite elements and rigid bodies)
under high velocities. Starting with the peridynamic theory, an approach for the threedimensional inverse isoparametric mapping based on Taylor’s series is then summarized.
The inverse isoparametric mapping is effective to determine whether a slave node is within a
solid element by determining the natural coordinates of the slave node. A contact algorithm
for the calculation of the contact point and the penetration depth is then explained. For
transient analyses by the explicit time integration, a penalty method is employed to enforce
displacement constraints on the contact surfaces. The contact forces applied on the slave
node and the master segment are proportional to the penetration depth.
The contact-impact problems are investigated numerically. First, the material behavior
of a three-dimensional peridynamic model is examined. The results indicate that the macroscopic material behavior of a peridynamic bar is different from the material model defined
at the bond level for peridynamics. Then the impact of a rigid projectile into a rigid wall
is studied. By comparing the kinetic energy before and after impact, the proposed numerical scheme for the contact-impact procedure between two rigid bodies is validated. The
last numerical example is the perforation of a steel plate by a blunt-nosed cylindrical rigid
projectile. The plate, where severe damage develops during perforation, is modeled using
peridynamics. Indentation is formed if the initial velocity of the projectile is smaller than
the ballistic limit velocity, and the value of ballistic limit velocity is determined by running
a set of numerical experiments. Compared with the ballistic limit velocity to perforate the
plate investigated by experiments in [24], the numerically determined value of the ballistic
134

limit velocity is larger due to the differences in the macroscopic material properties and the
geometry of the target plate. For the projectile with an initial velocity beyond the ballistic
limit velocity, a comparison of numerical results to the analytical solutions of the residual
velocity is conducted, and good agreement is observed. The overall physical behavior of
the perforation process is captured by the numerical simulations employing the proposed
contact-impact procedure.

135

Chapter7
Modeling of Porous Brittle Solids using Peridynamics
7.1

Introduction

Heterogeneous materials consist of more than one materials, insoluble in one another [26].
A finite and heterogeneous solid made of m subdomains is shown in Figure 7.1. Heterogeneous materials have been widely used throughout human history from concrete, the earliest
application of which was by Romans [119], to composites used by the aircraft industry.
Typical heterogeneous materials include, for example, metal alloys, polycrystalline, composites, and porous media. The macroscopic material behaviors of heterogeneous materials are
determined by the size, shape, properties, and spatial arrangement of constituents [96].

Ωm

Ωn

Ω1

Ω2

Figure 7.1: Finite and heterogeneous solid [169].

As a type of heterogeneous material, porous materials comprise a solid phase and pores
136

[105]. It has been found that the porosity has a significant effect on the elastic modulus of
materials. Spriggs [163] proposed an empirical equation describing the effect of porosity on
the elastic modulus of solids. Wang [177, 178] obtained a theoretical relationship between
the porosity and Young’s modulus by considering a cubic stacking pattern and Young’s
modulus of porous alumina with different pore structures. Phani and Niyogi [137] derived a
semi-empirical equation for the porosity dependence of Young’s modulus of brittle solids.
Similarly, the porosity also has significant effects on the strength of materials. In a commentary on the investigations by Ryshkewitch [152], Duckworth [43] suggested an expression
for the influence of porosity on the strength. Considering the combined effect of porosity
and grain size, Knudsen [95] developed an empirical equation. Hasselman [74] proposed
relations between effects of porosity on the strength and Young’s modulus of elasticity of
polycrystalline materials. Vardar et al. [175] investigated the effect of spherical pores on the
strength of a ceramic material, and compared the experimental results of brittle strength
and the predictions using the Weibull probabilistic approach. Krstic [99, 100] proposed an
analytical model describing the strength degradation of brittle solids containing spherical
pores and annular flaws, and developed a unified approach to determine the effect of microstructure on fracture of brittle materials. An evaluation of strength-porosity relationships
was presented by Dutta [45], and a semi-empirical relationship for the assessment of strength
by Young’s modulus was proposed. Wagh et al. [176] proposed an analytical model incorporating connected-grain model to describe the dependence of flexural strength of ceramics on
the porosity. Nielsen [129] studied relations between the stiffness and the strength of porous
materials considering the pore geometry. Kearsley and Wainwright [89] studied the effect
of porosity on the strength of foamed concrete. R¨ssler and Odler [149] investigated the
o
relationships between the porosity and the strength for a series of cement paste specimens.
137

Kumar and Bhattacharjee [102] experimentally studied the in situ concrete, and presented
an empirical model relating in situ strength of concrete with the porosity, pore size distribution, and binder content. Weiler et al. [184] used the X-ray tomography to examine the
magnesium alloy, and found that the tensile properties are mainly determined by the local
areal fraction of porosity.
In many engineering applications, it is necessary to determine the macroscopic characteristics of heterogeneous materials. However, the complexities in the material properties
and microstructural characteristics make it difficult to study heterogeneous materials. The
simplest method is to homogenize properties of a heterogeneous material as an average over
the properties of the constituents [96]. With the development of computational methods, research progresses for the analyses of heterogeneous materials have been made. For example,
Sumi and Wang [169] used the finite element method to study the growth of non-collinear
cracks in a two-dimensional heterogeneous solid. Due to the length scales involved in heterogeneous materials, it is not possible to generate meshes that can accurately represent the
microstructure. Therefore, multi-scale modelings are promising alternative approaches. For
example, the spatially periodic representative volume element has been applied to study the
large-strain mechanical responses of voided polycarbonate [162]. Kouznetsova et al. [96, 97]
presented a micro-macro strategy for the modeling of heterogeneous materials at large deformations.
In this chapter, peridynamics is employed to study brittle materials containing pores.
First, analytical and empirical equations describing the effect of porosity on Young’s modulus and strength of solid materials are summarized. A few numerical examples are then
presented. Specimens containing randomly distributed cubic voids with varying amounts
of porosity are studied. The degradation of Young’s modulus and strength with increasing
138

porosity in the numerical results is compared with analytical and empirical models. Brittle solids containing spherical voids are studied for different porosities, and a comparison
between numerical results and analytical solutions is conducted.

7.2

Theory

For porous brittle solids, the porosity P can be described using the volume fraction porosity
as
P =

Void volume
.
Total volume

(7.1)

The empirical equation, proposed by Spriggs [163], for the prediction of Young’s modulus at
the porosity value of P has been widely used as

E = E0 e−bP ,

(7.2)

where E0 is the zero-porosity Young’s modulus, and b is an empirical constant. Considering
the densification of spherical particles in a simple cubic array, Wang [177, 178] derived the
theoretical solution of Young’s modulus of porous materials. The relative density X, which
equals 1 − P , is described by [177]

X(θ) =







π
12

9 −3− 4
cos2 θ
cos3 θ



 (tan2 θ − 1)1/2 +


,

for θ

45◦

√
sin−1 ( 2cosθ)
1
[sin3 φ(π − 4α + 2sin2α)]dφ,
4cos3 θ θ

,
for θ

45◦
(7.3)

where θ is the angle of coalescence, and α = cos−1 (cosθ/sinφ). θ = 0◦ if the sphere is
inscribed within the unit cell, and θ = 54.74◦ for the 100% density. For the ideal case in
139

which the load is applied on a perfect cubic array, the normalized Young’s modulus is given
as [177]






E
(θ) =

E0




−1
dφ
1 ln tan(π/4−θ/2) + π/2
, for θ 45◦
π
tan(θ/2)
π/2−θ sinφ(π−4α+2sin2α)
−1
−1 √
dφ
2 θ − 1)1/2 + 4cosθ sin ( 2cosθ)
, for θ
(tan
θ
sinφ(π−4α+2sin2α)

1
4cosθ

.
45◦
(7.4)

Figure 7.2(a) shows the value of the relative density that is a function of the angle of coalescence, and Figure 7.2(b) shows Young’s modulus of porous materials as a function of the
relative density by numerically evaluating Equations (7.3) and (7.4), respectively. Based on
the theoretical solution, Wang [177] proposed an approximated solution as

2
E = E0 e−(bP +cP ) ,

(7.5)

where b and c are material constants.
The dependence of strength on the porosity has been represented by a few empirical
equations. Based on experimental data by Ryshkewitch [152], Duckworth [43] described the
relationship as
σf = σ0 e−bP .

(7.6)

where σ0 is the zero-porosity strength, and b is a constant. Dutta et al. [45] proposed an
empirical equation as
σf = σ0 (1 − P )m ,

(7.7)

where m is a constant. Krstic [99] derived the strength of a solid containing spherical voids

140

Relative density

1.0
0.8
0.6
0.4
0.2
0.0
0

10

20

30
θ (°)

40

50

60

(a)

E/E0

1.00

0.10

0.01
1.0

0.9

0.8
0.7
Relative density

0.6

0.5

(b)
Figure 7.2: (a) Relative density versus angle of coalescence and (b) theoretical solution of
Young’s modulus of porous materials.

141

as

σf =

πγE0
D(1 − ν 2 )(1 + s/R)
(4 − 5ν)
4P (1 − ν 2 )
3
+
× 1+
× 2(1 + s/R)3 +
2
π
2(7 − 5ν)
2(7 − 5ν)(1 + s/R)

−1/2

,

(7.8)
where γ is the fracture surface energy of a crack which equals Gf /2 for ideally brittle materials [5], D and R are the diameter and radius of spherical pores respectively, and s is the
annular flaw size as shown in Figure 7.3.

s D s

Figure 7.3: Spherical pores with annular flaws [99].

7.3

Numerical implementation

For numerical simulations, voids are generated on the discretization grid of a peridynamic
domain. First, the Fortran intrinsic function random_seed is used to initialize the random
seed based on the system’s time. Then the Fortran intrinsic function random_number is
142

invoked to return an array consisting three numbers from the uniform distribution over the
range of 0 to 1. By multiplying the number of discretization nodes in x, y, z directions and
rounding to integers, we can obtain the position to generate a void in the discretization grid.
A spherical void can be generated by deleting peridynamic nodes within the spherical radius
with the center determined by random numbers. A number of iterations are performed till
the volume fraction is close to the desired value of porosity within a tolerance. The algorithm
to generate spherical voids is summarized in Table 7.1.
1.
2.
3.
4.
5.
6.

Discretize the peridynamic domain.
Call random_seed to initialize the random seed.
Call random_number to generate three random numbers in the range of 0 to 1.
Determine the center to generate a spherical void in the discretization grid.
Determine the distance to the nearest void and specimen boundaries.
Go to step 3 if the distance determined in step 5 is less than the predefined
value.
7. Delete all nodes within the radius of the spherical voids.
8. Calculate the volume fraction.
9. Go to step 3 to generate additional spherical pores if the volume fraction is
less than the desired value of porosity.
Table 7.1: Algorithm to generate randomly distributed spherical voids.

7.4
7.4.1

Case studies
Cubic voids

Peridynamics is employed to study material behaviors of a porous brittle material. The
material properties of the non-porous body are taken to be Young’s modulus E0 = 70 GPa
and the critical bond stretch s0 = 0.05. Consider a three-dimensional model of dimensions
25.5 mm × 25.5 mm× 25.5 mm. The size of the horizon is set to δ = 1.0 mm, and the
143

peridynamic model is discretized with the grid spacing ∆x = 0.5 mm, which ends up with
51×51×51 peridynamic nodes. The horizon to the grid spacing ratio is selected to be 2 since
it is good to match the macroscopic material behavior in the elastic range. The numerical
simulations are conducted on a workstation with an NVIDIA C1060 GPU. A constant velocity of 10 mm/s is applied on boundary regions at both ends of the specimen for tensile
loading, and the width of each boundary region is 2.5 mm. For the numerical calculation,
the time step dt is set to be 5 × 10−8 s.
In order to introduce cubic voids into the peridynamic model, the random number generator is utilized to select a number of peridynamic nodes so that the volume fraction of
selected nodes equal to the desired porosity as described by Equation (7.1). The flowchart
is similar to the algorithm described in Table 7.1. Voids can be introduced by either deleting the selected peridynamic nodes or deleting all peridynamic bonds connecting to those
nodes. In the algorithm to generate randomly distributed cubic voids, it is prescribed that
no adjacent peridynamic nodes are deleted simultaneously. Therefore, only cubic voids are
generated in the peridynamic model. Figure 7.4 shows the specimen that contains 5% cubic
voids.
Porosity

Young’s modulus (GPa)

Tensile strength (MPa)

0%
5%
10%
15%
20%

71.7742
67.1853
62.1410
57.7252
52.7812

281.3537
236.0355
210.6670
188.7893
165.6249

Table 7.2: Young’s modulus and tensile strength of specimens containing cubic voids.

The cross-sectional stresses are plotted in Figure 7.5 for peridynamic models with varying
amounts of porosity. For the nonporous model, the resultant Young’s modulus by calculating
144

Figure 7.4: Cubic voids in a specimen of dimensions of 25 mm×25 mm×25 mm. Porosity
P = 5% .
the tangent of the stress-strain curve is E = 71.87 GPa, which is close the exact Young’s
modulus of the material. The strength of the nonporous model is σf = 280.97 MPa. As the
porosity P increases, both Young’s modulus and strength decrease as shown in Figure 7.5.
Young’s moduli and strengths for different porosities are summarized in Table 7.2. The effect
of porosity on Young’s modulus is plotted in Figure 7.6(a). For the comparison, the empirical
equation given in Equation (7.2), the approximate solution presented in Equation (7.5), and
the theoretical solution given in Equation (7.4) are utilized. As indicated in Figure 7.6(a),
the numerical results are in good agreement with empirical, approximate, and theoretical
solutions. Similarly, the effect of porosity on the strength is plotted in Figure 7.6(b). The
numerical results are compared with the solutions by fitting empirical models given in Equations (7.6) and (7.7) to the numerical data. Close agreement between the numerical results
and empirical models is observed in Figure 7.6(b).

145

400

σx (MPa)

300
200

P= 0%
P= 5%
P=10%
P=15%
P=20%

100
0
-100
0.000

0.001

0.002

0.003

0.004

0.005

εx
Figure 7.5: Cross-sectional stress versus strain of specimens containing cubic voids.

146

1.0

E/E0

0.9
0.8
0.7
Peridynamics
Spriggs (empirical)
Wang (approximate)
Wang (theoretical)

0.6
0.5
0.4
0.00

0.05

0.10
Porosity

0.15

0.20

(a)

Peridynamics
Dutta et al.
Duckworth

1.0

σf/σ0

0.9
0.8
0.7
0.6
0.5
0.00

0.05

0.10
Porosity

0.15

0.20

(b)
Figure 7.6: (a) Normalized Young’s modulus and (b) normalized tensile strength of specimens
containing cubic pores.

147

7.4.2

Spherical voids

In the following, peridynamic models that contain spherical pores are examined for different
porosities. The dimensions of the peridynamic model are 55 mm in length, 35 mm in width,
and 35 mm in thickness. The grid spacing ∆x is 0.5 mm, and the horizon δ is set to be 2∆x.
Spherical pores are randomly generated within the boundary of peridynamic models, and
specimens containing 5%, 10%, 15%, and 20% porosities are demonstrated in Figure 7.7.
The diameter of spherical pores is set to D = 11∆x. Young’s modulus of nonporous body
is E0 = 70 Gpa, and Poisson ratio is 0.25. For the computation, time step dt is set to be
5 × 10−8 s. A constant velocity of 10 mm/s is applied on boundary regions at both ends of
specimens for tensile loading, and the width of each boundary region is 2.5 mm.
Figure 7.8 shows the cross-sectional stress versus strain imposed by the boundary displacement for models of different porosities. The resultant Young’s modulus of the nonporous model is 71.65 GPa, and the strength is σf = 266.56 MPa. As the porosity increases,
both Young’s modulus and strength decrease as shown in Figure 7.8. Young’s moduli and
strengths of peridynamic models for different porosities are summarized in Table 7.3. In Figure 7.9(a), the normalized Young’s modulus is compared with the theoretical solution given
in Equation (7.4) and the fitted curves by Spriggs’ empirical relation given in Equation (7.2)
and Wang’s approximated solution given in Equation (7.5). As indicated in Figure 7.9(a),
Young’s moduli of peridynamic models for different porosities show close agreement to the
theoretical solutions and the fitted curves using empirical and approximate models. The
tensile strength of peridynamic models containing spherical pores is compared to the theoretical solution given in Equation (7.8). Note that a flaw size of s is assumed at pore edges
in Equation (7.8). Theoretical solutions considering flaw sizes of 0, ∆x and 2∆x are plotted
148

in Figure 7.9(b). Since the horizon δ = 2∆x is utilized for the peridynamic model, we might
anticipate that s = ∆x might yield the best agreement to the theoretical solution, which is
confirmed in Figure 7.9(b).
Porosity

Young’s modulus (GPa)

Tensile strength (MPa)

0%
5%
10%
15%
20%

71.6015
66.0444
63.1007
57.8934
54.2525

272.1414
167.2611
159.1419
145.0795
138.0373

Table 7.3: Young’s modulus and tensile strength of specimens containing spherical voids
within the boundary.

Similar to the case that spherical voids are generated within the boundary as shown in
Figure 7.7, the spherical voids intersecting with boundaries can also be generated as shown
in Figure 7.10. The stress-strain curves with different porosities in the peridynamic model
are plotted in Figure 7.11, which indicates that both Young’s modulus and strength decrease
as the porosity increases. The resultant Young’s modulus and strength are summarized in
Table 7.4. Figure 7.12(a) shows the comparison of the numerical results to the analytical and
empirical models in the degradation of Young’s modulus, and good agreement is observed.
The trend of decreasing strength with increasing porosity is plotted in Figure 7.12(b). The
numerical results are compared with the theoretical values by solving Equation (7.8). As
demonstrated in Figure 7.12(b), the numerical results are close to the theoretical values
considering the annular flaw size in the range of s = 1∆x to s = 2∆x.
By comparing Tables 7.2 and 7.3 which summarize Young’s moduli and tensile strengths
for solids containing cubic and spherical voids respectively, the effect of the shape of voids on
material properties can be found. The difference in Young’s modulus of models containing
cubic and spherical pores is inconsiderable. On the other hand, tensile strengths of models
149

Porosity

Young’s modulus (GPa)

Tensile strength (MPa)

0%
5%
10%
15%
20%

71.6015
66.8800
62.9406
59.0915
54.3961

272.1414
164.3560
138.5337
136.4608
129.4621

Table 7.4: Young’s modulus and tensile strength of specimens containing spherical voids
intersecting with boundaries.
containing spherical pores decrease. Similarly, by comparing Tables 7.3 and 7.4, we can
notice that the difference in Young’s modulus of specimens containing spherical pores intersecting with boundaries is insignificant compared to the results of spherical pores placed
within the boundaries. On the other hand, the strength decreases if spherical pores having
intersections with the boundaries as indicated by comparing Figures 7.9(b) and 7.12(b).

150

(a)

(b)

(c)

(d)
Figure 7.7: Spherical pore distributions in the specimen and distributions at the cross section
along the longitudinal direction. Porosities (a) P = 5%, (b) P = 10% , (c) P = 15%, and
(d) P = 20%. Pores are generated within the specimen boundaries.
151

400

σx (MPa)

300
200

0% porosity
5% porosity
10% porosity
15% porosity
20% porosity

100
0
-100
0.000

0.001

0.002

0.003

0.004

0.005

εx
Figure 7.8: Cross-sectional stress versus strain of specimens containing spherical pores within
the boundary.

152

1.0

E/E0

0.9
0.8
0.7
Peridynamics
Spriggs (empirical)
Wang (approximate)
Wang (theoretical)

0.6
0.5
0.4
0.00

0.05

0.10
Porosity

0.15

0.20

(a)

σf (MPa)

200
150
100
peridynamics
s=0∆x
s=1∆x
s=2∆x

50
0
0.00

0.05

0.10
Porosity

0.15

0.20

(b)
Figure 7.9: (a) Normalized Young’s modulus and (b) tensile strength of specimens containing
spherical pores within the boundary.

153

(a)

(b)

(c)

(d)

Figure 7.10: Spherical pore distributions in the specimen. Porosities (a) P = 5%, (b)
P = 10% , (c) P = 15%, and (d) P = 20%. Pores intersect with specimen boundaries.

400

σx (MPa)

300
200

0% porosity
5% porosity
10% porosity
15% porosity
15% porosity

100
0
-100
0.000

0.001

0.002

0.003

0.004

0.005

εx
Figure 7.11: Cross-sectional stress versus strain of specimens containing spherical pores
intersecting with boundaries.
154

1.0

E/E0

0.9
0.8
0.7
Peridynamics
Spriggs (empirical)
Wang (approximate)
Wang (theoretical)

0.6
0.5
0.4
0.00

0.05

0.10
Porosity

0.15

0.20

(a)

σf (MPa)

200
150
100
peridynamics
s=0∆x
s=1∆x
s=2∆x

50
0
0.00

0.05

0.10
Porosity

0.15

0.20

(b)
Figure 7.12: (a) Normalized Young’s modulus and (b) tensile strength of specimens containing spherical pores intersecting with boundaries.

155

7.5

Summary

In this chapter, peridynamics is employed to study brittle materials containing pores. Starting with analytical and empirical equations describing the effect of porosity on Young’s
modulus and the strength of solid materials, an algorithm is then introduced for generating
cubic and spherical voids in discretized peridynamic models using the pseudorandom number
generator.
Several numerical examples are given next. Specimens containing randomly distributed
cubic voids with different porosities are studied. As the value of porosity increases, Young’s
modulus decreases, and the trend of degradation shows close agreement with the empirical
models and the analytical solutions derived by considering the simple cubic stacking pattern.
The degradation of strength with increasing porosity conforms to the empirical models.
Brittle solids containing spherical voids are then studied with varying amounts of porosity.
Two cases, spherical pores distributed inside the specimen boundaries and pores having
intersections with the boundaries, are examined. A comparison of numerical results to the
empirical models and the analytical solutions in the degradation of Young’s modulus is
conducted, and the differences are within a few percent. The degradation of strength in the
solids containing spherical pores is compared with the analytical solutions considering the
effect of pore size and pore volume fraction. It is found that numerical solutions are close
to analytical values if the size of the horizon is considered as the annular flaw extending
from the void surface. In general, peridynamics is successfully implemented to study porous
brittle solids.

156

Chapter8
Conclusions and Recommendations
8.1

Conclusions

The primary focus of this research is development of the discretized peridynamics for solid
mechanics. The thrust of research on peridynamics provides an alternative theory to the
classical continuum mechanics that is directly applicable for numerical simulations of spontaneously formed discontinuities. The merit of peridynamics comes from the integration
used for the calculation of forces on material points. This is in contrast to the classical
theory in which partial derivatives are involved in the governing equation, and consequently
approaches based on the framework of classical continuum mechanics fail to be directly
applicable in the presence of singularities.
Chapter 3, Chapter 4, Chapter 5, Chapter 6, and Chapter 7, each has a different area
within the framework of bond-based peridynamics. First, a connection between the classical elasticity and discretized peridynamics is established. By introducing the peridynamic
stress corresponding to the classical (local) stress, calculating the corresponding peridynamic
Young’s modulus, and equilibrating the peridynamic Young’s modulus with the conventional
Young’s modulus, micromoduli for one- and three-dimensional models are derived. The micromoduli derived in this research are referred to as numerical micromoduli since the value
of micromodulus for a given horizon size is determined by numerical calculations on a dis157

cretization grid. The numerical micromoduli are different from the analytical micromoduli
presented in [20, 49, 50] in that a finite number of peridynamic bonds are considered within
the horizon of a node in the numerical micromoduli. It is found that the one-dimensional
analytical micromoduli can yield the numerical micromoduli by substituting discretization
parameters into the analytical micromoduli. For three-dimensional models, numerical micromoduli converge to the analytical micromoduli as the ratio of the horizon to the grid
spacing increases. In the bond-based peridynamic theory, Poisson ratio ν is a fixed value of
1/4. For numerical simulations of an elastic body which has Poisson ratio not equal to 1/4,
a pairwise compensation scheme is introduced in this research. Compensation forces ˆ are
f
superposed in the directions perpendicular to the direction of the bond force f , resulting in
additional transverse deformations to simulate the effect of Poisson ratio other than 1/4. The
numerical results applying the proposed pairwise compensation scheme match the classical
(local) solutions in the material responses including the stress, the strain, back-calculated
Young’s modulus, and Poisson ratio.
The research-purpose peridynamics code is implemented in the NVIDIA graphics processing unit (GPU) Tesla C1060 for the highly parallel computation. To the best of our
knowledge, it was the first implementation of peridynamics on a GPU. An algorithm separating nested loops into independent loops is introduced to eliminate the loop dependency
in the peridynamics code for the efficient parallelization. For three-dimensional problems, a
significant speedup over the serial calculation on a CPU is achieved. With the high computational efficiency of GPU, numerical studies are conducted to investigate the responses of
brittle and ductile material models. Stress-strain behaviors with different grid sizes and horizons are studied for a brittle material model. A comparison of stresses and strains between
finite element analyses and peridynamic solutions is performed for a ductile material. To
158

bridge material models at different scales, a multiscale procedure is proposed. First, material
properties are retrieved from the macroscopic responses of the peridynamic model. These
retrieved properties are then used to modify the material model for finite element analyses.
By applying the proposed procedure, the differences between finite element analyses and
peridynamic solutions reduce substantially.
An approach to couple the discretized peridynamics and FEM is proposed to take advantage of the generality of peridynamics and the computational efficiency of FEM. The coupling
of peridynamic and finite element subregions is achieved by means of interface elements. The
proposed method in this research is different from the coupling method implementing peridynamics in a conventional finite element analysis code [115], the approach using overlapping
regions [94], the submodeling approach [2, 133], and the morphing strategy [113]. Depending on how coupling forces are subdivided to FE nodes in an interface element, two types
of coupling schemes, the VL-coupling scheme and the CT-coupling scheme respectively, are
discussed. In the VL-coupling scheme, coupling forces are subdivided among all FE nodes
comprising the interface elements. In the CT-coupling scheme, interfaces between peridynamic subregions and finite element subregions are defined in the reference configuration,
and coupling forces are subdivided only among FE nodes at interfaces. For one-dimensional
simulations, it is found that the CT-coupling scheme yields the solution very close to the
analytical solution in all the domain including the finite element subregion, peridynamic subregion, and interfaces. For three-dimensional simulations, it is observed that the CT-coupling
scheme is effective in resolving displacements normal to the interface of peridynamic and FE
subregions. On the other hand, the VL-coupling scheme is capable to preserve Poisson effect
at the interfaces. The proposed coupling approach is used to model the mixed mode fracture
in a double-edge-notched concrete specimen, and numerical predictions of crack paths show
159

good agreement with the experimental result of crack patterns presented in [130].
A numerical scheme for the contact-impact procedure ensuring compatibility between a
peridynamic domain and a non-peridynamic domain is developed. The contact between a
discretized peridynamic domain and a non-peridynamic domain is modeled as a node-tosurface type. A penalty method is used to enforce displacement constraints, and dynamic
analyses are performed by an explicit algorithm without iterations, which is different from
the contact method used in [115, 107]. In the numerical examples, the impact between two
rigid bodies is presented to validate the contact algorithm. The ballistic perforation of a
steel plate is investigated numerically. Good agreement between the numerical simulations
and the analytical model is found in the results of residual velocities. The physical process
of perforation is well captured in the simulations using the proposed contact-impact scheme.
Peridynamics is applied to study porous brittle materials. An algorithm is developed
to generate randomly distributed cubic voids and spherical voids for a given porosity. The
material behaviors at the macroscopic level including the resultant Young’s modulus and the
strength are studied for varying amounts of porosity. The comparison of numerical results to
the empirical models and the analytical solutions in the degradation of Young’s modulus is
conducted, and the differences are within a few percent. The degradation of strength in the
solids containing spherical pores is compared with the analytical solution in which the effect
of pore size and pore volume fraction are considered. It is found that the numerical solution
is close to the analytical value if the horizon size is considered as the annular flaw extending
from the void surface. In general, peridynamics is successfully implemented to study porous
brittle solids.

160

8.2

Contributions

The contributions in this dissertation are summarized in the following:
• Connection between the classical elasticity and discretized peridynamics is established
in terms of peridynamic stress.
• Numerical micromoduli for one- and three-dimensional discretized models are derived.
• A pairwise compensation scheme is introduced for numerical simulations of materials
with Poisson ratio not equal to 1/4.
• The research-purpose peridynamics code is implemented in a GPU for highly parallel
computation, and a significant speedup is achieved.
• A multiscale procedure is proposed to bridge material models at different scales.
• An approach to couple the discretized peridynamics and the finite element method is
proposed, and two types of coupling schemes are investigated.
• The mixed-mode crack pattern is investigated using the proposed coupling approach.
• A numerical scheme for the contact-impact procedure ensuring compatibility between
a peridynamic domain and a non-peridynamic domain is developed.
• The physical process of perforation is well captured using the proposed contact-impact
scheme.
• Peridynamics is employed to study porous brittle materials.
161

8.3

Recommendations for future research

Compared with the classical theory, peridynamics is a relatively new development. While this
research has served to advance the development of peridynamics in establishing a connection
between the discretized peridynamic models and classical continuum mechanics, the high
performance computation of peridynamics, bridging material models at the peridynamic
bond level and the macroscopic level for finite element analyses, coupling peridynamics with
the conventional finite element method, the contact-impact procedure between a peridynamic
domain and a non-peridynamic domain, and applications on porous brittle materials, there
is considerable further research to be undertaken. The following are some aspects for further
research of peridynamics.
First, the focus in this research is the bond-based peridynamics. Since constitutive models
are defined at the bond level, the macroscopic material behavior might be different from the
material behavior at the bond level especially for plastic materials. In Chapter 4, a multiscale
approach is proposed to define a constitutive model for the finite element analysis based on
the macroscopic material behavior of a peridynamic model. But it is difficult to recast a
material model at the bond level based on the constitutive model in the classical theory, which
is dependent on the stress tensor. A generalized formulation of the bond-based peridynamics
is proposed by Silling et al. [156], which is referred to as the state-based peridynamics. A
force state T, which is similar to the stress tensor, is proposed, and constitutive models in
the classical theory are more likely to be implemented in the state-based peridynamic theory.
The force density acting on the material point x is defined by T − T′ , where T′ is the force
state on the material point x′ . Note that the direction of the force state T might be different
from the direction of the force state T′ in the state-based peridynamic theory as shown in
162

Figure 8.1.

T’
T
x

x’

ξ

Figure 8.1: State-based peridynamics [156].

The coupling approach of peridynamics with the finite element method has been proven
to be effective in Chapter 5. The numerical examples are validated under the quasi-static
condition. It is interesting to study coupling models subjected to the dynamic loading. Since
the stress waves have influence on the fracture behavior, particular attention is deserved
on the stress waves at interfaces of coupling models. Furthermore, the coupling approach
presented in Chapter 5 and the contact-impact procedure presented in Chapter 6 can be
combined together to further reduce the computational burden. For example, a hybrid
model discretized by peridynamic regions and finite element regions subjected to the impact
of another body can further take advantage of the generality of peridynamics in the presence
of continuities and the efficiency of FEM.
Most research efforts in the area of peridynamics have been concentrated on metals.
In Chapter 5, the mixed mode fracture in a concrete specimen is studied, and numerical
predictions of crack patterns are compared with experimental results. In the peridynamic
163

model, the concrete specimen is modeled using an ideally brittle material model. With the
framework of porous materials presented in Chapter 7, a more realistic model of concrete
using peridynamics can be developed. In the modeling of porous materials, randomly distributed voids are generated by deleting peridynamic bonds. Instead, voids can be treated
as aggregates by introducing corresponding material properties, and bonds acting on voids
can be considered as interactions between aggregates and cement if weak bond forces are
introduced accordingly.

164

APPENDICES

165

AppendixA
Derivation of Three-dimensional Micromodulus
Applying the volume reduction scheme, we have

VJ = βJ VJ = βJ (∆x)3 .

(A.1)

The peridynamic bond force is a function of the bond stretch s as

f = c3 s.

(A.2)

The projection of bond force f into the x axis is given as

fx = f

|xJ − xI |
|x − xI |
= c3 s J
.
ξ
ξ

166

(A.3)

The peridynamic stress is expressed as

1
σx =
AI
=

=

=

1
AI

NL
J=1
NL

(fx βJ VJ )VI
|x − xI |
c3 s J
βJ VJ
ξ

J=1
NL

1
(∆x)2
NL
J=1

VI
(A.4)

|x − xI |
βJ (∆x)3 (∆x)3
c3 s J
ξ

J=1

|x − xI |
c3 s J
βJ
ξ

(∆x)4 .

Since the elastic body is subjected to the isotropic expansion, the stresses in y− and
z−directions are identical to the x−directional stress σx , we have

σy = σz = σx .

(A.5)

By applying Hooke’s law of linear elasticity, the peridynamic Young’s modulus Epd can be
obtained as
σx ν(σy + σz )
−
s
s
σx 2σx
=
−
s
4s
σx
=
2s 

Epd =

NL



|x − xI |
1 
c3 s J
βJ (∆x)4 
2s
ξ
J=1


NL
c
|xJ − xI |
= 3
βJ (∆x)4  .
2
ξ

=

J=1

167

(A.6)

Setting peridynamic Young’s modulus equal to the exact Young’s modulus of the material,
we have
Epd = E.

(A.7)

The three-dimensional micromodulus c3 is obtained as

c3 = 


2E
NL
J=1

|xJ − xI |
βJ
ξ

168

.

(∆x)4 

(A.8)

AppendixB
Derivation of Modified Peridynamic Young’s
Modulus
ˆ
The modified peridynamic Young’s modulus Epd should be determined based on Epd to
keep the strain identical. The longitudinal strain in x−direction in Figure 3.11 equals the
strain in Figure 3.10, which is written as

σx
1
σx
−
=
ˆ
Epd
Epd 4

σy
ˆ
σz
ˆ
+
ˆ
ˆ
Epd Epd

.

(B.1)

The strains in y− and z−directions also need to be identical to the strains in the original
peridynamic model shown in Figure 3.10. Therefore, the equilibrium of the strain εy is
expressed as
−ν

σy
ˆ
σx
1
=
−
ˆpd 4
Epd
E

σx
σz
ˆ
+
ˆpd Epd
ˆ
E

.

(B.2)

Since the compensation forces are identical in all lateral directions perpendicular to the
pairwise force f , the lateral stresses are equal as

σy = σz .
ˆ
ˆ

169

(B.3)

Substituting Equation (B.3) into Equation (B.1), we have
ˆ
σx
1 σy
σx
−
.
=
ˆ
ˆ
Epd
Epd 2 Epd

(B.4)

ˆ
By multiplying Epd on both sides, Equation (B.4) is written as
ˆ σx = σx − 1 σy .
Epd
ˆ
Epd
2

(B.5)

Substituting Equation (B.3) into Equation (B.2), we have

−ν

ˆ
σx
1 σx
3 σy
−
.
=
ˆ
ˆ
Epd
4 Epd 4 Epd

(B.6)

ˆ
By multiplying Epd on both sides, Equation (B.6) is written as

ˆ
−Epd ν

σx
1
3
ˆ
= σy − σx .
Epd
4
4

(B.7)

Substituting Equation (B.5) into Equation (B.7), we have

1
−ν σx − σy
ˆ
2

3
1
= σy − σx .
ˆ
4
4

(B.8)

The transverse lateral stress is expressed as

σy =
ˆ

1 − 4ν
σx .
3 − 2ν

170

(B.9)

By substituting Equation (B.9) and Equation (B.3) into Equation (B.1), the modified peridynamic Young’s modulus is obtained as

ˆ
Epd =

5
E .
6 − 4ν pd

171

(B.10)

AppendixC
Nonlinear Least Squares Method
Consider a data set consisting of n points (xi , yi ), i = 1, 2, · · · , n. The model function is
given as f (x, β), where β consists of m adjustable parameters. A residual [1] is defined as

ri = yi − f (xi , β),

(C.1)

and the least squares method is to the find the minimum of
n

ri 2 .

S=

(C.2)

i=1

Since there is no analytical solution to a nonlinear least squares problem, the solution of the
adjustable parameter β needs to be determined by successive iterations as

βj k+1 = βj k + ∆βj .

172

(C.3)

At each iteration, the model function can be linearized using the first-order Taylor’s series
about β k as
m

f (xi , β) = f k (xi , β) +
j=1
m

= f k (xi , β) +

∂f (xi , β)
(βj − βj k )
∂βj
(C.4)

Jij ∆βj ,
j=1

where Jij is the Jacobian matrix. The residual given in Equation (C.1) can be written as
m

ri = yi − f k (xi , β) −

Jij ∆βj
j=1

(C.5)

m

= ∆yi −

Jij ∆βj .
j=1

To find the minimum of the sum of squares of ri , the gradient of Equation (C.2) is set to
zero as
∂S
=2
∂βj

n

ri
i=1
n

= −2

∂ri
∂βj

∂f (xi , β)
ri
∂βj
i=1

n

= −2

i=1

= 0.

Jij ∆yi −

m
j=1



(C.6)

Jij ∆βj 

Rearranging Equation (C.6), we obtain m linear equations
n

m

n

Jij Jik ∆βk =
i=1 k=1

Jij ∆yi ,
i=1

173

(C.7)

which can be written in the matrix form as

JJT ∆β = JT ∆y.

174

(C.8)

AppendixD

Newton-Raphson Method for Multi-dimensio
Nonlinear Systems
A multi-dimensional nonlinear system with N variables is given as

Fi (x1 , x2 , · · · , xN ) = 0

i = 1, 2, · · · , N.

(D.1)

Using a vector x to denote values of xi and F to denote functions Fi , each function can be
expanded using Taylor’s series as [144]
N

Fi (x + δx) = Fi (x) +
J=1

∂Fi
δxj + O(δx2 ).
∂xj

(D.2)

The matrix of partial derivatives is the Jacobian matrix as

Jij =

∂Fi
.
∂xj

(D.3)

The Equation (D.2) can be rewritten in matrix form as

F(x + δx) = F(x) + J · δx + O(δx2 ).

175

(D.4)

By keeping only the first order term and setting F(x + δx) = 0, a set of linear equations to
calculate the increment of x can be obtained

J · δx = −F.

(D.5)

xnew = xold + δx.

(D.6)

The solution of x is updated by

176

BIBLIOGRAPHY

177

BIBLIOGRAPHY

[1] http://en.wikipedia.org/wiki/Least_squares.
[2] A. Agwai, I. Guven, and E. Madenci. Damage prediction for electronic package drop
test using finite element method and peridynamic theory. In Electronic Components
and Technology Conference, 2009. ECTC 2009. 59th, pages 565–569. IEEE, 2009.
[3] A. Agwai, I. Guven, and E. Madenci. Predicting crack initiation and propagation using
XFEM, CZM and peridynamics: A comparative study. In Electronic Components and
Technology Conference (ECTC), 2010 Proceedings 60th, pages 1178–1185. IEEE, 2010.
[4] E. Anderheggen, D. Ekchian, K. Heiduschke, and P. Bartelt. A contact algorithm
for explicit dynamic fem analysis. In Computational Techniques, Proc. 1st Int. Conf,
pages 271–283. Southampton, UK, 1993.
[5] T. L. Anderson. Fracture mechanics: Fundamentals and applications. CRC, 2005.
[6] C. E. Anderson Jr, D. L. Littlefield, and J. D. Walker. Long-rod penetration, target
resistance, and hypervelocity impact. International Journal of Impact Engineering,
14(1-4):1–12, 1993.
[7] C. E. Anderson Jr, D. L. Orphal, T. Behner, and D. W. Templeton. Failure and penetration response of borosilicate glass during short-rod impact. International Journal
of Impact Engineering, 36(6):789–798, 2009.
[8] E. Askari, J. Xu, and S.A. Silling. Peridynamic analysis of damage and failure in
composites. In 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada.
Reston, VA: AIAA, 2006.
[9] M. E. Backman and W. Goldsmith. The mechanics of penetration of projectiles into
targets. International Journal of Engineering Science, 16(1):1–99, 1978.
[10] G. I. Barenblatt. The mathematical theory of equilibrium cracks in brittle fracture.
Advances in applied mechanics, 7(55-129):104, 1962.
178

[11] K. J. Bathe. Finite element procedures. Prentice Hall, Englewood Cliffs, NJ, 2006.
[12] T. Behner, C. E. Anderson, D. L. Orphal, V. Hohler, M. Moll, and D. W. Templeton.
Penetration and failure of lead and borosilicate glass against rod impact. International
Journal of Impact Engineering, 35(6):447–456, 2008.
[13] T. Belytschko and T. Black. Elastic crack growth in finite elements with minimal
remeshing. International Journal for Numerical Methods in Engineering, 45(5):601–
620, 1999.
[14] T. Belytschko and J. I. Lin. A three-dimensional impact-penetration algorithm with
erosion. International Journal of Impact Engineering, 5(1-4):111–127, 1987.
[15] T. Belytschko and M. O. Neal. Contact-impact by the pinball algorithm with penalty
and lagrangian methods. International Journal for Numerical Methods in Engineering,
31(3):547–572, 1991.
[16] G. Ben-Dor, A. Dubinsky, and T. Elperin. Ballistic impact: recent advances in analytical modeling of plate penetration dynamics–a review. Applied Mechanics Reviews,
58:355, 2005.
[17] R. F. Bishop, R. Hill, and N. F. Mott. The theory of indentation and hardness tests.
Proceedings of the Physical Society, 57:147, 1945.
[18] S. J. Bless and C. E. Anderson Jr. Penetration of hard layers by hypervelocity rod
projectiles. International Journal of Impact Engineering, 14(1-4):85–93, 1993.
[19] S. J. Bless, Z. Rosenberg, and B. Yoon. Hypervelocity penetration of ceramics. International Journal of Impact Engineering, 5(1-4):165–171, 1987.
[20] F. Bobaru, M. Yang, L. F. Alves, S. A. Silling, E. Askari, and J. Xu. Convergence,
adaptive refinement, and scaling in 1D peridynamics. International Journal for Numerical Methods in Engineering, 77(6):852–877, 2009.
[21] T. Børvik, O. S. Hopperstad, T. Berstad, and M. Langseth. A computational model
of viscoplasticity and ductile damage for impact and penetration. European Journal
of Mechanics-A/Solids, 20(5):685–712, 2001.
[22] T. Børvik, O. S. Hopperstad, S. Dey, E. V. Pizzinato, M. Langseth, and C. Albertini.
Strength and ductility of weldox 460 E steel at high strain rates, elevated temperatures
and various stress triaxialities. Engineering fracture mechanics, 72(7):1071–1087, 2005.
179

[23] T. Børvik, O. S. Hopperstad, M. Langseth, and K. A. Malo. Effect of target thickness
in blunt projectile penetration of weldox 460 E steel plates. International journal of
impact engineering, 28(4):413–464, 2003.
[24] T. Børvik, M. Langseth, O. S. Hopperstad, and K. A. Malo. Ballistic penetration of
steel plates. International journal of impact engineering, 22(9-10):855–886, 1999.
[25] M. J. Buehler and H. Gao. Dynamical fracture instabilities due to local hyperelasticity
at crack tips. Nature, 439(7074):307–310, 2006.
[26] V. A. Buryachenko. Micromechanics of heterogeneous materials. Springer Verlag, 2007.
[27] C. A. Calder and W. Goldsmith. Plastic deformation and perforation of thin plates
resulting from projectile impact. International Journal of Solids and Structures,
7(7):863–868, 1971.
[28] N. J. Carpenter, R. L. Taylor, and M. G. Katona. Lagrange constraints for transient finite element surface contact. International Journal for Numerical Methods in
Engineering, 32(1):103–128, 1991.
[29] A. B. Chaudhary and K. J. Bathe. A solution method for static and dynamic analysis of
three-dimensional contact problems with friction. Computers & Structures, 24(6):855–
873, 1986.
[30] J. S. Chen, C. Pan, C. T. Wu, and W. K. Liu. Reproducing kernel particle methods
for large deformation analysis of non-linear structures. Computer Methods in Applied
Mechanics and Engineering, 139(1-4):195–227, 1996.
[31] C. Chinnaswamy, B. Amadei, and T. H. Illangasekare. A new method for finite element transitional mesh generation. International journal for numerical methods in
engineering, 31(7):1253–1270, 1991.
[32] J. D. Cinnamon, S. E. Jones, J. W. House, and L. L. Wilson. A one-dimensional
analysis of rod penetration. International journal of impact engineering, 12(2):145–
166, 1992.
[33] R. Citarella and F. G. Buchholz. Comparison of crack growth simulation by DBEM and
FEM for SEN-specimens undergoing torsion or bending loading. Engineering Fracture
Mechanics, 75(3-4):489–509, 2008.
[34] R. D. Cook, D. S. Malkus, M. E. Plesha, and R. J. Witt. Concepts and applications
of finite element analysis. A1bazaar, 2007.
180

[35] J. V. Cox. An extended finite element method with analytical enrichment for cohesive crack modeling. International Journal for Numerical Methods in Engineering,
78(1):48–83, 2009.
[36] K. Dayal and K. Bhattacharya. Kinetics of phase transformations in the peridynamic
formulation of continuum mechanics. Journal of the Mechanics and Physics of Solids,
54(9):1811–1842, 2006.
[37] S. De, J. W. Hong, and K. J. Bathe. On the method of finite spheres in applications:
towards the use with ADINA and in a surgical simulator. Computational Mechanics,
31(1):27–37, 2003.
[38] R. De Borst. Some recent developments in computational modelling of concrete fracture. International journal of fracture, 86(1):5–36, 1997.
[39] R. De Borst. Numerical aspects of cohesive-zone models. Engineering fracture mechanics, 70(14):1743–1757, 2003.
[40] P. N. Demmieand and S. A Silling. An Approach to Modeling Extreme Loading of
Structures using Peridynamics. Bulletin of the American Physical Society, 2006.
[41] M. Di Prisco, L. Ferrara, F. Meftah, J. Pamin, R. De Borst, J. Mazars, and J. M.
Reynouard. Mixed mode fracture in plain and reinforced concrete: some results on
benchmark tests. International journal of fracture, 103(2):127–148, 2000.
[42] J. F. Doyle. Determining the contact force during the transverse impact of plates.
Experimental Mechanics, 27(1):68–72, 1987.
[43] W. Duckworth. Discussion of paper by e. ryshkewitch,’compressive strength of porous
sintered alumina and zirconia’. Journal of the American Ceramic Society, 36:68, 1953.
[44] D. S. Dugdale. Yielding of steel sheets containing slits. Journal of the Mechanics and
Physics of Solids, 8(2):100–104, 1960.
[45] S. K. Dutta, A. K. Mukhopadhyay, and D. Chakraborty. Assessment of strength by
young’s modulus and porosity: a critical evaluation. Journal of the American Ceramic
Society, 71(11):942–947, 1988.
[46] D. G. B. Edelen, A. E. Green, and N. Laws. Nonlocal continuum mechanics. Archive
for Rational Mechanics and Analysis, 43(1):36–44, 1971.
181

[47] N. Elabbasi, J. W. Hong, and K. J. Bathe. On the reliable solution of contact problems
in engineering design. International Journal of Mechanics and Materials in design,
1(1):3–16, 2004.
[48] M. Elices, G. V. Guinea, J. Gomez, and J. Planas. The cohesive zone model: advantages, limitations and challenges. Engineering fracture mechanics, 69(2):137–163,
2002.
[49] E. Emmrich and O. Weckner. The peridynamic equation of motion in non-local elasticity theory. In III European Conference on Computational Mechanics. Solids, Structures
and Coupled Problems in Engineering, Lisbon,, Springer, volume 19, 2006.
[50] E. Emmrich and O. Weckner. On the well-posedness of the linear peridynamic model
and its convergence towards the navier equation of linear elasticity. Communications
in Mathematical Sciences, 5(4):851–864, 2007.
[51] F. Ercolessi. A molecular dynamics primer. Spring College in Computational Physics,
ICTP, Trieste, pages 24–25, 1997.
[52] M. J. Forrestal, B. S. Altman, J. D. Cargile, and S. J. Hanchak. An empirical equation
for penetration depth of ogive-nose projectiles into concrete targets. International
Journal of Impact Engineering, 15(4):395–405, 1994.
[53] M. J. Forrestal, D. J. Frew, J. P. Hickerson, and T. A. Rohwer. Penetration of concrete targets with deceleration-time measurements. International Journal of Impact
Engineering, 28(5):479–497, 2003.
[54] M. J. Forrestal, V. K. Luk, Z. Rosenberg, and N. S. Brar. Penetration of 7075-T651
aluminum targets with ogival-nose rods. International journal of solids and structures,
29(14-15):1729–1736, 1992.
[55] M. J. Forrestal, K. Okajima, and V. K. Luk. Penetration of 6061-T651 aluminum
targets with rigid long rods. Journal of applied mechanics, 55:755, 1988.
[56] M. J. Forrestal and A. J. Piekutowski. Penetration experiments with 6061-T6511
aluminum targets and spherical-nose steel projectiles at striking velocities between 0.5
and 3.0 km/s. International journal of impact engineering, 24(1):57–67, 2000.
[57] M. J. Forrestal, D. Y. Tzou, E. Askari, and D. B. Longcope. Penetration into ductile
metal targets with rigid spherical-nose rods. International journal of impact engineering, 16(5-6):699–710, 1995.
182

[58] J. T. Foster, S. A. Silling, and W. W. Chen. Viscoplasticity using peridynamics.
International Journal for Numerical Methods in Engineering, 81(10):1242–1258, 2009.
[59] J. T. Foster, S. A. Silling, and W. W. Chen. An energy based failure criterion for
use with peridynamic states. International Journal for Multiscale Computational Engineering, 9:112, 2011.
[60] J. W. Foulk, D. H. Allen, and K. L. E. Helms. Formulation of a three-dimensional
cohesive zone model for application to a finite element algorithm. Computer Methods
in Applied Mechanics and Engineering, 183(1-2):51–66, 2000.
[61] D. J. Frew, M. J. Forrestal, and S. J. Hanchak. Penetration experiments with limestone
targets and ogive-nose steel projectiles. Journal of applied mechanics, 67:841–845, 2000.
[62] D. J. Frew, S. J. Hanchak, M. L. Green, and M. J. Forrestal. Penetration of concrete targets with ogive-nose steel rods. International Journal of Impact Engineering,
21(6):489–497, 1998.
[63] H. Gao, Y. Huang, and F. F. Abraham. Continuum and atomistic studies of intersonic
crack propagation. Journal of the Mechanics and Physics of Solids, 49(9):2113–2132,
2001.
[64] T. C. Gasser and G. A. Holzapfel. 3D crack propagation in unreinforced concrete:
A two-step algorithm for tracking 3D crack paths. Computer Methods in Applied
Mechanics and Engineering, 195(37-40):5198–5219, 2006.
[65] W. Gerstle, N. Sau, and S. A. Silling. Peridynamic modeling of concrete structures.
Nuclear Engineering and Design, 237(12-13):1250–1258, 2007.
[66] W. Goldsmith. Impact: the theory and physical behaviour of colliding solids. Dover
Pubns, 1960.
[67] J. T. Gomez and A. Shukla. Multiple impact penetration of semi-infinite concrete.
International journal of impact engineering, 25(10):965–979, 2001.
[68] F. I. Grace. Long-rod penetration into targets of finite thickness at normal impact.
International journal of impact engineering, 16(3):419–433, 1995.
[69] Y. D. Ha and F. Bobaru. Studies of dynamic crack propagation and crack branching
with peridynamics. International Journal of Fracture, 162(1-2):229–244, 2010.
183

[70] Y. D. Ha and F. Bobaru. Characteristics of dynamic brittle fracture captured with
peridynamics. Engineering Fracture Mechanics, 78(6):1156 – 1168, 2011.
[71] C. Hahn and P. Wriggers. An explicit multi-body contact algorithm. PAMM, 3(1):280–
281, 2003.
[72] J. O. Hallquist. LS-DYNA theory manual. Livermore Software Technology Corporation, 2006.
[73] M. J. Harris. Fast fluid dynamics simulation on the GPU. GPU gems, 1:637–665, 2004.
[74] D. P. H. Hasselman. Relation between effects of porosity on strength and on young’s
modulus of elasticity of polycrystalline materials. Journal of the American Ceramic
Society, 46(11):564–565, 1963.
[75] M. W. Heinstein, F. J. Mello, S. W. Attaway, and T. A. Laursen. Contact–impact
modeling in explicit transient dynamics. Computer methods in applied mechanics and
Engineering, 187(3-4):621–640, 2000.
[76] S. A. Hill. Determination of an empirical model for the prediction of penetration hole
diameter in thin plates from hypervelocity impact. International journal of impact
engineering, 30(3):303–321, 2004.
[77] V. Hohler, A. J. Stilp, and K. Weber. Hypervelocity penetration of tungsten sinteralloy rods into aluminum. International journal of impact engineering, 17(1-3):409–418,
1995.
[78] B. L. Holian and R. Ravelo. Fracture simulations using large-scale molecular dynamics.
Physical Review B, 51(17):11275–11288, May 1995.
[79] J. W. Hong and K. J. Bathe. On analytical transformations for efficiency improvements
in the method of finite spheres. Computational fluid and solid mechanics, pages 1990–
1994, 2003.
[80] J. W. Hong and K. J. Bathe. Coupling and enrichment schemes for finite element and
finite sphere discretizations. Computers & Structures, 83(17-18):1386–1395, 2005.
[81] W. Hu, Y. D. Ha, and F. Bobaru. Modeling dynamic fracture and damage in a fiberreinforced composite lamina with peridynamics. International Journal for Multiscale
Computational Engineering, 9(6), 2011.
184

[82] W. Hu, Y. D. Ha, and F. Bobaru. Peridynamic model for dynamic fracture in unidirectional fiber-reinforced composites. Computer Methods in Applied Mechanics and
Engineering, 217:247–261, 2012.
[83] R. Huang, N. Sukumar, and J. H. Pr´vost. Modeling quasi-static crack growth with the
e
extended finite element method Part II: Numerical applications. International journal
of solids and structures, 40(26):7539–7552, 2003.
[84] T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, and W. Kanoknukulchai.
A finite element method for a class of contact-impact problems. Computer Methods in
Applied Mechanics and Engineering, 8(3):249–276, 1976.
[85] A. Ibrahimbegovic and A. Delaplace. Microscale and mesoscale discrete models for
dynamic fracture of structures built of brittle material. Computers & structures,
81(12):1255–1265, 2003.
[86] G. R. Irwin. Linear fracture mechanics, fracture transition, and fracture control. Engineering Fracture Mechanics, 1(2):241–257, 1968.
[87] G. R. Irwin and D. E. Wit. A summary of fracture mechanics concepts. Journal of
Testing and Evaluation, 11:56–65, 1983.
[88] S. E. Jones, P. P. Gillis, and J. C. Foster. On the penetration of semi-infinite targets
by long rods. Journal of the mechanics and physics of solids, 35(1):121–131, 1987.
[89] E. P. Kearsley and P. J. Wainwright. The effect of porosity on the strength of foamed
concrete. Cement and concrete research, 32(2):233–239, 2002.
[90] B. Kilic. Peridynamic theory for progressive failure prediction in homogeneous and
heterogeneous materials. PhD thesis, The University of Arizona, 2009.
[91] B. Kilic, A. Agwai, and E. Madenci. Peridynamic theory for progressive damage
prediction in center-cracked composite laminates. Composite Structures, 90(2):141–
151, 2009.
[92] B. Kilic and E. Madenci. Prediction of crack paths in a quenched glass plate by using
peridynamic theory. International Journal of Fracture, 156(2):165–177, 2009.
[93] B. Kilic and E. Madenci. Structural stability and failure analysis using peridynamic
theory. International Journal of Non-Linear Mechanics, 44(8):845–854, 2009.
185

[94] B. Kilic and E. Madenci. Coupling of peridynamic theory and the finite element
method. Journal of mechanics of materials and structures, 5,No.5:707–733, 2010.
[95] F. P. Knudsen. Dependence of mechanical strength of brittle polycrystalline specimens
on porosity and grain size. Journal of the American Ceramic Society, 42(8):376–387,
1959.
[96] V. Kouznetsova, W. A. M. Brekelmans, and F. P. T. Baaijens. An approach to micromacro modeling of heterogeneous materials. Computational Mechanics, 27(1):37–48,
2001.
[97] V. Kouznetsova, M. G. D. Geers, and W. A. M. Brekelmans. Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational
homogenization scheme. International Journal for Numerical Methods in Engineering,
54(8):1235–1260, 2002.
[98] J. Kozicki and J. Tejchman. Modelling of fracture process in concrete using a novel
lattice model. Granular Matter, 10(5):377–388, 2008.
[99] V. D. Krstic. Porosity dependence of strength in brittle solids. Theoretical and applied
fracture mechanics, 10(3):241–247, 1988.
[100] V. D. Krstic. Effect of microstructure on fracture of brittle materials: Unified approach.
Theoretical and applied fracture mechanics, 45(3):212–226, 2006.
[101] J. Kruger and R. Westermann. Linear algebra operators for GPU implementation of
numerical algorithms. In ACM SIGGRAPH 2005 Courses, page 234. ACM, 2005.
[102] R. Kumar and B. Bhattacharjee. Porosity, pore size distribution and in situ strength
of concrete. Cement and concrete research, 33(1):155–164, 2003.
[103] P. Lall, S. Shantaram, and D. Panchagade. Peridynamic-models using finite elements
for shock and vibration reliability of leadfree electronics. In Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm), 2010 12th IEEE Intersociety
Conference on, pages 1–12. IEEE, 2010.
[104] R. B. Lehoucq and S. A. Silling. Force flux and the peridynamic stress tensor. Journal
of the Mechanics and Physics of Solids, 56(4):1566–1577, 2008.
[105] R. W. Lewis, B. A. Shrefler, and R. W. Lewis. The finite element method in the static
and dynamic deformation and consolidation of porous media. John Wiley, 1998.
186

[106] S. Li and W. K. Liu. Meshfree and particle methods and their applications. Applied
Mechanics Reviews, 55:1, 2002.
[107] D. J. Littlewood. Simulation of dynamic fracture using peridynamics, finite element
modeling, and contact. In Proceedings of the ASME 2010 International Mechanical
Engineering Congress and Exposition (IMECE), 2010.
[108] G. R. Liu and Y. T. Gu. Coupling of element free galerkin and hybrid boundary
element methods using modified variational formulation. Computational Mechanics,
26(2):166–173, 2000.
[109] G. R. Liu and M. B. Liu. Smoothed particle hydrodynamics: a meshfree particle method.
World Scientific Pub Co Inc, 2003.
[110] W. Liu and J. W. Hong. Three-dimensional lamb wave propagation excited by a phased
piezoelectric array. Smart Materials and Structures, 19:085002, 2010.
[111] W. K. Liu, Y. Chen, S. Jun, J. S. Chen, T. Belytschko, C. Pan, R. A. Uras, and
C. T. Chang. Overview and applications of the reproducing kernel particle methods.
Archives of Computational Methods in Engineering, 3(1):3–80, 1996.
[112] Livermore Software Technology Corporation. LS-DYNA keyword user’s manual, 2007.
[113] G. Lubineau, Y. Azdoud, F. Han, C. Rey, and A. Askari. A morphing strategy to
couple non-local to local continuum mechanics. Journal of the Mechanics and Physics
of Solids, 2012.
[114] P. Lundberg, R. Renstr¨m, and B. Lundberg. Impact of metallic projectiles on ceramic
o
targets: transition between interface defeat and penetration. International Journal of
Impact Engineering, 24(3):259–275, 2000.
[115] R. W. Macek and S. A. Silling. Peridynamics via finite element analysis. Finite
Elements in Analysis and Design, 43(15):1169–1178, 2007.
[116] S. Manavski and G. Valle. CUDA compatible GPU cards as efficient hardware accelerators for Smith-Waterman sequence alignment. BMC bioinformatics, 9(Suppl 2):S10,
2008.
[117] S. Mariani and U. Perego. Extended finite element method for quasi-brittle fracture.
International Journal for Numerical Methods in Engineering, 58(1):103–126, 2003.
187

[118] G. T. Mase, R. E. Smelser, and G. E. Mase. Continuum mechanics for engineers.
CRC, 1999.
[119] P. K. Mehta and P. J. M. Monteiro. Concrete: structure, properties, and materials.
Prentice-Hall, 1986.
[120] J. M. Melenk and I. Babuska. The partition of unity finite element method: basic
theory and applications. Computer methods in applied mechanics and engineering,
139(1-4):289–314, 1996.
[121] G. Meschke and P. Dumstorff. Energy-based modeling of cohesive and cohesionless
cracks via X-FEM. Computer methods in applied mechanics and engineering, 196(2124):2338–2357, 2007.
[122] N. Mo¨s and T. Belytschko. Extended finite element method for cohesive crack growth.
e
Engineering Fracture Mechanics, 69(7):813–833, 2002.
[123] N. Mo¨s, J. Dolbow, and T. Belytschko. A finite element method for crack growth
e
without remeshing. International Journal for Numerical Methods in Engineering,
46(1):131–150, 1999.
[124] J. J. Monaghan. Smoothed particle hydrodynamics. Reports on Progress in Physics,
68:1703, 2005.
[125] H. Moslemi and A. R. Khoei. 3D adaptive finite element modeling of non-planar curved
crack growth using the weighted superconvergent patch recovery method. Engineering
Fracture Mechanics, 76(11):1703–1728, 2009.
[126] L. E. Murr, S. A Quinones, E. Ferreyra T, A. Ayala, O. L Valerio, F. H¨rz, and R. P.
o
Bernhard. The low-velocity-to-hypervelocity penetration transition for impact craters
in metal targets. Materials Science and Engineering: A, 256(1):166–182, 1998.
[127] V. Murti and S. Valliappan. Numerical inverse isoparametric mapping in remeshing
and nodal quantity contouring. Computers & structures, 22(6):1011–1021, 1986.
[128] V. Murti, Y. Wang, and S. Valliappan. Numerical inverse isoparametric mapping in
3d fem. Computers & structures, 29(4):611–622, 1988.
[129] L. F. Nielsen. Strength and stiffness of porous materials. Journal of the American
Ceramic Society, 73(9):2684–2689, 1990.
188

[130] M. B. Nooru-Mohamed, E. Schlangen, and J. G. M. van Mier. Experimental and
numerical study on the behavior of concrete subjected to biaxial tension and shear.
Advanced Cement Based Materials, 1(1):22–37, 1993.
[131] NVIDIA. NVIDIA CUDA Programming Guide, version 2.3.1 edition, 2009.
[132] E. T. Ooi and Z. J. Yang. Modelling dynamic crack propagation using the scaled
boundary finite element method. International Journal for Numerical Methods in
Engineering, 2011.
[133] E. Oterkus. Peridynamic theory for modeling three-dimensional damage growth in
metallic and composite structures. 2010.
[134] M. L. Parks, R. B. Lehoucq, S. J. Plimpton, and S. A. Silling. Implementing peridynamics within a molecular dynamics code. Computer Physics Communications,
179(11):777–783, 2008.
[135] M. L. Parks, S. J. Plimpton, R. B. Lehoucq, and S. A. Silling. Peridynamics with
LAMMPS: A user guide. Technical report, Technical Report SAND 2008-1035, Sandia
National Laboratories, 2008.
[136] B. Patz´k and M. Jir´sek. Adaptive resolution of localized damage in quasi-brittle
a
a
materials. Journal of engineering mechanics, 130:720, 2004.
[137] K. K. Phani and S. K. Niyogi. Young’s modulus of porous brittle solids. Journal of
materials science, 22(1):257–263, 1987.
[138] E. H. Phillips, R. L. Davis, and J. D. Owens. Unsteady Turbulent Simulations on
a Cluster of Graphics Processors. In Proceedings of the 40th AIAA Fluid Dynamics
Conference, number AIAA-2010-5036. 40th AIAA Fluid Dynamics Conference, 2010.
[139] A. J. Piekutowski, M. J. Forrestal, K. L. Poormon, and T. L. Warren. Penetration
of 6061-T6511 aluminum targets by ogive-nose steel projectiles with striking velocities
between 0.5 and 3.0 km/s. International journal of impact engineering, 23(1):723–734,
1999.
[140] P. Pivonka, J. Oˇbolt, R. Lackner, and H. A. Mang. Comparative studies of 3Dz
constitutive models for concrete: application to mixed-mode fracture. International
journal for numerical methods in engineering, 60(2):549–570, 2004.
189

[141] J. Planas, M. Elices, G. V. Guinea, F. J. G´mez, D. A. Cend´n, and I. Arbilla. Generalo
o
izations and specializations of cohesive crack models. Engineering fracture mechanics,
70(14):1759–1776, 2003.
[142] The Portland Group. PGI Fortran & C Accelerator Programming Model v1.2, 1.2
edition, 3 2010.
[143] The Portland Group. PGI User’s Guide, release 2010 edition, 2010.
[144] W. H. Press. Numerical recipes in FORTRAN: the art of scientific computing, volume 1. Cambridge Univ Pr, 1992.
[145] R. F. Recht and T. W. Ipson. Ballistic perforation dynamics. Journal of Applied
Mechanics, 30:384, 1963.
[146] J. R´thor´, S. Roux, and F. Hild. Mixed-mode crack propagation using a hybrid analyte
e
ical and extended finite element method. Comptes Rendus M´canique, 338(3):121–126,
e
2010.
[147] Z. Rosenberg and E. Dekel. On the deep penetration and plate perforation by rigid
projectiles. International Journal of Solids and Structures, 46(24):4169–4180, 2009.
[148] Z. Rosenberg, E. Dekel, V. Hohler, A. J. Stilp, and K. Weber. Hypervelocity penetration of tungsten alloy rods into ceramic tiles: experiments and 2-D simulations.
International journal of impact engineering, 20(6-10):675–683, 1997.
[149] M. R¨ssler and I. Odler. Investigations on the relationship between porosity, structure
o
and strength of hydrated portland cement pastes: I. effects of porosity. Cement and
Concrete Research, 15:320–330, 1985.
[150] C. L. Rountree, R. K. Kalia, E. Lidorikis, A. Nakano, L. Van Brutzel, and P. Vashishta.
Atomistic aspects of crack propagation in brittle materials: Multimillion atom molecular dynamics simulations. Annual Review of Materials Research, 32(1):377–400, 2002.
[151] G. Ruiz, A. Pandolfi, and M. Ortiz. Three-dimensional cohesive modeling of dynamic
mixed-mode fracture. International Journal for Numerical Methods in Engineering,
52(1-2):97–120, 2001.
[152] E. Ryshkewitch. Compression strength of porous sintered alumina and zirconia. Journal of the American Ceramic Society, 36(2):65–68, 1953.
190

[153] S. A. Silling. Reformulation of elasticity theory for discontinuities and long-range
forces. Journal of the Mechanics and Physics of Solids, 48(1):175–209, 2000.
[154] S. A. Silling. Linearized theory of peridynamic states. Journal of Elasticity, 99(1):85–
111, 2010.
[155] S. A. Silling and E. Askari. A meshfree method based on the peridynamic model of
solid mechanics. Computers & Structures, 83(17-18):1526–1535, 2005.
[156] S. A. Silling, M. Epton, O. Weckner, J. Xu, and E. Askari. Peridynamic states and
constitutive modeling. Journal of Elasticity, 88(2):151–184, 2007.
[157] S. A. Silling and R. B. Lehoucq. Convergence of peridynamics to classical elasticity
theory. Journal of Elasticity, 93(1):13–37, 2008.
[158] S. A. Silling and R. B. Lehoucq. Peridynamic theory of solid mechanics. Advances in
Applied Mechanics, 44:73 – 168, 2010.
[159] S. A. Silling and R. B. Lehoucq. Peridynamic theory of solid mechanics. Advances in
Applied Mechanics, 44:73–168, 2010.
[160] S. A. Silling, O. Weckner, E. Askari, and F. Bobaru. Crack nucleation in a peridynamic
solid. International Journal of Fracture, 162(1):219–227, 2010.
[161] S. A. Silling, M. Zimmermann, and R. Abeyaratne. Deformation of a peridynamic bar.
Journal of Elasticity, 73(1):173–190, 2003.
[162] R. J. M. Smit, W. A. M. Brekelmans, and H. E. H. Meijer. Prediction of the large-strain
mechanical response of heterogeneous polymer systems: local and global deformation
behaviour of a representative volume element of voided polycarbonate. Journal of the
Mechanics and Physics of Solids, 47(2):201–221, 1999.
[163] R. M. Spriggs. Expression for effect of porosity on elastic modulus of polycrystalline
refractory materials, particularly aluminum oxide. Journal of the American Ceramic
Society, 44(12):628–629, 1961.
[164] T. A. C. Stock and K. R. L. Thompson. Penetration of aluminum alloys by projectiles.
Metallurgical and Materials Transactions B, 1(1):219–224, 1970.
[165] W. J. Stronge. Impact mechanics. Cambridge University Press, 2004.
191

[166] R. Subramanian, S. J. Bless, J. Cazamias, and D. Berry. Reverse impact experiments
against tungsten rods and results for aluminum penetration between 1.5 and 4.2 km/s.
International journal of impact engineering, 17(4-6):817–824, 1995.
[167] N. Sukumar, N. Mo¨s, B. Moran, and T. Belytschko. Extended finite element method
e
for three-dimensional crack modelling. International Journal for Numerical Methods
in Engineering, 48(11):1549–1570, 2000.
[168] N. Sukumar and J. H. Pr´vost. Modeling quasi-static crack growth with the extended
e
finite element method Part I: Computer implementation. International Journal of
Solids and Structures, 40(26):7513–7537, 2003.
[169] Y. Sumi and Z. N. Wang. A finite-element simulation method for a system of growing
cracks in a heterogeneous material. Mechanics of Materials, 28(1-4):197 – 206, 1998.
[170] A. Tate. A theory for the deceleration of long rods after impact. Journal of the
Mechanics and Physics of Solids, 15(6):387–399, 1967.
[171] A. Tate. Further results in the theory of long rod penetration. Journal of the Mechanics
and Physics of Solids, 17(3):141–150, 1969.
[172] E. A. Taylor, K. Tsembelis, C. J. Hayhurst, L. Kay, and M. J. Burchell. Hydrocode
modelling of hypervelocity impact on brittle materials: depth of penetration and conchoidal diameter. International journal of impact engineering, 23(1):895–904, 1999.
[173] T. G. Trucano and D. E. Grady. Impact shock and penetration fragmentation in porous
media. International journal of impact engineering, 17(4-6):861–872, 1995.
[174] J. F. Unger, S. Eckardt, and C. Konke. Modelling of cohesive crack growth in concrete
structures with the extended finite element method. Computer methods in applied
mechanics and engineering, 196(41-44):4087–4100, 2007.
¨
[175] O. Vardar, I. Finnie, D.R. Biswas, and RM Fulrath. Effect of spherical pores on the
strength of a polycrystalline ceramic. International Journal of Fracture, 13(2):215–223,
1977.
[176] A. S. Wagh, J. P. Singh, and R. B. Poeppel. Dependence of ceramic fracture properties
on porosity. Journal of materials science, 28(13):3589–3593, 1993.
[177] J. C. Wang. Young’s modulus of porous materials: Theoretical derivation of modulusporosity correlation. Journal of materials science, 19(3):801–808, 1984.
192

[178] J. C. Wang. Young’s modulus of porous materials: Young’s modulus of porous alumina with changing pore structure. Journal of Materials Science, 19:809–814, 1984.
10.1007/BF00540452.
[179] T. L. Warren and M. J. Forrestal. Effects of strain hardening and strain-rate sensitivity on the penetration of aluminum targets with spherical-nosed rods. International
journal of solids and structures, 35(28):3737–3753, 1998.
[180] T. L. Warren and K. L. Poormon. Penetration of 6061-T6511 aluminum targets by
ogive-nosed VAR 4340 steel projectiles at oblique angles: experiments and simulations.
International journal of impact engineering, 25(10):993–1022, 2001.
[181] T. L. Warren, S. A. Silling, A. Askari, O. Weckner, M. A. Epton, and J. Xu. A nonordinary state-based peridynamic method to model solid material deformation and
fracture. International Journal of Solids and Structures, 46(5):1186–1195, 2009.
[182] T. L. Warren and M. R. Tabbara. Simulations of the penetration of 6061-T6511
aluminum targets by spherical-nosed VAR 4340 steel projectiles. International journal
of solids and structures, 37(32):4419–4435, 2000.
[183] O. Weckner and R. Abeyaratne. The effect of long-range forces on the dynamics of a
bar. Journal of the Mechanics and Physics of Solids, 53(3):705–728, 2005.
[184] J. P. Weiler, J. T. Wood, R. J. Klassen, E. Maire, R. Berkmortel, and G. Wang.
Relationship between internal porosity and fracture strength of die-cast magnesium
am60b alloy. Materials Science and Engineering: A, 395(1):315–322, 2005.
[185] M. L. Wilkins. Mechanics of penetration and perforation. International Journal of
Engineering Science, 16(11):793–807, 1978.
[186] M. Wolfe. The pgi accelerator programming model on NVIDIA GPUs. Technical
report, The Portland Group, 2009.
[187] P. Wriggers. Computational contact mechanics. Springer Verlag, 2006.
[188] J. Xu, A. Askari, O. Weckner, and S. A. Silling. Peridynamic analysis of impact damage
in composite laminates. Journal of Aerospace Engineering, 21:187, 2008.
[189] M. Zimmermann. A continuum theory with long-range forces for solids. PhD thesis,
Massachusetts Institute of Technology, 2005.

193