DEVELOPMENT AND APPLICATIONS OF NEW PERIDYNAMIC MODELS
By
Tao Jia

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

ABSTRACT
DEVELOPMENT AND APPLICATIONS OF NEW PERIDYNAMIC MODELS
 
By
Tao Jia

Studies of solid mechanics are traditionally based on continuum mechanics which utilizes
differential equations. Differential equations, however, become troublesome when damage takes
place. Instead of differential equations, this study presents a theory, so-called peridynamics,
based on integral equations. Similar to molecular dynamics, peridynamics assumes that the
domain of interest is organized by points. Each point interacts with every other point within a
horizon through a bond. Damage at a point takes place when a critical amount of bonds
associated with the point are broken. Besides the bond strength, peridynamics does not impose
any additional damage theory such as fracture mechanics used in continuum mechanics.
Peridynamics is still in its infant stage. New models need to be developed. In this study, a onedimensional model was firstly proposed and verified by an associated solution based on
continuum mechanics. The model was then used to simulate wave propagations in split
Hopkinson’s pressure bar (SHPB) and was validated by the experiment results. The peridynamic
model was then used for designing required shapers in SHPB application and greatly improves
the experiment efficiency. Secondly, a two-dimensional model was proposed and verified by an
associated solution based on continuum mechanics. Two computational algorithms were then
proposed and incorporated into peridynamic programming to significantly improve its
computational efficiency. The uses of the two-dimensional peridynamic model for simulating
dynamic damage progression were favorably validated by experiments. A four-parameter
peridynamic model was finally presented for investigating orthotropic materials. The model was

verified by analyses involving uniaxial tension and vibration. It was also validated with singleedge-notch testing results.

ACKNOWLEDGEMENTS

I would like to thank my professor, Dr. Dahsin Liu, for both the financial support of my research
and the directing of the research work. His insight allowed me to find suitable topics for research
and complete this dissertation in a timely manner. I would also like to thank Dr. Ronald Averill,
Dr. Thomas Pence and Dr. Zhengfang Zhou for agreeing to serve on my dissertation committee.
Lastly, I would especially like to thank my parents for their love and support through the years.

iv 

TABLE OF CONTENTS

LIST OF TABLES ………………………………………………………………………………vii
LIST OF FIGURES……………………………………………………………………………..viii
Chapter 1 Introduction………………………………………………………………………….1
1.1 Background……………………………………………………………………………1
1.2 Literature review………………………………………………………………………2
1.3 Motivations and objectives……………………………………………………………5
1.4 Outline of dissertation..………………………………………………………………..6
References…………………………………………………………………………………8
Chapter 2 Basic Formulations……...………………………………………………………….11
2.1 Peridynamic model of a continuum………………………………………………….11
2.2 Constitutive modeling………………………………………………………………..12
2.3 Failure criterion………………………………………………………………………13
2.4 Numerical method……………………………………………………………………14
Figures……………………………………………………………………………………15
References………………………………………………………………………………..19

Chapter 3 One-Dimensional Wave Propagation Analysis……..…………………………….21
3.1 Identification of stiffness……………………………………...………………….….21
3.2 Definition of stress in peridynamics………………………………………................22
3.3 Comparison with analytical solution………………………………………………...23
3.3.1 Analytical solution…………………………………………………………23
3.3.2 Numerical studies…………………………………………………………..26
3.4 Comparison with experiment results……..………………………………………….28
3.4.1 Square wave input………………………………………………………….29
3.4.2 Using filter for computational results…………………………………...…33
3.4.3 Using the incidence wave as the computational input…………………......35
3.5 Simulation with a striker……………………………………………………………..35
3.6 Validation with results in literature………………………………………….……….39
3.7 Modeling shaper in SHPB……………………………………………………….…..41
3.8 Conclusions………………………………………………………………………..…42
Figures…………………………………………………………………………………....44
References……………………………………………………………………………..…85
Chapter 4 Two-Dimensional Peridynamic Studies…………………………………………...87
4.1 Bar model…………………………………………………………………………….87
4.2 Definition of stress in peridynamics…………………………………………………89
4.3 Beam model………………………………………………………………………….91
4.4 Computational efficiency………………………..…………………………………...94
v 

4.4.1 Using matrix computation………………………………………………….94
4.4.2 Using parallel computing…………………………………………………..99
4.5 Comparison with vibration theory………………………………………………….101
4.5.1Vibration theory………………………………………………………..….101
4.5.2 Numerical studies………………………………………………………....102
4.6 Convergence to the one-dimensional results…………………………………...…..104
4.7 Axisymmetric problems…………………………………………………………….105
4.8 Failure theory……………………………………………………………………….107
4.9 Comparison of crack propagation velocity with experiments……………………...108
4.10 Dynamic crack branching…………………………………………………………109
4.11 Edge-cracked plate under impulsive loading……………………………………...110
4.12 Conclusions………………………………………………………………………..111
Figures…………………………………………………………………………………..113
References………………………………………………………………………………149

Chapter 5 Orthotropic Model Analysis..………………………………………………...…..152
5.1 Bar model for orthotropic materials………………………………………………...153
5.2 Beam model for orthotropic materials……………………………………………...160
5.3 Calculation of stress from peridynamics……………………………………………161
5.4 Laminated plate under static loading…………………………………………….…162
5.5 Free vibration of a laminated beam………………………………………………...164
5.6 Comparison of crack propagation velocity with experiments…………………..….167
5.7 Dynamic fracture mode in unidirectional composites…………………………..….169
5.8 Conclusions………………………………………………………………………....171
Figures…………………………………………………………………………………..173
References………………………………………………………………………………195
Chapter 6 Conclusions and Recommendations…………………………...…….…………...198
6.1 Conclusions…………………………………………………………………………198
6.2 Recommendations…………………………………………………………………..201
References………………………………………………………………………………204

vi 

LIST OF TABLES

Table 3.1 Wave period (

) based on various of

and

………………………………..……27

Table 3.2 Dimensions and material properties of the SHPB in Fig. 3.7…...................................28
Table 3.3 Dimensions and material properties of the split bar in Fig. 3.28……………………...36
Table 3.4 Dimensions and material properties of the SHPB in Fig. 3.33………………………..40
Table 3.5 Dimensions and material properties of the SHPB used ……………….……………...41
Table 4.1 Twelve pairs of

and …………………………………………………………….103

Table 5.1 Simplified Eqn. 5.9 in terms of independent terms………………………………….155
Table 5.2 Simplified Eqn. 5.12 in terms of independent terms………………………………...157
Table 5.3 Material propteties of E-Glass/Epoxy………………………………………………..163
Table 5.4 Material propteties of Kevlar/Epoxy………………………………………………...164
Table 5.5 Material propteties of graphite/epoxy……..…………………………..………..……168
Table 5.6 Material propteties of carbon/epoxy …………………………………...……….…...170

vii 

LIST OF FIGURES

Figure 2.1 Horizon of point …………………………….....................................……………..15
Figure 2.2 Reference configuration and current configuration………………………..…………16
Figure 2.3 Finite summation domain in numerical computation…………………...……………17
Figure 2.4 Flow chart of a peridynamic program…………………………………...…………..18
Figure 3.1 Calculation of stress…………………………………………………………...……..44
Figure 3.2 Wave propagation in a one-dimensional bar …………………………………….…..45
Figure 3.3 Analytical solution for wave propagation……………………………...…………….46
Figure 3.4 Number of nodes within a horizon …………………………………………………..47
Figure 3.5 Comparison between classical mechanics and peridynamics………………….…….48
Figure 3.6 Convergence of time step size………………………………………………………..49
Figure 3.7 Set up of Split Hopkinson pressure bar (SHPB)……………………………….….....50
Figure 3.8 Wave propagation in the SHPB………………………………………………………51
Figure 3.9 Square wave input for computational simulation ……………………………..……..52
Figure 3.10 Comparison in Study One with a square wave input ………………………………53
Figure 3.11 Study Two: Two Bar ……………………………………………………………….54
Figure 3.12 Study Two: comparison between experiment results and numerical solution, strain
history in incidence bar (top) and strain history in transmission bar (bottom)…………………..55
Figure 3.13 Study Two: comparison between experiment results and numerical solution with
grease, strain history in incidence bar (top) and strain history in transmission bar (bottom)…....56
Figure 3.14 A bond crossing an interface ……………………………………..………………...57
Figure 3.15 Study Three: aluminum specimen with the same cross-section area as the split
bars ………………………………………………………………………………………………58
Figure 3.16 Study Three: Aluminum specimen with the same-cross section area as the split bars,
incidence bar strain history (top) and transmission bar strain history(bottom)………….………59
viii 

Figure 3.17 Study Four: Steel specimen with a smaller cross-section area ……………………..60
Figure 3.18 Study Four: Steel specimen with a smaller cross-section area, incidence bar strain
history (top) and transmission bar strain history (bottom)…………………………...…………..61
Figure 3.19 Study One: computational results with filter compared with experiment results......62
Figure 3.20 Study Two: Computational results with filter compared with experiment results …63
Figure 3.21 Study Three: Computational results with filter compared with experiment results...64
Figure 3.22 Study Four: Computational results with filter compared with experiment results….65
Figure 3.23 Incidence wave from Study Two ……………………………………..…………….66
Figure 3.24 Study One: Incident wave as computational input …………………………………67
Figure 3.25 Study Two: Incident wave as computational input ………………………………...68
Figure 3.26 Study Three: Incident wave as computational input ……………………………….69
Figure 3.27 Study Four: Incident wave as computational input ………………………………...70
Figure 3.28 Split bar system for validating impact process …………………………………..…71
Figure 3.29 Simulation of impact process……………………………………………………….72
Figure 3.30 Comparison between experiment and simulation with velocity input……………...73
Figure 3.31 Impact of two cylindrical bars right before impact …………………..…………….74
Figure 3.32 Impact of two cylindrical bars after impact ………………………………………...75
Figure 3.33 SHPB system in [6]…………………………………………………………………76
Figure 3.34 Stress-strain relation of the aluminum specimen from [6]………………………….77
Figure 3.35 Comparing computational results with experiment results from [6]…………...…...78
Figure 3.36 Computational results with filter compared with experiment results……...………..79
Figure 3.37 Incidence wave after filter from Ref. [8]……………………………………..……..80
Figure 3.38 SHPB with a shaper…………………………………………………………………81

ix 

Figure 3.39 Split bar system in Ref. [8]………………………………... ……………………….82
Figure 3.40 Stress-strain relation for the copper shaper given in Ref. [8]…………………..…...83
Figure 3.41 Simulation results compared with experiment results from Ref. [8]………………..84
Figure 4.1 Two-dimensional domain under radial deformation………………………………..113
Figure 4.2 Definition of stress in peridynamics………………………………………………...114
Figure 4.3 Calculation of stress on point x……………………………………………………..115
Figure 4.4 Bond force in beam model ……………………………………...………………….116
Figure 4.5 Original configuration and current configuration…………………………………...117
Figure 4.6 Flow chart of algorithm in peridynamic matrix calculation ………………………..118
Figure 4.7 Flow chart for parallel computing…………………………………………………..119
Figure 4.8 Code before parallel modification and code after parallel modification……………120
Figure 4.9 Vibration of a circular plate…………………………………………………………121
Figure 4.10 Strain history on a circular plate...............................................................................122
Figure 4.11 Convergence of numerical solution to analytical solution ………………………..123
Figure 4.12 Convergence of numerical solution to analytical solution ………………………..124
Figure 4.13 Convergence of numerical solution to vibration theory ……………………..……125
Figure 4.14 Peridynamic calculation of strains at point A compared with results from vibration
theory …………………………………………………………………………………………..126
Figure 4.15 Peridynamic calculation of strain at point B compared with results from vibration
theory ……………………………………………………………………………………...…...127
Figure 4.16 Peridynamic calculation of strain at point C compared with results from vibration
theory …………………………………………………………………………………...……...128
Figure 4.17 Rectangular plate for free vibration study ………………………………………...129
Figure 4.18 One-dimensional result compared with two-dimensional result for aspect ratio of
0.2.................................................................................................................................................130

x 

Figure 4.19 One dimensional result compared with two-dimensional result for aspect ratio of
1…………………………………………………………………………………………………131
Figure 4.20 One-dimensional result compared with two-dimensional result for aspect ratio of
10………………………………………….…………………………………………………….132
Figure 4. 21 One-dimensional result compared with two-dimensional result for aspect ratio of
50..................................................................................................................................................133
Figure 4. 22 One-dimensional result compared with two-dimensional result for aspect ratio of
100………………………………………...…………………………………………………….134
Figure 4.23 A two-dimensional axisymmetric model ……………………………...…………..135
Figure 4.24 Angular integration boundary of Eqn. 4.75………………………………………..136
Figure 4.25 Comparison of peridynamic calculation of strain at Point A using axisymmetric
model with results from the vibration theory ………………………………………………..…137
Figure 4.26 Comparison of peridynamic calculation of strain at Point B using axisymmetric
model with results from the vibration theory …………………………………………..………138
Figure 4.27 Comparison of peridynamic calculation of strain at Point C using axisymmetric
model with results from the vibration theory ………………………………..…………………139
Figure 4.28 Calculation of the critical stretch…………………………………………………..140
Figure 4.29 Schematic drawing of the sample for crack propagation test……………………...141
Figure 4.30 Comparison of crack propagation simulated by peridynamics with experiment
result ……………………………………………………………………………...…………….142
Figure 4.31 Comparison of crack propagation speed simulated by peridynamics with experiment
result ……………………………………………………………………………………...…….143
Figure 4.32 A plate with single notch for crack branching study………………………………144
Figure 4.33 Computational crack path (a) damage and (b) strain energy density simulated by
peridynamics …………………………………………………………………………………...145
Figure 4.34 A two-notch plate under impulsive loading……………………………………….146
Figure 4.35 Computational crack path of the top notch ………………………………………..147
Figure 4.36 Horizontal displacement at 33.2

from computational results………………......148

xi 

Figure 5.1 A composite plate with fiber in ° direction and a bond in °
direction ………………………………………………………………………………………..173
Figure 5.2 Coordinate system for calculating strain energy density at point

…………...…...174

Figure 5.3 Beam model for orthotropic materials………………………………………………175
Figure 5.4 Pulling test in a composite laminate ………………………...……………………...176
Figure 5.5 Comparison of
of 0° laminate calculated from peridynamics (top) and composite
theory (bottom)……………………………………………………………………………..…..177
of 45° laminate calculated from peridynamics (top) and the
Figure 5.6 Comparison of
composite theory (bottom)…………………………………………………………….………..178
Figure 5.7 Comparison of
of 60° laminate calculated from peridynamics (top) and the
composite theory (bottom)………………………………………………………………….…..179
Figure 5.8 Comparison of
of 0° laminate calculated from peridynamics (top) and the
composite theory (bottom)……………………………………………………………….……..180
Figure 5.9 Comparison of
of 45° laminate calculated from peridynamics (top) and the
composite theory (bottom)…………………………………………………………….………..181
Figure 5.10 Comparison of
of 60° laminate calculated from peridynamics (top) and the
composite theory (bottom)……………………………………………………….……………..182
Figure 5.11 Free vibration of a laminated beam with fibers in x direction…………………….183
Figure 5.12 Peridynamic results compared with composite beam theory results with
5.
Top: horizontal displacement of point A. Bottom: vertical displacement of point B …….184
Figure 5.13 Peridynamic results compared with composite beam theory results with
Top: horizontal displacement

of point A. Bottom: vertical displacement

20.

of point B ….....185

Figure 5.14 An unidirectional composite plate with single edge notch under three point
bending………………………………………………………………………………………….186
Figure 5.15 Comparison of crack propagation velocity between peridynamics and
experiment………………………………………………………………………………………187
Figure 5.16 Compact tension test for a unidirectional composite
plate …………………………………………………………………………………………….188

xii 

Figure 5.17 Simulation results for
50	
and
0°: (a) vertical displacement, (b) strain
energy density, (c) local damage……………………………………………………………….189
Figure 5.18 Figure 5.18 Simulation results for
70	 and
30°: (a) vertical displacement,
(b) strain energy density, (c) local damage………………………………...…………………...190
Figure 5.19 Simulation results for
70	 and
45°: (a) vertical displacement, (b) strain
energy density, (c) local damage………………………………………………..……………...191
Figure 5.20 Simulation results for
90	 and
60°: (a) vertical displacement, (b) strain
energy density, (c) local damage………………………………………………………..……...192
Figure 5.21 Simulation results for
100.5	 and
90°: (a) vertical displacement, (b)
strain energy density, (c) local damage…………………………………………………..……..193
Figure 5.22 Simulation results for
40	 and
30°: (a) vertical displacement, (b) strain
energy density, (c) local damage……………………………………………..………………...194

xiii 

Chapter 1
Introduction

1.1 Background
The governing equation of continuum mechanics is a partial differential equation as given in the
following equation
					

(1.1)

All solutions to solid mechanics problems, analytical as well as numerical, are based on this
governing equation. The divergence of stresses in Eqn. 1.1 requires that the displacement field to
be continuous. This requirement cannot be continuously satisfied when cracks occur since new
cracks come with new boundaries which require new boundary conditions for solving the partial
differential equation. Numerical techniques based on continuum mechanics such as finite
element method are commonly used for investigating dynamic crack propagation in structures
with complex geometry. However, prior knowledge of crack position and orientation is required.
The simulation of dynamic crack propagation can become even more challenging, if not
impossible, for fiber composites due to their high inhomogeneity and high anisotropy.

Silling [1] proposed an alternative formulation, so-called peridynamics. Its equation of motion
can be expressed as follows
	
where

is mass density,

is acceleration,

(1.2)
is bond force and

is body force. The bond force

is a function of displacements. The integration in Eqn. 1.2 is applied to the entire domain of

1

study. In peridynamics, the domain of study is assumed to be organized by points. Each point is
connected to all points within its horizon. Eqn. 1.2 can be applied to continuous as well as
discontinuous domains such as those with cracks. This ability renders peridynamics more useful
for progressive damage analysis than commonly used finite element method which requires
constant remeshing of the domain. Besides, peridynamics does not require any additional
damage theory such as fracture mechanics since it is essentially based on the bond strength
between points.

1.2 Literature review
The peridynamic theory proposed by Silling [1] assumes that all points in a domain of interest
are connected by bond forces. A two-point bar model is used in Ref. [1] and the bond force
between the two points is dependent on the positions of the two points. The properties of the bar
are related to the properties of the material of the domain. A shear crack propagation is also
presented in Ref. [1] to show the benefit of the peridynamic theory. Silling [2] mentioned that
peridynamics is similar to molecular dynamics in accounting of all points within a horizon of a
point for calculating the internal forces of the point. However, the peridynamic domain is
continuous and the grid can be of any size based on numerical consideration. This is a major
difference from molecular dynamics in which the points can only be molecules or atoms. Hence,
the molecular dynamics cannot be used to solve any solid mechanics problem in reasonable time
due to the overwhelmingly large number of molecules involved in any structure. The bar-based
peridynamic model is similar to the models proposed by Kunin [3] and Rogula [4] in the
investigation of continuum properties of crystals based on interatomic interactions.

2

Several one-dimensional peridynamic solutions can be found in Refs. [5-8]. The studies include
the effects of infinite horizon and comparisons of peridynamic solutions with those from
continuum mechanics. Silling, Simmermann and Abeyaratne [5] studied an infinite long bar
subject to two concentrated loads. They found the displacements away from the loading points
are the same between peridybnamic analysis and classical mechanics. However, the
peridynamics-based displacements oscillate when points are close to the loading. Weckner and
Abeyaratne [6] studied dynamic deformation of an infinite bar with various initial functions.
They found the influence of horizon on the results. The displacements were discontinuous even
though they started with continuous functions. However, the displacements converged to the
classical mechanics solution when horizon approached zero. Silling [9,10] discussed the
numerical implementations of peridynamics and EMU code [9]. Bobaru, Yang, Alvels, Silling,
Askari and Xu [11] discussed the convergence issues in a one-dimensional static problem. They
investigated an infinitely long bar with two concentrated forces by using analytical method and
then compared the solution with that from numerical solution. In peridynamic solution,
displacement was not continuous at points where forces were applied. There were singularity
points (infinite jump). Three convergence techniques were studied in the paper: (1) fixing
horizon and increasing the number of nodes in each horizon, (2) decreasing horizon and
increasing the number of nodes per horizon and (3) decreasing the horizon and fixing the number
of nodes per horizon. In case (1), the numerical solution converged to the analytical solution with
a finite displacement jump. In case (2), the numerical solution converged to the peridynamic
solution with singularity. In case (3), singularity disappeared when the horizon became zero. The
peridynamic solution converged to the solution of classical mechanics when the horizon reduced
to zero. However, it is almost impossible to have a zero horizon since the grid size will become

3

zero. Adaptive refinement of mesh was also discussed to improve the computational efficiency
by Silling [12].

Applications of peridynamics for isotropic materials can be found in Refs. [2,9,12-15]. Silling
and Bobaru [2] applied peridynamics to one-dimensional fiber structures and two-dimensional
membrane structures. They investigated central crack propagation, burst of balloon, and tearing
of rectangular sheet. They also simulated deformation of a fiber with different parts of the fiber
interacting with each other through van der Waals force. Gerstle, Sau and Silling [13] studied
dynamic crack propagation without initial crack in a two-dimensional plate. Silling [9] used
EMU code to simulate crack propagation in a two-notch plate under drop weight. These studies
showed the advantage of using long-range force in peridynamics for dealing with contact issues.
Weckner, Askari and Xu [14] simulated three-dimensional crack branching in a unitized metal
structure. Bobaru [15] studied a three-dimensional nano-fiber network using a similar method as
that used in Ref. [2]. It was found that there was no need of a separate failure theory in
peridynamic analysis.

The application of peridynamics on fiber-reinforced composite materials and the validation of
peridynamics with experiments can be found in Refs. [16-18]. Silling, Xu and Askari [16,17]
modeled fiber-reinforced composite materials by using two kinds of bonds, fiber bonds and
matrix bonds. They then studied central crack and delamination problems. Xu, Askari and
Weckner [19] compared the results from Refs. [16,17] with experiment results. Bobaru and
Silling [18] modeled composite materials with micro-structure based on the method mentioned in
Ref. [2].

4

The Poisson’s ratio in the bar model introduced by Silling [1] and Zimmerman [20] was a fixed
value. Gerstle, Sau and Silling [21], however, included pairwise moments in a model to account
for various Poisson’s ratios. A state-based peridynamic model was proposed in Refs. [22-26]. In
this model, there were interactions between two points in two parts. The first part depended on
the positions of the two points, similar to the bond-based model. The second part depended on all
points in the horizon. Warren, Silling, Askari, Weckner, Epton and Xu [23] studied damage with
state-based peridynamics. Silling and Lehoucq [25] defined a stress tensor in peridynamics.
When the horizon goes to zero, the stress tensor converged to Piola-Kirchhoff stress tensor.
Demmie and Silling [24] simulated gas using only the second part of the state-based peridynamic
model.

1.3 Motivations and objectives
Computational methods, such as finite element method, based on Eqn. 1.1 require a continuous
displacement field. The information concerning crack location and dimension also need to be
known and included in the numerical model prior to computation to avoid imposing exhausted
boundary conditions and crack surfaces everywhere in the domain of study. These prerequisites
are almost impossible for problems involving dynamic crack propagation in fiber-reinforced
composite materials which can have very complex damage configurations. The integral equation
based peridynamic theory, on the contrary, does not have any prerequisite since it is essentially
based bond force between points. The modeling of the damage process of fiber composites can
be benefited from peridynamics since crack can occur whenever bonds break and there is no
need of any additional damage theory other than the bond strength.

5

As a novel theory, peridynamics is still in its infant stage. More advanced models are needed.
Analytical verifications and experiment validations are also required for justifying the models.
Silling’s bar model [10] assumes that two points are linked by a bar and there is only one
material property – the modulus - and the Poisson’s ratio is limited to be ¼ in three-dimensional
problems. Silling, Xu and Askari’s fiber-reinforced composite material model [16,17] is also
based on the bar model. There are only two independent material properties in their models
although there are four independent material properties for two-dimensional orthotropic
materials and six independent material properties for three-dimensional orthotropic materials.

A new beam model for two-dimensional problems with two independent material properties will
be presented in this dissertation research. This model should be capable of simulating isotropic
materials with various Poisson’s ratios. An orthotropic material model will also be proposed
based on the new beam model. There will be four independent material properties in this model
to accommodate the four independent material properties in orthotropic materials. With this
material model and bond failure criterion, it should be sufficient to simulate dynamic damage
crack propagation in fiber composite materials.

1.4 Outline of dissertation
The basic formulation of peridynamics is introduced in Chapter 2. In Chapter 3, one-dimensional
equation of motion is investigated. A simple bar model is proposed and the associated material
property is identified from energy method. The convergence of the numerical method is then
studied. The one-dimensional model can be used to study split Hopkinson’s pressure bar (SHPB)

6

for verification and validation of the bar model. The bar model can also be used for studying
two-dimensional problems with a fixed Poisson’s ratio. In Chapter 4, a beam model is proposed
to simulate two-dimensional isotropic materials. The beam model is able to capture various
Poisson’s ratios. The verification of numerical studies requires a comparison with an analytical
solution based on classical mechanics. A failure theory is then established with the use of energy
method. The failure theory is subsequently used for investigating three problems concerning
dynamic damage propagation. The results are validated by experiment results from literature. A
model for anisotropic materials is proposed in Chapter 5, in which there are four independent
material properties. This model is verified by a static problem and a dynamic vibration problem.
The dynamic damage propagation in composites is then simulated by the proposed fourparameter model and the results are validated by experiment results from published works.

7

REFERENCES

8

References
[1] Silling, S. A., Reformulation of elasticity theory for discontinuities and long-range forces, J.
Mech. Phys. Solids, 48, 175-209 (2000)
[2] Siling. S. A., Bobaru, F., Peridynamic modeling of membranes and fibers. Int. J. Non-linear
Mech., 40, 395-409 (2005)
[3] Kunin, I. A., Elastic Media with Microstructure I: One-dimensional Models, pp. 27. Springer,
Berlin (1982)
[4] Rogula, D., Nonlocal Theory of Material Media, pp. 137-149 and 243-278. Springer, Berlin
(1982)
[5] Siling, S. A., Simmermann, M., Abeyaratne, R., Deformation of a peridynamic bar. J. Elast.,
73, 173-190 (2003)
[6] Weckner, O., Abeyaratne, R., The effect of long-range forces on the dynamics of a bar. J.
Mech. Phys. Solids, 53, 705-728 (2005)
[7] Weckner, O., Emmrich, E., Numerical simulation of the dynamics of a nonlocal,
inhomogeneous, infinite bar. J. Comput. Appl. Mech., 6, 311–319 (2005)
[8] Meerich, E., Weckner, O. Analysis and Numerical Approximation of an Integro-differential
Equation Modeling Non-local Effects in Linear Elasticity, Mathematics and Mechanics of Solids,
12, 363-384, 2007
[9] Silling, S.A., Dynamic fracture modeling with a meshfree peridynamic code, Computational
Fluid and Solid Mechanics, 641–644, Elsevier, Amsterdam (2003)
[10] Silling, S.A., Askari, E., A meshfree method based on the peridynamic model of solid
mechanics. Comput. Struct., 83, 1526–1535 (2005)
[11] Bobaru, F., Yang, M., Alvels, L. F., Silling, S. A., Askari, E. Xu, J., Convergence, Adaptive
Refinement, and Scaling in 1D Peridynamics, International Journal for Numerical Methods in
Engineering, 77, 852-877 (2009)
[12] Bobaru, F., Alves, L., Silling, S., Askari, A., Dynamic Crack Branching and Adaptive
Refinement in Peridynamics, 8th World Congress on Computational Mechanics, June 30-July 5,
2008, Venice, Italy
[13] Gerstl, W., Sau, N., Silling, S., Peridynamic modeling of plain and reinforced concrete
structures. In: 18th International Conference on Structural Mechanics and Reactor Technology
(SMiRT 18), Beijing, China, SMiRT18-B01-2(2005)

9

[14] Weckner, O., Askari, A., Xu, J., Damage and Failure Analysis Based on PeridynamicsTheory and Applications, 48th AIAA/ASME/ASCE/AHA/ASC Structures, Structural Dynamics,
and Materials Conference, April 23-26, 2007, Honolulu, Hawaii
[15] Bobaru, F., Influence of van der Waals Forces on Increasing the Strength and Toughness in
Dynamic Fracture of nanofiber Networks: a Peridynamic Approach, Modeling and Simulation in
Materials Science and Engineering, 15, 397-417 (2007)
[16] Silling, S. A., Xu, J. Askari, E., Peridynamic modeling of impact damage. In: Moody, F.J.
Problems Involving Thermal-Hydraulics, Liquid Sloshing, and Extreme Loads on Structures,
489, 197-205, American Society of Mechanical Engineers, New York (2004)
[17] Askari, A., Xu, J., Silling, S., Peridynamic analysis of damage and failure in composite. In:
44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, AIAA-2006-88 (2006)
[18] Bobaru, F., Silling, S.A., Peridynamic 3D models of nanofiber networks and carbon
nanotub-reinforced composites. In: Materials Processing and Design: Modeling, Simulation, and
Applications– NUMIFORM 2004 – American Institute of Physics Conference Proceedings, 712,
1565–1570 (2004)
[19] Xu, J., Askari, A., Weckner, O., Damage and Failure Analysis of Composite laminates under
Biaxial Loads, 48th AIAA/ASME/ASCE/AHA/ASC Structures, Structural Dynamics, and
Materials Conference, April 23-26, 2007, Honolulu, Hawaii
[20] Zimmermann, M., A Continuum Theory with Long-Range Forces for Solids, PhD
dissertation, Massachusetts Institute of Technology (2005)
[21] Gerstle, W., Sau, N., Silling, S., Peridynamic modeling of concrete structures, Nucl. Eng.
Des., 237, 1250-1258 (2007)
[22] Silling, S., Epton, M., Weckner, O., Xu, J., Askari, E., Peridynamic States and Constitutive
Modeling, Journal Elasticity, 88, 151-184 (2007)
[23] Warren, T., Silling, A. S., Askari, Weckner, O., Epton, A., Xu, J., A Non-ordinary Statebased Peridynamic Method to Model Solid Material Deformation and Fracture, International
Journal of Solids and Structures, 46, 1186-1195 (2009)
[24] Demmie, P., Silling, A., An Approach to Modeling Extreme Loading of Structures Using
Peridynamics, Journal of Mechanics of Materials and Structures, 2, No. 10 (2007)
[25] Bobaru, F., Yang, M., Alvels, L. F., Silling, S. A., Askari, E. Xu, J., Convergence, Adaptive
Refinement, and Scaling in 1D Peridynamics, International Journal for Numerical Methods in
Engineering, 77, 852-877 (2009)
[26] Lehoucq, R., Silling, S., Force Flux and the Peridynamic Stress Tensor, Journal of the
Mechanics and Physics of Solids, 56, 1566-1577 (2008)
10

Chapter 2
Basic Formulations

2.1 Peridynamic model of a continuum
The acceleration of any particle at 	in the reference configuration at time 	can be expressed by
the following equation [1]
																					
where

,

,

is displacement,

, ,

	

is body force density,

square which exerts on particle

by particle

,

is mass density and

(2.1)
is force per volume

. Fig. 2.1 shows the relation.

The relative position of these two particles in the reference configuration is denoted by	 , i.e.
																																																																											

(2.2)

while the relative displacement is denoted by	 , i.e.
																																																																						

Accordingly,

,

,

(2.3)

represents the current relative position vector between the two particles in

the current configuration as shown in Fig. 2.2.

It is convenient to assume that there exists a positive number , called horizon, also shown in
Fig. 2.1, such that the force

vanishes when | | is greater than , i.e.

																																																													| |

⇒	

,

0								∀

The domain of the integration in Eqn. 2.1 should be the spherical neighborhood of
11

(2.4)

with a

radius equal to	 .

Bases on the conservation of angular momentum, the following equation should hold
																																																							

,
,

That is
where

,

,

0								∀	 ,
∙

/|

(2.5)
|

(2.6)

is a scalar function. This is to say, the force vector between the two particles is

parallel to their current relative position vector.

The natural boundary condition and external forces are applied through the body force density	 .
They can be made nonzero within a boundary layer. The thickness of the boundary layer is
commonly assumed to be equal to the horizon .

2.2 Constitutive modeling
If the material of interest is isotropic, 	from Eqn. 2.6 does not depend on the direction of	 . For
simplicity, it can be assumed that the scalar bond force

depends only on the bond stretch,	 ,

defined by
																																											

|

| | |

(2.7)

| |

where is positive when the bond is in tension and negative when the bond is in compression.

Accordingly,
																																										

∙

(2.8)

where is stiffness which depends on both material and the geometrical property of the material.

12

A material is said to be micro-elastic [1] if the force function is derivable from a scalar micropotential

:

																																											

,

,

																							∀	 ,

(2.9)

The micro-potential is the energy in a single bond and has the dimension of energy per unit
volume square. The energy per unit volume at a given point, i.e. the strain energy density, then
becomes
																																										

,

	

(2.10)

The domain of integration is the same as that of Eqn. 2.1. The factor of ½ is imposed because
each point of a bond owns only half the energy in the bond.

2.3 Failure criterion
The bond between two particles will break if they are stretched beyond a critical value	

defined

as the critical stretch. Once the bond fails, there is no tensile force between the two points. The
force function

,

can be modified as history dependent:

																																														

,

∙ ∙

,

(2.11)

where 	is a history-dependent scalar function, i.e.
,

	

,
1								 	
0				 	
,

The damage at a point can be defined by a function	
																																														

,

1

, 	

		
		

	
	

	0
	0

′
′

	
	

(2.12)

, , where
(2.13)

13

The domain of integration is the same as that of Eqn. 2.1. Note that 0

1,	 where 0

represents virgin material and 1 represents complete disconnection of the point from all points
with which it initially interacted. A value of

should be defined as the criterion of the beginning

of the crack (discontinuity in displacement). A value of 0.3 is used in [2].

2.4 Numerical method
The force function

is usually very complex for integration. Numerical integration hence should

be used. To begin with, the integral form of Eqn. 2.1 is replaced by a finite summation
∑

,

where superscript n is the integration step number,		

(2.14)
is the volume of node p and subscripts

denotes the node number as shown in Fig. 2.3. The summation is taken over all nodes p such that
																																																|

After obtaining

|

(2.15)

at each point, the displacement of the next integration step

can be found

from the formula for finite difference, i.e.

									

∆

∆

∆

∆

Fig. 2.4 shows the flow chart of the numerical program.

14

(2.16)

Figure 2.1 Horizon of point .

15

Figure 2.2 Reference configuration and current configuration.

16

Figure 2.3 Finite summation domain in numerical computation.

17

Figure 2.4 Flow chart of a peridynamic program.

18

REFERENCES

19 
 

References

[1] Silling, S. A., Reformulation of elasticity theory for discontinuities and long-range forces, J.
Mech. Phys. Solids, 48, 175-209 (2000)
[2] Silling, S.A., Askari, E., A meshfree method based on the peridynamic model of solid
mechanics. Comput. Struct, 83, 1526–1535 (2005)

20 
 

Chapter 3
One-Dimensional Wave Propagation Analysis

3.1 Identification of stiffness
Based on the peridynamic theory established in Chapter 2, this chapter investigates onedimensional wave propagation. To obtain the peridynamic equation of motion for the onedimensional problem, consider a bar with a length L and a constant cross-sectional area A. Eqn.
2.1 is integrated over the cross-section and divided by A, resulting in
																																																				

,

,

,

′

																																																											(3.1)

,

where

≡

and

′∙

unit length per unit volume, e.g.

(3.2)

,
,

,
,

(3.3)

,

∬

,

,

′

(3.4)

′ are used. The bond force f is expressed in terms of force per
/

.

A simple model of the bond force can be expressed as
∙

(3.5)

where c is stiffness and is strech. To find , it is necessary to find the strain energy density from
peridynamics and set it equal to that used in classical mechanics.

Consider a bar under a uniform stretch of . Based on the bond micro-potential of Eqn. 2.9,
integrate f to find

21

∙| |

																																																									

∙

∙| |

(3.6)

∙| |

This is the energy in one bond. From Eqn. 2.10, the strain energy density based on peridynamic
analysis becomes

																															

∙

	 

The strain energy density from classical mechanics is

∙|  |

	 

∙

∙

(3.7)

. Setting it equal to that obtained

from peridynamic analysis, i.e. Eqn. 3.7, the stiffness c in Eqn. 3.5 can be found to be
(3.8)

3.2 Definition of stress in peridynamics
One way to verify a peridynamic model is to compare it with the associated term in classical
mechanics. Most classical mechanics results are presented by stresses. Therefore, to verify
peridynamics, it is useful to define stress in peridynamics and compare it with the classical
mechanics counterpart. This definition is for comparison purpose only since stress is not
commonly used in peridynamics. As shown in Fig. 3.1, the stress at a point

is defined as the

force connecting the two half-domains together at point , i.e. the summation of all positive bond
forces crossing the point .

∙ 	

	

∙

(3.9)

where stiffness is defined in Eqn. 3.8. Hence, the stress in peridynamics, shown in Eqn. 3.9, is
exactly the same as that in classical mechanics.
22

3.3 Comparison with analytical solution
Consider the wave propagation in a one-dimensional bar due to an initial uniform strain

, as

depicted in Fig. 3.2. Both ends of the bar are assumed to be free of load.

3.3.1 Analytical solution
The equation of motion for the bar [1] is

0
where u is displacement, t is time, E is Young’s modulus,

(3.10)
is density and x is position.

The boundary conditions are

0,

(3.11)

,
where

0
0

(3.12)

is the length of the bar.

The initial conditions are
,0
,0

(3.13)
0

(3.14)

Using the method of separation of variables, the displacement can be expressed as a product of
position

and time , i.e.
,

(3.15)
23

Substituting Eqn. 3.15 into Eqn. 3.10, it yields
(3.16)

Since the left-hand side of the equation is dependent on only and the right-hand side of the
equation is dependent on x only, both sides of Eqn. 3.16 then must be equal to a constant, such as
.

Rewrite Eqn. 3.16 as
																																																																									

(3.17)

0
where

0

(3.18)

cos

(3.19)

.

A general solution of Eqn. 3.17 is
sin

0. Similarly, a general solution of Eqn.

In order to satisfy the initial condition of Eqn. 3.14,
3.18 is
sin

cos

In order to satisfy the boundary condition given in Eqn. 3.11,
the boundary condition given in Eqn. 3.12, sin

0. Hence,

24

(3.20)

0. Similarly, in order to satisfy

,

0, 1, 2, 3, …

,

0, 1, 2, 3, …

Redefine
(3.21)

The solution to Eqn. 3.10 can be expressed as
∑

,
where

cos

in Eqn. 3.19 has been replaced by

cos

(3.22)

0, 1, 2, 3, … and a new coefficient

(

has

been installed.

Eqn. 3.13 can be rewritten as
∑

In order to find

cos

(3.23)

in Eqn. 3.23, multiply both sides of the equation with cos

and integrate

them from 0 to , i.e.
cos

∙ cos

	

cos

	

(3.24)

Hence,

1
∑

1
∑

1
1

1

cos
1

25

(3.25)

∙ cos
sin

∙ cos

(3.26)

(3.27)

in which

∙

	 .

and

This analytical solution will be used to justify numerical solutions based on peridynamics. It will
also be used to study the convergence of parameters used in peridynamic analysis. Fig. 3.3 shows
the strain history of a bar under an initial strain of 0.0001. The bar is made of a material with
Young’s modulus

193	

8027	

and density

/

. The ripples at the corners are

due to Gibbs phenomenon [2], which states that the truncated Fourier series of a discontinuous
signal will exhibit high frequency ripples and overshoot near the discontinuities. The overshoot
does not die out as the frequency increases, but approaches a finite limit.

3.3.2 Numerical studies
Three parameters are involved in the peridynamic numerical analysis. They are the time step
the size of horizon

and the number of nodes within the horizon

,

. Fig. 3.4 shows a one-

dimensional peridynamic model and the associated horizon. When solving the peridynamic
problem and comparing its results with those obtained from classical mechanics analysis, it is
found that the magnitudes of waves from both analyses are always the same. However, the wave
speed and wave period vary with
combinations of

and . Table 3.1 shows the wave periods based on nine

and .

Based on theory of vibration, the wave period in the bar is equal to
2
where L is the bar length and

/

0.41	

/ 	 is wave speed.
26

(3.28)

Table 3.1 Wave period (

) based on various of

1 mm

5 mm

10mm

1

0.2912

0.2893

0.2873

5

0.3761

0.3747

0.3732

10

0.3827

0.3915

and

0.3903

	

Comparing the peridynamic results given in Table 3.1 with that from vibration theory, the
numerical solution becomes closer to the analytical solution when the number of nodes within
the horizon ( ) increases, i.e., the distance between nodes decreases. Besides, when

is fixed,

the peridynamic solution becomes closer to the classical mechanics solution if the horizon
decreases. This result is thought to be due to the fact that classical mechanics is based on contact
force. When

decreases, the bond force in peridynamics becomes closer to the contact force in

classical mechanics and the peridynamics solution converges to the classical mechanics solution.
It then seems to be ideal to use small

and large

in peridynamic analysis. However, this will

pose a significant increase in computational time. Fig. 3.5 shows comparison between classical
mechanics and peridynamic analysis. Numerical studies based on
0.5	

2.2	

and

(with m = 4) give peridynamic results within 1.2% of discrepancy from that of classical

mechanics. The convergence of time step based on these parameters is also shown in Fig. 3.6. As
long as

and

Consequently,

are fixed, the result does not vary too much with the change of time step.
1

10

,

0.5	

and

numerical investigations.

27

2.2	

will be used in subsequent

3.4 Comparison with experiment results
For studying one-dimensional wave propagation, the commonly used split Hopkinson’s pressure
bar (SHPB) [3] is taken as an example in this study. SHPB consists of three bars, a striker bar, an
incidence bar and a transmission bar as shown in Fig. 3.7. It is commonly used for characterizing
material constitutive relations. Table 3.2 shows an example of dimensions and material
properties of bars used in SHPB. The results from the experiments based on SHPB can be used
to justify the peridynamic analysis.

Table 3.2 Dimensions and material properties of the SHPB in Fig. 3.7
Striker

Incidence bar

Transmission bar

Length

191 mm

991 mm

724 mm

Diameter

13 mm

13 mm

13 mm

Young’s modulus

193 GPa

193 GPa

193 GPa

Mass density

8027

/

8027

/

8027

/

In performing the SHPB test, a specimen is placed in between the incidence bar and the
transmission bar. The striker is then accelerated to collide with the incidence bar, creating a strain
wave. The strain wave will propagate through the two bars and rebound from any surface
between any two bars. A pair of strain gages are installed on each of the incidence bar and
transmission bar to record the strain wave history. From the strain wave, the dynamic stressstrain curve of the specimen can be identified.

When the striker bar is fired onto the incidence bar, it creates a constant pressure over a period of

28

time. The strain wave then propagates toward the specimen with partial transmission into the
specimen and partial reflection from the specimen, which has different material properties and
cross-sectional areas than the split bars as shown in Fig. 3.8. The reflected strain

	is related to

the strain in the specimen 	by the following relation
/

where

,

and

(3.29)

are Young’s modulus, mass density and length of the incidence bar,

respectively. The transmitted strain

is related to the stress in the specimen

by the following

relation
(3.30)
where

and

are Young’s modulus and cross-sectional area of the transmission bar and

is

the original cross-sectional area of the specimen. The stress-strain relation of the material can
then be found from Eqn. 3.29 and Eqn. 3.30.

The goal of this study is to validate peridynamic models and associated computational results
with experiment results from SHPB tests, i.e. validation of theory by experiment. Hence,
materials with known properties will be used as input for peridynamic analysis. The strain
history in the incidence bar and the transmission bar will be computed based on peridynamics
and the results will be compared with the experiment results.

3.4.1 Square wave input
A perfect impact between a striker bar and an incidence bar can result in a square-shaped pulse
with a stress magnitude

and a duration
29

2 , where , ,

and

are velocity

before impact, mass density, wave velocity and length of the striker bar. In peridynamic analysis,
the constant pressure

2 , shown in Fig. 3.9, are used as input in

and duration

the computational program.

Study One - one bar
In this study, no specimen is used and the transmission bar is isolated from the incidence bar,
resulting in only wave propagation in the incidence bar. Experiment results and computational
results are shown in Fig. 3.10 for comparison.

Study Two - two bars
In this study, no specimen is used and the incidence bar and the transmission bar are put in
contact with each other as shown in Fig. 3.11. The testing results are compared with peridynamic
simulation as shown in Fig.3.12. They are the strain histories of the incidence bar and the
transmission bar. One obvious observation is that the majority of the wave enters the
transmission bar with only a small portion of approximately 12% reflected. At the end of wave
propagation, the two bars separate from each other.

Based on theoretical analysis, the wave period in the transmission bar is

																																																							2
where

	

2.95

is the length of transmission bar and

measurement gives an average wave period of 2.99

10

(3.31)

is the wave velocity. The experiment
10

, which is 1.4% higher than the

theoretical value, while the peridynamic simulation gives an average of 3

30

10

s, which is

1.7% higher than the theoretical value.

In the peridynamic simulation of SHPB, the bond force crossing the interface between two bars
can only be negative, otherwise the two bars will separate. At each time step of calculation, the
stretch of the bond between the two bars is calculated. If the stretch is greater than zero, it
implies there is a separation between two bars. Hence, the bond force is set to zero. When the
bond stretch is negative, the bond force can then be calculated by Eqn. 3.5. Using this technique,
there will be no interaction between the incidence bar and the transmission bar when they are
separated. The comparison between the simulation and experiment for the transmission bar
shows good agreement. The periods also match well even after several periods. However, the
comparison for the incidence bar shows a discrepancy. Oscillations on the incidence bar do not
look like reflection waves. They are smaller than the reflection wave shown in Fig. 3.12.

In the experiment of SHPB, a layer of grease may be applied to the interface between two bars to
promote the continuity of wave propagation. The effect of the grease layer on the wave
propagation should then be considered in the investigation. The testing result given in Fig. 3.13
shows almost no reflection in the incidence bar, implying the complete wave transmission
through the grease layer to the transmission bar.

Study Three – specimen with different property
The wave propagation through an interface can be affected by the material change and the crosssectional area change across the interface. Both the changes can cause wave reflection on the
interface. This study investigates the effect of material change on the wave propagation. A study
31

model given in Fig. 3.14 is used for the study. In identifying the interfacial bond force, it is
necessary to determine a composite stiffness c, given in Eqn. 3.5 and a composite cross-sectional
area A, given in Eqn. 3.2. The composite stiffness c can be expressed in terms of the individual
stiffness on

and

on both sides of the interface with the coordinate

and the displacement

of the interface.

The bond force for material with stiffness equal to

is defined as

																																																																	
while that for material with stiffness

(3.32)

is defined as

																																																																

(3.33)

Both of them should be equal to the composite bond force, i.e.

																																																									

(3.34)

Solving for , it yields

																																																							

(3.35)

In order to demonstrate the effect due to the change of material across an interface, an aluminum
specimen with a cross-sectional area identical to that of SHPB (Fig. 3.15) is used. The specimen
is made of AL6061 with a Young’s modulus of 69 GPa and a mass density of 2700

32

/

.

Results are shown in Fig. 3.16 for comparison. The simulation results match reasonably well
with the experiment counterparts. The technique for modeling an interfacial bond joining
different materials is thus validated.

Study Four – specimen with smaller cross-sectional area
This study investigates the effect of different cross-sectional areas across an interface on the
wave propagation. A specimen made of Steel 347, identical to that used for the SHPB is used.
The diameter of the specimen is 10

and the length of it is 26

. Since the specimen

diameter is smaller than the incidence bar diameter, the bond stiffness of the specimen should be
modified with the area ratio between them, i.e.
	

(3.37)

Once the specimen bond stiffness is defined, the composite bond stiffness can be found from
Eqn. 3.35. The result is given in Fig. 3.18.

3.4.2 Using filter for computational results
All the four computational results given above match reasonably well with their experiment
counterparts. However, it is necessary to point out that the computational results always show
square waves with overshoots around the corners when a square wave as shown Fig. 3.9 is used
as input. The experiment results, however, show trapezoidal waves with smooth, low shoulders.
The overshoots are likely due to numerical process while the low shoulders are likely due to the
bandwidth limit of data acquisition system. The sharp rise of the square wave is of a very high
frequency signal. It is almost impossible to have a perfectly square wave in experiment since all
33

electrical devices work as low-pass filters and they have a limited bandwidth. Based on this
realty, a filter is added to the following computational investigations.

The Fourier transform for a discrete-time function [2] is
∑

(3.38)

The inverse Fourier transform, back to time domain, is
(3.39)
where ω is the frequency. If the high frequency part of function

is cut off and Eqn. 3.39 is

added to transform back to the time domain, a low-pass filter can be achieved. In computation,
Matlab function

	[4] shown below may be used
(3.40)

where

is the strain history vector and

Once a cutoff value is set, the first
inverse Fourier transform function

is the magnitude of the -th frequency term.

components of

can be used to find a new strain history by

[5], i.e
,

(3.41)

Simulation results for Study One, Study Two, Study Three and Study Four with filters are
compared with experiment results. Fig. 3.19, Fig. 3.20, Fig. 3.21 and Fig. 3.22 show the
comparisons. The cut-off frequency used is 100 MHz. As can be seen from the four sets of
comparison, the computational results with filers are of trapezoidal shape and the overshoots are
smaller. The numerical solutions match more favorably with the experiment results with the use
34

of filter.

3.4.3 Using the incidence wave as the computational input
Another source to cause the discrepancies between the computational simulations and the
experiment results is likely due to the imperfect impact between the striker bar and the incidence
bar. The computational study is based on a perfect impact while the experiment study is not. One
way to eliminate the discrepancy is to use the incidence wave from experiment as input for
computation. Take Study Two as an example, the incidence wave can be extracted from the
experiment shown in Fig. 3.23. It then can be used as the input wave for the computational study.
The computational result is given in Fig. 3.25. The computational result is almost the same as the
experiment result. With the incidence waves from experiments as the computational inputs,
computational simulations for Study One, Study Two, Study Three and Study Four are shown in
Fig. 3.24, Fig. 3.25, Fig. 3.26 and Fig. 3.27, respectively. The results between experiments and
computations match very well.

3.5 Simulation with a striker
All simulations given above use a mathematical stress function as the input pressure. The
function is either a perfect square function or an incidence wave from an experiment. The
simulations do not model the impact process between the striker bar and the incidence bar. To
simulate and validate the true impact process, an impact bar system shown in Fig. 3.28 is
considered. Only a striker bar and an incidence bar are needed for investigating the incident
wave created by the impact. Their dimensions and material properties are listed in Table 3.3. A
velocity sensor is installed right before the incidence bar to monitor the striker velocity right

35

before the impact.

Table 3.3 Dimensions and material properties of the split bar in Fig. 3.28
Striker

Incidence bar

Length

491

2743

Diameter

25.4

25.4

Young’s modulus

227.26

227.26

7610.5

Mass density

/

7610.5

/

The simulation begins when the striker just touches the incidence bar without causing any
deformation to the incidence bar. It is a two-bar problem with the striker bar having an initial
velocity while the incidence bar stays static.

From the velocity sensor, the striker has a velocity of 4.12

/ right before impact. The strain

history in the incidence bar after the impact is recorded by the strain gages mounted on it. In
simulation, the velocity of 4.12
the striker bar (491

/ is the only input for the computation, i.e. all nodes within

) have an initial velocity of 4.12

/ . Fig. 3.30 shows the comparison

between the experiment and the simulation. Results from simulation with a stress input, similar
to Study One, is also presented in Fig. 3.30. The two simulation results are exactly the same.

Another way to verify the aforementioned results is to study the impact of two bars using
classical mechanics. Consider a short cylindrical bar, with a length
and a mass density of

, a cross-sectional area of

, is travelling to a long cylindrical bar with a cross-sectional area of

36

and a mass density of

. Right before the impact, the short bar has a velocity of

and both

bars are stress free as shown in Fig. 3.31. After impact (Fig. 3.32), a compressive wave travels in
both bars. Two conditions must be satisfied at the interface between the two bars. The first
condition is the continuity of velocity. If the interface velocity after the impact is , the particle
velocity at the interface of the short bar is

. The following relation holds for

,

and

.
(3.42)

The particle velocity at the interface of the long bar is

, i.e.
(3.43)

Combining Eqn. 3.42 and Eqn. 3.43, it yields
(3.44)

Another condition which must be satisfied is the balance of forces at the interface, i.e.
(3.45)
where

and

are axial forces on the short bar and long bar, respectively. The relation between

the particle velocity and the axial force can be derived from wave theory as follows
(3.46)
where

,

and

.

(3.47)
where

and

and

.
37

Combining Eqn. 3.44, Eqn. 3.45, Eqn. 3.46 and Eqn. 3.47, the following relations can be
concluded
(3.48)
(3.49)

In SHPB, the striker bar is the short bar while the incidence bar is the long bar. The two bars
have identical mass density, Young’s modulus and cross-sectional area. From Eqn. 3.44, Eqn.
3.48 and Eqn. 3.49, it yields
(3.50)

Combining Eqn. 3.46, Eqn. 3.47 and Eqn. 3.50, the following stresses can be obtained
(3.51)
where

and

. Then the strain on the incidence bar becomes
/

Consider an SHPB system of
bars

25.4	

227.26	

(3.52)

,

7610.5	

/

and the diameter of the

. Also, from the velocity sensor, the striker velocity before the impact is 4.12

/ . Substituting these values into Eqn. 3.52, the magnitude of strain is found to be 3.77
10 . The magnitude of strain wave from computation is 3.79
0.5% exists.

38

10 . A difference of only

In conclusion, based on the assumption of perfect impact and perfect signal acquisition, the study
based on striker impact is the same as the square function input as shown in Fig. 3.9. Therefore,
there is no need to include the striker impact in the simulation, hence significant computational
effort can be save. Another observation is that peridynamics is convenient for simulating the
impact process without using the complex theory involved in contact mechanics. The long range
force covers the contact between two bodies automatically.

3.6 Validation with results in literature
The purpose of the previous study is to validate the elastic model of peridynamics. The initial
velocity of the striker bar is controlled to be small so the deformation of the specimen is always
within elastic range. In SHPB experiment, the deformations of the bars are always within elastic
range.

Most of the existing work on SHPB is to find material properties of a specimen, i.e. stress-strain
curve, in elastic as well as in plastic range. To simulate an SHPB test from the existing works, a
plastic or nonlinear bond function must be used for the specimen.

An SHPB system is given in Fig. 3.33. The dimensions and material properties are listed in Table
3.4. An aluminum specimen with a length of 25.38	

and a diameter of 12.6	

is used in

the SHPB study. From the reflection wave and the transmission wave, the stress-strain relation of
the specimen can be identified and is shown in Fig. 3.35.

39

Table 3.4 Dimensions and material properties of the SHPB in Fig. 3.33
Striker

Incidence bar

Al Specimen

Transmission bar

Length

203 mm

838 mm

25.38 mm

419 mm

Diameter

19 mm

19 mm

12.6 mm

19 mm

191.4 GPa

Fig. 3.35

191.4 GPa

Young’s modulus 191.4 GPa
Mass density

7858	

/

7858

/

2700

/

7858

/

Fig. 3.35 is then used as the input material property for peridynamic simulations. The
computational results are then compared with the experiment results from Ref. [6]. Matlab
curve-fitting toolbox [7] is used to convert the curve in Fig. 3.34 for a piece-wise linear function.
For every strain, the corresponding stress calculated from the fitted function is divided by the
strain to obtain an equivalent Young’s modulus which is subsequently used to find the bond
stiffness. Fig. 3.35 compares experiment results with computational results. They show a good
correlation except some discrepancies around the corner. Fig. 3.36 shows the filtered
computational results with the experiment results. The agreement between the two is
significantly improved. This study demonstrates the capability of peridynamics in modeling
plastic materials.

3.7 Modeling shaper in SHPB
In performing SHPB tests, it is necessary to maintain constant strain rate during the loading
process. This requires a slow loading process such as a triangular incidence wave as shown in
Fig. 3.37. A square incidence wave, as shown previously for computational purpose, can cause a

40

rapid loading in the beginning of testing and hence is not suitable since it can damage the
specimen prematurely. A so-called shaper must be installed between the striker bar and the
incidence bar as shown in Fig. 3.38 to slow down the initial loading rate.

The selection of the shaper materials and its dimensions, however, requires significant
experience and trial and error. A computer simulation based peridynamics may help to ease the
selection process.

Frew [8] used a copper shaper in SHPB tests. To study the shaping effect, only the striker and the
incidence bar are needed. Table 3.5 shows dimensions and material properties of the individual
bars. The strain wave propagated into the incidence bar is altered due to the plastic deformation
of the shaper during impact. Frew [8] proposed the following stress-strain model for the copper
shaper while Fig. 3.40 shows the plot of it, i.e. Eqn. 3.53,
(3.53)
where

625	

0.32 and

,

Young’s modulus is 117	

4.25 . The unloading process is elastic and the

.

Table 3.5 Dimensions and material properties of the SHPB used
striker

shaper

incidence bar

length

152mm

1.6mm

2130mm

diameter

12.7mm

4.8mm

12.7mm

Young’s modulus

200 GPa

Fig. 3.40

200 GPa

Mass density

8100

/

8910
41

/

8100

/

Similar to those presented in Section 3.6, the stresses from Fig. 3.40 are divided by the strain to
gain an equivalent Young’s modulus

. The deformed cross-sectional area of the shaper, , is

calculated from the following equation
(3.54)
where

is the deformed shaper length,

length. The equivalent strain

is the original shaper area and

is the original shaper

is modified by the change of the cross-sectional area as shown

below
(3.55)

The bond stiffness involved in the shaper can then be calculated by using Eqn. 3.35. It should be
pointed out that the significant change of area should be considered in Eqn. 3.53 because the
deformation of the shaper can be quite large, deviating from one-dimensional assumption.

Simulation of the shaped incidence wave is compared with the experiment results in Fig. 3.41.
They agree with each other very well. Overall, this study shows that peridynamics is able to
simulate the shaping effect required for SHPB operations. It provides an efficient way to select a
shaper material with associated dimensions without numerous trial-and-errors.

3.8 Conclusions
A one-dimensional peridynamic model is proposed and its engineering applications are
evaluated. The numerical solutions from peridynamic analyses are well compared with the
associated analytical solutions in both convergence and accuracy studies. The feasibility of using
42

peridynamics in modeling Split Hopkinson Pressure Bar (SHPB) is also validated by the
experiments. Since SHPB tests performed in this dissertation are essentially focused on
validating the peridynamic theory, the tests are specifically designed for comparison purpose,
some assumptions and requirements of SHPB technique are not be completely satisfied. For
example, it is required in SHPB tests that the specimen diameter be smaller than that of the
striker, incidence and transmission bars. Study Three, however, does not follow this guideline.
The experiment data from Study Three is able to validate Eqn. 3.35 but may be not able to get
the accurate material properties from Eqn. 3.29 and Eqn. 3.30. Peridynamics is convenient for
simulating impact process with its long range force. All fundamental elements for mechanical
analysis are included in the formulation of peridynamics. No additional theory, such as contact
mechanics, is required. Because of its numerical nature, peridynamic method can be used to
simulate the shaper of the SHPB. This can help accelerate the use SHBP in dynamic material
characterizations significantly.

43

Figure 3.1 Calculation of stress.

44 
 

Figure 3.2 Wave propagation in a one-dimensional bar.

45 
 

Figure 3.3 Analytical solution for wave propagation.

46 
 

Figure 3.4 Number of nodes within a horizon.

47 
 

Figure 3.5 Comparison between classical mechanics and peridynamics.

48 
 

Figure 3.6 Convergence of time step size.

49 
 

Figure 3.7 Set up of Split Hopkinson pressure bar (SHPB).

50 
 

Figure 3.8 Wave propagation in SHPB.

51 
 

Figure 3.9 Square wave input for computational simulation.

52 
 

Figure 3.10 Comparison in Study One with a square wave input.

53 
 

Figure 3.11 Study Two: Two Bar.

54 
 

Figure 3.12 Study Two: comparison between experiment results and numerical solution, strain
history in incidence bar (top) and strain history in transmission bar (bottom).

55 
 

Figure 3.13 Study Two: comparison between experiment results and numerical solution with
grease, strain history in incidence bar (top) and strain history in transmission bar (bottom).

56 
 

Figure 3.14 A bond crossing an interface.

57 
 

Figure 3.15 Study Three: aluminum specimen with the same cross-section area as the split bars.

58 
 

Figure 3.16 Study Three: Aluminum specimen with the same-cross section area as the split bars,
incidence bar strain history (top) and transmission bar strain history(bottom).

59 
 

Figure 3.17 Study Four: Steel specimen with a smaller cross-section area.

60 
 

Figure 3.18 Study Four: Steel specimen with a smaller cross-section area, incidence bar strain
history (top) and transmission bar strain history (bottom).

61 
 

Figure 3.19 Study One: computational results with filter compared with experiment results.

62 
 

Figure 3.20 Study Two: Computational results with filter compared with experiment results.

63 
 

Figure 3.21 Study Three: Computational results with filter compared with experiment results.

64 
 

Figure 3.22 Study Four: Computational results with filter compared with experiment results.

65 
 

Figure 3.23 Incidence wave from Study Two.

66 
 

Figure 3.24 Study One: Incident wave as computational input.

67 
 

Figure 3.25 Study Two: Incident wave as computational input.

68 
 

Figure 3.26 Study Three: Incident wave as computational input.

69 
 

Figure 3.27 Study Four: Incident wave as computational input.

70 
 

Figure 3.28 Split bar system for validating impact process.

71 
 

Figure 3.29 Simulation of impact process.
For interpretation of the references to color in this and all other figures, the reader is referred to
the electronic version of this dissertation.

72 
 

Figure 3.30 Comparison between experiment and simulation with velocity input.

73 
 

Figure 3.31 Impact of two cylindrical bars right before impact.

74 
 

Figure 3.32 Impact of two cylindrical bars after impact.

75 
 

Figure 3.33 SHPB system [6].

76 
 

Figure 3.34 Stress-strain relation of the aluminum specimen from [6].

77 
 

Figure 3.35 Comparing computational results with experiment results from [6].

78 
 

Figure 3.36 Computational results with filter compared with experiment results.

79 
 

Figure 3.37 Incidence wave after filter from Ref. [8].

80 
 

Figure 3.38 SHPB with a shaper.

81 
 

Figure 3.39 Split bar system in Ref. [8].

82 
 

Figure 3.40 Stress-strain relation for the copper shaper given in Ref. [8].

83 
 

Figure 3.41 Simulation results compared with experiment results from Ref. [8].

84 
 

REFERENCES

85 
 

References
[1] Thomson, W. T., Dahleh, M. D., Theory of Vibration and Applications, 5th edition, Prentice
Hall (1997)
[2] Oppenheim, A. V., Willsky, A. S., Nawab, S. H., Signals and Systems, 2nd edition, Prentice
Hall (1996)
[3] Chen, W. W., Song, B., Split Hopkinson (Kolsky) Bar: Design, Testing and Applications
(Mechanical Engineering Series), 1st edition, Springer (2010)
[4] http://www.mathworks.com/help/techdoc/ref/fft.html
[5] http://www.mathworks.com/help/techdoc/ref/ifft.html
[6] Li, G., Construction, calibration and applications of a Split Hopkinson Pressure Bar, thesis
for the Degree of M. S., Michigan State University (2002)
[7] http://www.mathworks.com/products/curvefitting/
[8] Frew, D. J., Forrestal, M. J. and Chen, W., Pulse shaping techniques for testing brittle
materials with a split Hopkinson pressure bar, Experimental Mechanics, 42, 93-106 (2002)

86 
 

Chapter 4
Two-Dimensional Peridynamic Studies

The one-dimensional bar-like model presented in the previous chapter is excellent for onedimensional analysis. This chapter investigates a beam-like model for two-dimensional analysis.
The two-dimensional governing equation of peridynamics [1] is
																																																								

	

(4.1)

where the integration domain is horizon .

4.1 Bar model
Similar to the bar model given in section 3.1, a simple model of bond force can be expressed as
∙

(4.2)

where c is a material property, such as stiffness. To identify , it is necessary to find strain energy
density based on peridynamics and set it equal to that based on classical mechanics.

Consider a plane stress problem with the following displacement field

0

																																																																																

(4.3)

as shown in Fig. 4.1.

For the point at the origin , the bond stretch is
(4.4)

87

From Eqn. 2.7, it is defined that
																																													

(4.5)

Substituting Eqn. 4.5 into Eqn. 4.2, it yields
∙

∙ /

(4.6)

Based on the definition of the micro-potential (Eqn. 2.9), it follows that
/2

(4.7)

This is the strain energy in one bond. Integrating

over the horizon of , the strain energy

density becomes
	

	

	

(4.8)

The strain energy density can also be calculated from theory of elasticity. From Eqn. 4.3, the
strain components [2] are
							

(4.9)
	

																																																																										
																																																																											

(4.10)

0	

(4.11)

Using Hook’s law [2], the stress elements are found to be
	

(4.12)
(4.13)

88

The strain energy density from classical mechanics then becomes
(4.14)

Setting the two strain energy densities, i.e. Eqn. 4.8 and Eqn. 4.14, equal to each other, it yields
(4.15)

Rearranging Eqn. 4.15, the material property can be found as follows
(4.16)

4.2 Definition of stress in peridynamics
Similar to Section 3.2, the definition of stress is only sought for comparison purposes since there
is no need of stress in peridynamic simulations. In two-dimensional analysis, the traction [3] at
point

is defined as the total internal force applied on the segment [

right-hand side of

, 0] by the domain on the

as shwon in Fig. 4.2 and can be expressed by the following equation
	

where

(4.17)

is the domain on the right-hand side of point , ′ is a point in

integration element area of ′ and
all points on the

is the bond force applied on

apply forces not only on point

, which is the horizon. The

componet of

																																																											

is
	

,

is the

by ′. It should be noted that

but also beyond

along a layer of thickness

and is defined as
(4.18)

89

Consider a deformation with strains

and
cos

Then the

, the bond stretch should become
sin

(4.19)

component of the bond force is

																																																													

∙ ∙ cos

(4.20)

Substituing Eqn. 4.20 into Eqn. 4.18, it arrives
	 	
																		

	 ∙

∙

	

	

(4.21)

The integration domain shown in Fig. 4.3 can be explained as that 	
force applied by
because

∙

on point A. The distance from point

is the

to point A is

component of the

which covers from 0 to

has no interaction with points beyond . For point A, which is located a distance

from , the distance between a point in
largest value of

and point A is . The smallest value of

is . For every , only points from

domain. Points outside this range are not in

cos

	 to

cos

is

and the

are in

domain. The symbolic integration can be

calculated by Mathematica [4].

In classical mechanics,

is defined as follows
(4.22)

90

1/3. There is only one material

Eqn. 4.21 and Eqn. 4.22 are equal to each other when

property defined in the bar model based on peridynamic analysis as opposed to two independent
properties used in classical mechanics. Hence, the Poisson’s ratio for the peridynamic bar model
has to be 1/3.

4.3 Beam model
In order to accommodate two material properties, an upgraded peridynamic model resembling a
beam is proposed here. The interaction force between two points depends not only on the
deformation along the axial direction but also along the transverse direction as depicted in Fig.
4.4.

Consider a point

located away from

based on the local coordinate on

, with an angle . The bond forces between

and

,

, are
′

and r is the original length of the bond

′ /

(4.23)

′
where ′ is the axial force appiled on

′
′

′ /

(4.24)

by

,

′ is the transverse force applied on

by

- . The transformation equations between the local

coordinate and the global coordinate are
	

cos
sin

Consider a case with strains in

and

directions are

Eqn. 4.23 and Enq. 4.24 become
91

sin

(4.25)

cos

and

(4.26)

, i.e.

		and

.

′
′

1

2 cos

2 sin

2

2 sin

2 cos

1 cos

3

1 sin

1 sin

	

2

1 cos

1

cos

1

1

cos sin

2

sin
2

2

/

(4.27)

cos sin

/

3

(4.28)

The strain energy of this bond becomes
(4.29)

Integrating Eqn. 4.29, the strain energy density at

becomes

	
3
where

3

2

3

(4.30)

is the horizon.

On the other hand, from classical mechanics, the stresses are defined as
(4.31)
(4.32)

The strain enery density is
(4.33)

Set Eqn. 4.30 and Eqn. 4.33 equal to each other, the following solutions can be obtained

92

(4.34)

(4.35)
where

and

are the two independent material properties of the beam model. When

and

1/3

0, the beam model can be reduced to the bar model presented in Section 4.1.

From Eqn. 4.34 and Eqn. 4.35, Young’s modulus
by

and Poisson’s ratio

can also be expressed

and
(4.36)

1

From Eqn. 4.37, it can be concluded that

(4.37)

/

can vary from -1 to infinity. This implies that the

peridynamic beam model is capable of modeling any isotropic material which has

raging from

-1 to 0.5.

Similar to Section 4.2, the stresses can be calculated from the peridynamic beam model and
compared with those given in classical mechancis. Consider a deformation with strains
and

. In the beam model, the

component of bond force is
′ ∙ cos

′ ∙ sin

This is different from the bar model which gives
expressed as
93

(4.38)

′ ∙ cos . Eqn. 4.18 can then be

	

3

∙ cos

∙ sin

3

(4.39)

Substituing Eqn. 4.34 and Eqn. 4.35 into Eqn. 3.39, it yields
(4.40)

The interpretation of Eqn. 4.39 is similar to that of Eqn. 4.21. As can be seen, Eqn. 4.40 is
exactly the same as

in Eqn. 4.22 obtained from classical mechanics analysis.

4.4 Computational efficiency
From Fig. 2.4, it can be seen that there are three loops in peridynamic calculations, namely time
loop, domain loop and neighbor loop. With the loops, small time step and small distance between
nodes, peridynamic calculations will become very time consuming. For example, Study Three in
Chapter 3 has about 3500 nodes in the one-dimensional domain. It takes about two hours to
complete the calculation using the serial computational algorithm shown in Fig. 2.4. Hence, it
can be expected that two-dimensional peridynamic problems containing more nodes will further
slow down the calculation process if the serial computational algorithm shown in Fig. 2.4 is
used. This section studies methods to improve the compuational efficiency of peridyamics.

4.4.1 Using matrix computation
For problems involving infinitesimal deformation, the relationships between nodes may be
94

considered unchanged during the deformation process. The calculation of these relationships in
each time step shown in the flowchart given in Fig. 2.4 may be disregarded. For computational
efficiency, it is desired to calculate these relationships before the beginning of computation and
store them in a matrix like the stiffness matrix practiced in finite element analysis. Therefore,
only matrix multiplication is required in each time step.

For two-dimensional analysis, the peridynamic equation of motion can be rewritten in the
following discrete form based on node	
∑

	
where

	∆

(4.41)

is any node in the horizon of point shown in Fig. 4.5. Other notations are mass density

, acceleration (of node )

, force (between node and node	

, area (of node

body force (of node

. The following notations are also defined below:

nodal positions: [ ,

for node and [

node displacements: [

,

∆ and

,

] for node
,

] for node and [

] for node

,

relative displacement:

	

,

	

force between node and node	 :

∙ ∙

relative position:
bond stretch:

The

|

|
| |

component of

1

|

(4.42)

|

is
,

,

,

∙ ∙

95

|

|

(4.43)

To express
,

,

in a matrix form, the following truncated Taylor’s series may be used
,

0,0,0,0

0,0,0,0 ∙
2
2

0,0,0,0 ∙

0,0,0,0 ∙

component of

,

,

,

0,0,0,0 ∙

0,0,0,0 ∙

2

2

0,0,0,0 ∙

0,0,0,0 ∙
0,0,0,0 ∙

2

0,0,0,0 ∙

0,0,0,0 ∙

⋯

(4.44)

is
∙ ∙

Similarly, to express

0,0,0,0 ∙

2

0,0,0,0 ∙
0,0,0,0 ∙

The

0,0,0,0 ∙

|

(4.45)

|

in a matrix form, the following truncated Taylor’s series may be used
0,0,0,0

0,0,0,0 ∙

0,0,0,0 ∙

0,0,0,0 ∙

0,0,0,0 ∙
0,0,0,0 ∙
2

0,0,0,0 ∙

2

2

0,0,0,0 ∙
0,0,0,0 ∙

0,0,0,0 ∙

2

0,0,0,0 ∙

2

0,0,0,0 ∙

0,0,0,0 ∙
0,0,0,0 ∙

2

⋯

0,0,0,0 ∙
(4.46)

The equation of motion can then be written in a matrix form

96

∙
where

∙

∙

∙

(4.47)

is a displacement vector and is defined as

(4.48)

is the first-order term matrix (two-dimensional symmetric matrix) and is defined as
0,0,0,0

,

0,0,0,0

,

0,0,0,0

,

0,0,0,0

,

0,0,0,0

,

0,0,0,0

,

0,0,0,0

,

0,0,0,0

,

M is the second-order term matrix (three-dimensional matrix) and is defined as

97

(4.49)

,

1
2

,

0,0,0,0

, ,

,

0,0,0,0

,

,

0,0,0,0

,

1
2

, ,

,

0,0,0,0

0,0,0,0

,

,

0,0,0,0

,

,

1
2

,

,

,

0,0,0,0

0,0,0,0

,

,

1
2

,

,

1
2

,

0,0,0,0

0,0,0,0

0,0,0,0

, ,

0,0,0,0

,

,

0,0,0,0

0,0,0,0

,

98

1
2

, ,

,

0,0,0,0

0,0,0,0

,

,

0,0,0,0

,

,

,

,

,

1
2

0,0,0,0

0,0,0,0

,

0,0,0,0

,

(4.50)

Eqn. 4.49 and Eqn. 4.50 are calculated before the first time step. Once the K matrix and the M
matrix are built, each step can be carried out by Eqn. 4.47. The computational flowchart is shown
in Fig. 4.6. As compared with that shown in Fig. 2.4, two layers of loop in each time step are
eliminated. A significant amount of time can then be saved. For a one-dimensional problem with
about 3500 nodes as used in Study Three of Chapter 3, it takes about two hours to complete the
simulation on a PC with an Intel(R) Core(TM) i7-2860QM CPU based on the algorithm showing
in Fig. 2.4. However, it only requires approximately five seconds to complete the simulation on
the same PC using the matrix computation algorithm.

4.4.2 Using parallel computing
Another problem continues to hinder the computational process in the matrix computation is
associated with damage process. When damage takes place, bond stretches still have to be
calculated and checked one by one to identify the locations and directions of the damage. Once
99

damage is identified, the corresponding elements of

and M will have to be nulled to represent

the loss of bond stiffness due to the broken bonds. This process will also slow down the
peridynamic computation and hinders the convenience of the peridynamic method in damage
problems.

Since the computation of peridynamic nodes is independent of one another, there is a great
advantage to use parallel computing in peridynamic simulations. Therefore, the algorithm shown
in Fig. 2.4 is modified for parallel computing. As shown in Fig. 4.7, the computation of
individual nodes are carried out independently.

In this study, the Matlab parallel computing toolbox [5] is used to implement the parallel
computing. . The program is modified as shown in Fig. 4.8. As can be seen from Fig. 4.7, the
neighboring loop is also eliminated by using matrix computation in Matlab. Instead of processing
one neighboring node to another, it is possible to use Matlab matrix computation to manipulate
all neighboring nodes at the same time. The facility located at High Performance Computing
Center (HPCC) at MSU [6] is used to run the modified program. However, due to the limited
resources at the HPCC, eight workers are considered as a balanced choice between queue time
and computation time. When solving a two-dimensional problem with 10,201 nodes, it takes 215
seconds to compute 100 time steps based on serial computing while it only takes 27 seconds to
complete the same task based on eight-worker parallel computing. Ideally, if resource is
available, it is possible to use the same number of workers as the numbber of nodes to
significantly reduce the computational time.

100

4.5 Comparison with vibration theory
4.5.1 Vibration theory
, , , consider the free vibration of a circular

In a two-dimensional polar coordinate system

plate, shown in Fig. 4.9, under the following initial conditions
	

, ,0

(4.51)

, ,0

Since the circumference of the circle,

0

(4.52)

, is a free boundary, the boundary condition should

be
, ,

0

(4.53)

Apparently, this is an axisymmetric problem and the displacements can be redefined as
, ,

,

, ,

0

(4.54)
(4.55)

Hence, the strains are [2]
(4.56)
(4.57)
0

(4.58)

By substituting Eqn. 4.56, Eqn. 4.57 and Eqn. 4.58 into Hook’s law, the following stresses can be
found

101

(4.59)

	

(4.60)

0

(4.61)

The polar coordinates based equation of motion in r-direction is defined as
(4.62)

The equation of motion can be solved using numerical method. With
2800	

/

and

75	

,

0.3,

0.1	 , the strain history of point B on the circular plate is shown in

Fig. 4.10.

4.5.2 Numerical studies
There are three parameters in numerical calculation. They are time step
horizon

, number of nodes in

and the size of horizon . As discussed in Chapter 3, the smaller the horizon, the

closer the solution will be to that of classical mechanics. When the horizon is fixed, more nodes
within the horizon, i.e. smaller node size, will give higher solution accuracy. However, smaller
node size and higher node number in the horizon will require higher computational time. A
practical computation should have a balance between computational efficiency and numerical
accuracy.

Twelve pairs of

and

are investigated for comparison. They are shown in Table 4.1. Fig. 4.11

shows the computational results of case 10, case 11 and case 12 along with their analytical
102

10

solution. Although

is a large horizon, it still can be seen that the increase of the node

number in the horizon, i.e. the decrease of node size, improves the accuracy of the peridynamic
solution to the theoretical solution. Fig. 4.12 compares case 5 and case 6 with the solution from
vibration theory. When the number of nodes in a horizon is fixed at 4, two different horizons,
2.2	

2.5	

and

, are used. As shown in Fig. 4.12, the smaller the horizon, the

closer the peridynamic solution becomes to the solution from vibration theory solution.

Table 4.1 Twelve pairs of

and

1
1

3

4

5

1

10

2

1.6	

3

4

2.2	

5

2.5	

6

5	

7

8

9

10	

10

11

12

When

4 and

2.2	

, the size of time step

1

different time steps,

10 s and

1

is studied. As shown in Fig. 4.13, two

10 s, are used and the results are compared

with the results from the vibration theory. The result with
that from

1

10

1

10

10

is slightly better than

at the rising time of the radial strain. Other than that, the two results
1

are almost the same. However, the case with
case with

1

.
103

10

is about ten times slower than the

0.5

Based on these results,

,

2.2	

and

1

10

are used and the strain

histories of point A, point B and point C, as shown in Fig. 4.9, are compared with the vibration
theory results in Fig. 4.14, Fig. 4.15 and Fig. 4.16 respectively. They agree well.

4.6 Convergence to the one-dimensional results
To further verify the two-dimensional peridynamic model, it is also of interest to compare the
results from the two-dimensional peridynamic simulations with those from the one-dimensional
peridynamic simulations.

The free vibration of a rectangular plate as shown in Fig. 4.17 is considered. The length of the
0.4	

plate is
0.3 and

/ . The material properties are

and the aspect ratio is
2800	

/

75	

,

. The boundaries are all free and the initial conditions are
,

∙

,

(4.63)
∙

(4.64)

The vibration process is simulated by the two-dimensional peridynamic method and the strain
history of point A, shown in Fig. 4.17, is recorded and compared with the results from the onedimensional peridynamic analysis.

Figs. 4.18, 4.19, 4.20, 4.21 and 4.22 show the comparison for
and

0.2,

1,

10,

50

100, respectively. As can be seen, the results based on the two-dimensional theory

converge to those based on the one-dimensional theory when the aspect ratio becomes larger.

104

4.7 Axisymmetric problems
In this section, an axisymmetric model of peridynamics is proposed. Similar to classical
mechanics, the calculation of an axisymmetric problem can be simplified and the dimension can
be reduced. A two-dimensional axisymmetric problem can be reduced to a one-dimensional
problem. This will simplify the calculation and reduce the computational time significantly.

For a two-dimensional axisymmetric problem, as shown in Fig. 4.23, consider a point
The governing equation of point

using local radial coordinate

,

	
where

The

is horizon and

.

is expressed as follows
(4.65)

is the force applied on point P by point Q in global radial direction .

component of bond force

is
∙ ∙ sin

(4.66)

and is the stretch between the two points, which can be defined as
|

|
|

1

|

(4.67)

where
(4.68)
cos

sin

(4.69)
(4.70)

cos

sin

105

(4.71)

Substituting Eqn. 4.67, Eqn. 4.68, Eqn. 4.69 and Eqn. 4.70 into Eqn. 4.71, it yields
	

	

1

The two radial coordinate systems are used for convenience. They are local coordinate
global coordinate

,

(4.72)

,

and

. The relationship between them, as shown in Fig. 4.22, is
2

(4.73)

tan

	

(4.74)

After substituting Eqn. 4.73 and Eqn. 4.74 into Eqn. 4.72, the governing equation at point P in
global coordinate can be found as
,

	

(4.75)

where the integration bound can be found from the triangle shown in Fig. 4.24.

The innermost integration of Eqn. 4.75 is a function of

only. Let
,

(4.76)

Eqn. 4.75 then becomes
	

(4.77)

A two-dimensional problem is therefore reduced to a one-dimensional problem with the
106

assumption of axisymmetry.

The same problem solved in Section 4.5.1 is simulated using the axisymmetric peridynamic
model proposed above and the strains at points A, B and C are compared with the results from
Section 4.5.1 in Fig. 4.25, Fig. 4.26 and Fig. 4.27. Good matches are shown between the two
results.

4.8 Failure theory
The critical stretch

in Eqn. 2.12 can be found from the energy method similar to Eqn. 4.7.

Consider a fracture surface in a large homogeneous body. In order to completely seperate the
body into two halves, it will require breaking all the bonds crossing the two halves. Let
denote the work required to break a single bond,
	
where

(4.78)

.

Using Eqn. 4.2, Eqn. 4.78 becomes
																																																											

The work

(4.79)

required to break all the bonds per unit fracture surface area shown in Fig. 4.28 can

then be expected as follows
																																					

/
/

∙

The calculation of this intrgration is similar to that of Eqn. 4.21.

107

(4.80)

Solving for

, the critical bond stretch can be related to the energy per unit fracture area for

completely separating the two halves of the body, i. e.
(4.81)

4.9 Comparisons of crack propagation velocity with experiments
It is challenging to validate the dynamic crack propagation obtained from numerical simulation
since there is no analytical solution which can be readily compred with the numerical solution
[7]. Therefore, experiment results are often used to compare and validate the numerical models.

Many crack propagation problems were simulated by finite element method. Song [7] compared
three different finite element analysis (FEA) techniques with experiment results. The three
methods are the extended finite element method (XFEM), element deletion method and
interelement crack method. The element deletion method is unable to predict crack branching.
The XFEM and the interelement method showed similar crack velocity and crack paths but both
failed to predict a benchmark experiment without an adjustment of the energy release rate.

Crack tip propagation from peridynamic analysis are compared with exerperimental study here.
Boudet [8] studied the crack propagation in a PMMA plate which has a mass density
1200	
100

/

and a Young’s modulus
290

plate with a 10

2.5	

. The experiments were conducted on a

notch in the middle of a boundary, as shown in Fig. 4.29.

The velocity of the crack propagation was measured using several equally spaced conductive
strips deposited on one of the surfaces of the sample.

108

In the experiment, the specimen was first loaded by displacement control and the apllied
displacement was slowly increased every 10 seconds. Because 10 seconds is too long for
computaitonal method, an initial strain just below what is needed to trigger crack propagation is
applied as the initial condition. The starter crack is made by releasing all bonds crossing the
10

notch.

Fig. 4.30 shows the predicted crack length as a function of time. Using the predicted crack
length, crack propagation velocity can be calculated and is shown in Fig. 4.31. The maximum
crack propagation velocity is about 600 m/s, which is less than the Rayleigh wave speed of the
PMMA material, 930 m/s. The predicted crack propagation velocity is within the theoretical limit
for steady mode-I fracture [9]. The peridynamic results are compared with the experiment results
given in Ref. [8]. Both Fig. 4.30 and Fig. 4.31 show good agreements between the two results.

4.10 Dynamic crack branching
Experiments on crack branching were reported in literature [10-14]. In these experiments, a
crack starts to propagate from a notch. At a certain point, the crack branches into at least two
cracks. The angle of the branches, however, varies from one experiment to another. Sharon and
Fineberg [13] showed that the branching angles were in Gaussian distribution with an average of
30° . Numerical methods to predict crack branching were studied. The Yoffe calculation [15]
predicted an angle of 60° while molecular simulations [16] predicted an angle of 30° . An angle
of 18° was also predicted using an energy criterion with consideration of nonsingular terms in
the stress field [17].

109

The qualitative comparison between peridynamic prediction of crack branching and associated
experiment study is to be examined. The goal is to see if peridynamics is able to predict crack
branching without requiring additional theories.

Consider a pre-notched glass plate with

2450	

dimensions of the plate is 100	

with a 10	

	

100	

left boundary as shown in Fig. 4.32. A tensile stress

/

,

32	

and

0.2 . The

notch located in the middle of the
2	

is applied at both top and

bottom boundaries. The initial displacement and the initial velocity is zero. Fig. 4. 33(a) shows
the crack path at different times calculated by peridynamics. Fig. 4.33(b) is strain energy density
of the corresponding steps. Crack branching can be clearly seen in Fig. 4.33(a) and Fig. 4.33(b)
and there is always stress concentration at the crack tips.

4.11 Edge-cracked plate under impulsive loading
Kalthoff and Winkler [18] studied crack propagation of an edge-cracked plate under impulsive
loading and found that the crack propagated in a direction approximately 70° from the origianl
crack. In this section, the same problem is simulated by the proposed peridynamic model and the
results are to be compared with the experiment results from [18].

Consider a 100

200

maraging steel 18Ni1900 plate as shown in Fig. 4.34. There are

two parallel notches along the left edge of the plate. The material properties [19] are
8000	

/

,

190	

and

0.3. The two edge notches are impacted by a projectile

with an intial velocity of 36	 / . This impact results in a compressive wave propagating to the
interior of the plate. Once the compressive wave arrives at the notch tips, a mode-II brittle failure

110

occurs. The crack do not propagate in the same direction as the original crack but in a direction
that is about 70° away from the original direction.

In numerical simulation based on peridynamics, it is assumed that the projectile has the same
material properties as the specimen. Hence, one half of the initial velocity, 18	 / , is used as the
initial condition at the left edge [20, 21] and all boundaries are free. Fig 4.35 shows the simulated
crack path from peridynamic computation. The crack path is about 68° from the original notch.
Fig. 4.36 shows the horizontal displacement of the plate at

33.2	

. The displacement

discontinuity shows the crack path which is 68° from the original notch. The peridynamic
computation successfully predicts the experimentlly observed crack propagation angle.

4.12 Conclusion
A novel two-parameter beam model for peridynamic analysis is proposed in this chapter. The
numerical solution of peridynamics has been verified by an analytical solution and it converges
to the analytical solution. To improve the computational efficiency in two-dimensional
simulations, two computational techniques are implimented in this work. One uses matrix
computation and the other parallel computation. The first method works well for elastic problems
without damage while the second method is proved to be excellent for simulating dynamic
damage propagation. Failure theory of peridynamics is also studied and applied to simulate
several dynamic damage propagation problems. Three experiment studies from literature [8,
10,18] are used to validate the two-dimensional peridynamic model and algorithm.

Being a reformulation of continuum mechanics, peridynamics covers fracture mechanics

111

automatically. This poses a great potential for peridynamics since it is convenient for simulating
dynamic damage propagation. Compared with the commonly used finite element method, there
is no need of remeshing in peridynamics since peridynamics is mesh free. Peridynamics does not
require additional external theories for crack growth since it is controlled by bond strength.
Besides, peridynamics does not require tracking individual cracks since cracks occur when the
bonds are damaged. As shown in this work, peridynamics is also capable of simulating elastic
deformations without damage.

112

Figure 4.1 Two-dimensional domain under radial deformation.

113 
 

Figure 4.2 Definition of stress in peridynamics.

114 
 

Figure 4.3 Calculation of stress at point .

115 
 

Figure 4.4 Bond force in beam model.

116 
 

Figure 4.5 Original configuration and current configuration.

117 
 

Figure 4.6 Flow chart of algorithm in peridynamic matrix calculation.

118 
 

Figure 4.7 Flow chart for parallel computing.

119 
 

Original program: 
 
… 
for k=1:nt                 % time loop  
 
for i=1:n       % domain 
loop 
 
 
 
… 
 
 
end 
end 
… 

Current program: 
 
… 
matlabpool(8)       %open matlab pool: 8 workers 
for k=1:nt                      % time loop 
 
parfor i=1:n     % domain loop 
 
 
 
… 
 
 
end 
end 
… 

Figure 4.8 Code before parallel modification and code after parallel modification.

120 
 

Figure 4.9 Vibration of a circular plate.

121 
 

Figure 4.10 Strain history on a circular plate.

122 
 

Figure 4.11 Convergence of numerical solution to analytical solution.

123 
 

Figure 4.12 Convergence of numerical solution to analytical solution.

124 
 

Figure 4.13 Convergence of numerical solution to vibration theory.

125 
 

Figure 4.14 Peridynamic calculation of strains at point A compared with results from vibration
theory.

126 
 

Figure 4.15 Peridynamic calculation of strain at point B compared with results from vibration
theory.

127 
 

Figure 4.16 Peridynamic calculation of strain at point C compared with results from vibration
theory.

128 
 

Figure 4.17 Rectangular plate for free vibration study

129 
 

Figure 4.18 One-dimensional result compared with two-dimensional result for aspect ratio of 0.2.

130 
 

Figure 4.19 One dimensional result compared with two-dimensional result for aspect ratio of 1.

131 
 

Figure 4.20 One-dimensional result compared with two-dimensional result for aspect ratio of 10.

132 
 

Figure 4. 21 One-dimensional result compared with two-dimensional result for aspect ratio of 50.

133 
 

Figure 4. 22 One-dimensional result compared with two-dimensional result for aspect ratio of
100.

134 
 

Figure 4.23 A two-dimensional axisymmetric model.

135 
 

Figure 4.24 Angular integration boundary of Eqn. 4.75.

136 
 

Figure 4.25 Comparison of peridynamic calculation of strain at Point A using axisymmetric
model with results from the vibration theory.

137 
 

Figure 4.26 Comparison of peridynamic calculation of strain at Point B using axisymmetric
model with results from the vibration theory.

138 
 

Figure 4.27 Comparison of peridynamic calculation of strain at Point C using axisymmetric
model with results from the vibration theory.

139 
 

Figure 4.28 Calculation of the critical stretch.

140 
 

Figure 4.29 Schematic drawing of the sample for crack propagation test.

141 
 

Figure 4.30 Comparison of crack propagation simulated by peridynamics with experiment result.

142 
 

Figure 4.31 Comparison of crack propagation speed simulated by peridynamics with experiment
result.

143 
 

Figure 4.32 A plate with single notch for crack branching study.

144 
 

(a) damage

(b) strain energy density

200

214

224

257

289

Figure 4.33 Computational crack path (a) damage and (b) strain energy density simulated by
peridynamics.

145 
 

Figure 4.34 A two-notch plate under impulsive loading.

146 
 

Figure 4.35 Computational crack path of the top notch.

147 
 

Figure 4.36 Horizontal displacement at 33.2	

from computational study.

148 
 

REFERENCES

149 
 

References
[1] Silling, S. A., Reformulation of elasticity theory for discontinuities and long-range forces, J.
Mech. Phys. Solids, 48, 175-209 (2000)
[2] Timoshenko, S. P. Goodier, J. N., Theory of Elasticity, Third edition, McGraw-Hill
Publishing Company (1970)
[3] Silling, S.A., Askari, E., A meshfree method based on the peridynamic model of solid
mechanics. Comput. Struct, 83, 1526–1535 (2005)
[4] http://www.wolfram.com/mathematica/
[5] http://www.mathworks.com/products/parallel-computing/
[6] http://www.icer.msu.edu/?q=hpcc
[7] Song, J. H., Wang, H., Belytschko, T., A comparative study on finite element methods for
dynamic fracture, Computational Mechanics, 42, 239-250 (2008)
[8] Boudet, J. F., Ciliberto, S., Steinberg, V., Dynamics of Crack Propagation in Brittle Materials,
J. Phys., 6, 1493-1516 (1996)
[9] Yoffe, E. H., The moving Griffith crack, Philosophical Magazine, 42, 739-750 (1951)
[10] Ramulu, M., Kobayashi A. S., Mechanics of crack curving and branching-a dynamic
fracture analysis, Int. J. Fract., 27, 187-201 (1985)
[11] Ravi-Chandar, K., Dynamic fracture of nominally brittle materials, Int. J. Fract., 90, 83-102
(1998)
[12] Sharon, E., Gross, S. P., Fineberg J., Local crack branching as a mechanism for instability in
dynamic fracture, Phys. Rev. Lett., 74, 5096-5099 (1995)
[13] Sharon, E., Fineberg, J., Microbranching instability and the dynamic fracture of brittle
materials, Phys. Rev. B, 54, 7128-7139 (1996)
[14] Fineberg, J., Sharon, E., Cohen, G., Crack front waves in dynamic fracture, Int. J. Fract.,
119, 247-261 (2003)
[15] Marder, M., Liu, X., Instability in lattice fracture, Phys. Rev. Lett., 71, 2417 (1993)
[16] Abraham, F. F., Brodbeck, D., Rafey, R. A., Rudge, W. E., Instability dynamics of fracture:
A computer simulation investigation, Phys. Rev. Lett., 73, 272 (1994)

150 
 

[17] Michopoulos, J. G., Theocaris, P. S., The Directional Instability of the Running Crack: A
Deterministic and Stochastic Approach, Int. J. Engng. Sci., 29, 13-36 (1991)
[18] Kalthoff, J. F., Winkler, S., Failure mode transition at high rates of shear loading,
International conference on impact loading and dynamic behavior of materials, 1, 185-195 (1987)
[19] Decker, R.F., Source book on maraging steels, American Society for Metals (1979)
[20] Lee, Y.J., Freund, L.B., Fracture initiation due to asymmetric impact loading of an edge
cracked plate, J. Appl. Mech, 57, 104-111 (1990)
[21] Kalthoff, J.F., Modes of dynamic shear failure in solids, Int. J., Fract., 101, 1-31 (2000)

151 
 

Chaper 5
Orthotropic Model Analysis

It has been shown in Chapter 4 that peridynamic simulation of damage process does not require
any knowledge of the damage location and orientation prior to the simulation. This is
fundamentally different from finite element analysis which requires knowledge of damage
location and orientation in advance to impose special finite element mesh, such as initial damage
elements and cohesive zone layers [1], for damage simulations. This prerequisite becomes even
more challenging when inhomogeneous and anisotropic composite materials are of interest. In
addition, peridynamic simulation does not require remeshing at the end of each damage
processing step since it is a mesh free method. On the contrary, finite element analysis does.
Based on these difference, peridynamics should be more suitable for simulating dynamic damage
process in composite materials which have different properties in different locations and different
orientations.

Quite some simulations of composite damage process have been available in the literatures.
Dwivedi [1] modeled the propagation of single-edge notch (SEN) in 0° laminated plate using
cohesizve zone method. Xu [2] and Hu [3] proposed a two-parameter discrete peridynamic
model for composite damage simulations, in which there were two kinds of bonds: fiber bond
and matrix bond. Two material properties,

and

, were associated with the two types of

bonds. Only the bonds along the fiber direction were associated with the material property
while all other bonds with the material property

. This model required remeshing for different

fiber directions. For example, a 0°-90° grid mesh could only be used for a 0° or 90° laminae. For
152

a 45° lamina, a grid mesh consisting of 45° and 135°	was required. The two-parameter model
was an aproximation of the four material properties invovled in orthotropic materials. They were
mainly associated with two Young’s moduli,

and

. Its capability of modeling shear behavior

is unknown.

In this study, a continuous orthrotropic material model is proposed. It is based on continuous
trigonometric functions. With the continuous material property functions, it is not necessary to
have bond in fiber direction and therefore, this model is mesh independent.

5.1 Bar model for orthotropic materials
This model is based on the bar model presented in Seciton 4.1. The peridynamic equation of
motion [4] in two-dimensional domain can be expressed as
																																																								
where

	

	

	

(5.1)

is external force. The force boundary condition can be included in the external force.

For bar model, the bond function

is
∙

(5.2)

where is bond stretch. Contrary to the isotropic material model, bond material property

is

assumed to be a trigonometric function
cos
where

,

and

bond direction and

cos

(5.3)

are constants and can be identified from composite material properties.
is the fiber direction as shown in Fig. 5.1.

153

is

Similar to the analysis in the previous chapters,

,

and

can be identified from comparing

the strain energy densities based on peridynamic analysis and those based on classical
mechanics.

Consider a composite plate with the fibers oriented in ° direction and subjected to the following
strain field
(5.4)
(5.5)
(5.6)
The three components are independent of one another.

For a bond in

direction, and connected to a point
∙

cos

in the domain, the bond force should be

sin

sin cos

(5.7)

From Eqn. 4.7, the strain energy in the bond becomes
/2

(5.8)

Substituting Eqn. 5.3 and Eqn. 5.7 into Eqn. 5.8 and integrating
energy density at the point

should be

	
cos 2
8

over the horizon, the strain

	

	
cos 4

sin 2

2 3

16
4

sin 4

154

8
/768

3

2

3
(5.9)

Eqn. 5.9 can be simplified to find the coefficient of each independent term, as shown in Table
5.1. The simplification is achieved based on Mathematica [5].

Table 5.1 Simplified Eqn. 5.9 in terms of independent terms
Coefficients

Independent terms
/96

cos

/96

cos

/48
/96
3

cos
cos

4

5

/96

4

cos

/96

cos

/48

cos

/96

cos

3

8

48

/768

35

40

5

8

16

/384

5

8

16

/768

48

/768

3

4

/96

sin cos

5

4

/96

sin cos

/48

sin cos

/48

sin cos

155

On the other hand, from the theory of composite materials [6], the strain energy density under
,

and

is
(5.10)

The stresses in Eqn. 5.10 can be calculated by

(5.11)

where,
2

4

(5.12)

2

4

(5.13)
4

(5.14)

2

(5.15)

2

(5.16)

2

(5.17)

cos

(5.18)

sin

(5.19)

(5.20)

(5.21)

(5.22)

156

(5.23)
,

,

and

are the four material properties of orthotropic materials. Substituting Eqns.

5.11 - 5.23 into Eqn. 5.10, it yields
2
2

2
cos

1
sin

cos

2

2

2

2

2

2

2

cos

2

2

sin

1

2

2

sin

2
2

2
cos sin
sin 2
(5.24)

After simplification processes, the coeffiencet of each independent terms can be found. They are
listed in Table 5.2.

Table 5.2 Simplified Eqn. 5.12 in terms of independent terms
Coefficients

Independent
terms

4

4

2

cos

4

2

cos

2
4
2
4

4

cos

2

157

Table 5.2 (cont'd)

4

4

2

cos

2
2

1

2

2

2

cos

cos

4

4

2

cos

4

4

2

cos

2

2

2

2
2

1

2

4

4

sin cos

2

sin cos

2

4

4

2

2

sin cos

sin cos

158

Eqn. 5.9 should be equal to Eqn. 5.24. Hence, it is possible to set the coefficients of identical
terms equal to each other and identify

,

and

.

12

4

⁄

(5.25)

3

5

48

5

⁄

4

(5.26)

5

20

/

5

4
(5.27)

(5.28)

The composite model has three parameters and it can model shearing deformation. The bond
stiffness is a continuous function so the stiffness at any direction can be calculated. It can be
found that the bond becomes stronger as the angle between the bond and the fiber becomes
smaller. Similar to the bar model in Section 4.1, the Poisson’s ratio of this model is fixed and is
related to the remaining three independent material properties. For example, the Poisson’s ratio
of Kevlar/Epoxy is 0.34 while the Poisson’s ratio calculated from this model is 0.39.

159

5.2 Beam model for orthotropic materials
To accomondate the four material properties of orthotropic materials, a beam model based on
Section 4.2 is proposed here.

,

From the composite theory, there are four independent material properties,

,

and

.

The composite stiffness varies with fiber orientation. Bond functions are
(5.29)
(5.30)

They are identifcal to Eqn. 4.23 and Eqn. 4.24. However,

and

for a bond in

direction will

be dependent on fiber orientations as follows,
cos

cos

(5.31)
(5.32)

where

is the orientation of the fiber, as shown in Fig. 5.3. The coefficients

,

,

and

are

four independent material properties. They are related to the four material properties defined in
the composite theory.

Consider a composite plate with fiber in

direction and subject to the following strain field
(5.33)
(5.34)
(5.35)

The three strains are independent from each other.
160

Similar to Eqn. 4.29, the strain energy in one bond becomes
(5.36)

Integrate Eqn. 5.36 to find the strain energy density for a point
	

16

2
24
24

γ
6

γ

3

2 24

8
γ

γ

cos 4

4

9

16
γ

cos 2

8

9
γ

8

sin 4

12
12

24
γ

sin 2
(5.37)

The Strain energy density based on the composite theory is the same as Eqn. 5.24. Similar to
Section 5.1, by setting Eqn. 5.24 equal to Eqn. 5.37 and comparing the coefficients of the
independent terms, the following equations are obtained
(5.38)

(5.39)

(5.40)

																																				

(5.41)

5.3 Calculation of stresses from peridynamics
161

Similar to Section 4.2, stresses can be defined and calculated from peridyamics. They can be
used for some special comparison but not necessary in peridynamic simulations.

Consider a case with

and

. The stress

can be calcualted by using Eqn. 4.39

and is given below
	
2 3

4

8

∙ cos

∙ sin

3

16

48
cos 2

cos 4

(5.42)

Substituting Eqns. 5.38 - 5.41 into Eqn. 5.42, it yields

3

3
4

6
4

4
4

4

4

2

While from the composite theory,

4

cos 4

2

cos 2
	

(5.43)

can be expressed as
cos

4

4

2

cos

4
sin

sin

(5.44)

With further simplification, it can be found that Eqn. 5.43 is identical to Eqn. 5.44.

5.4 Laminated plate under static loading
162

In this section, it is to verify the proposed peridynamic model with an anlytical solution. A
simple tensile test is performed on a laminated plate. The peridynamic results will be compared
with the results obtained from the composite theory.

Consider a 100	
tensile pressure of 10	

100	

laminated plate with fibers in

direction as shown in Fig. 5.4. A

is applied at the bottom and the top of the plate. The plate is made of

E-Glass/Epoxy and the material properties are shonw in Table 5.3.

Table 5.3 Material propteties of E-Glass/Epoxy
Longitudinal Young’s modulus,
Transverse Young’s modulus,
Poisson’s ratio,
Shear modulus,
Mass density,

41
10.4
0.28
4.3
1970

/

Based on the composite theory [6], the components of the complicance matrix are
1/

(5.45)

1/

(5.46)
/

1/

Strains from the composite theory can be calcuated as follows

163

(5.47)
(5.48)

(5.49)

The transformed compliance matrix is caculated as
0
0

	
0

(5.50)

0

where
2
2
and

cos

and

(5.51)

sin .

From Eqn. 5.49, the displacement field can be obtained from the composite theory.
Displacements from peridynamics are compared with those from the composite theory in Fig.
5.5-Fig. 5.10 with

0°, 45°	and 60°. As can be seen, the results from peridynamics and those

from composite theory are identical to each other.

5.5 Free vibration of a laminated beam
The free vibration of a laminated beam is investigated in this section by the Classical Laminated
Beam Theory and the solution will be used to verify that obtained from peridynamic model.

Consider a simply supported beam as shown Fig. 5.11. The length of the beam is
the thickness of the beam is

10	

, where

made of Kevlar/Epoxy with fibers oriented in

∙

and

is the aspect ratio of the beam. The beam is
direction. The material properties are shonw in

164

Table 5.4.

Table 5.4 Material propteties of Kevlar/Epoxy
Longitudinal Young’s modulus,

80

Transverse Young’s modulus,

5.5
0.34

Poisson’s ratio,

2.2

Shear modulus,

1380

Mass density,

/

From [7], the governing equations of the beam are
°

°
°

0
0	

(5.52)

(5.53)

where
/
/

(5.54)

/
/

0

/
/

(5.55)
/12	

(5.56)

A solution satisfying the governing equations is
cos

(5.57)

sin

(5.58)

165

Eqn. 5.57 and Eqn. 5.58 satisfy the simply supported boundary conditions automatically.
Substituting Eqn. 5.57 and Eqn. 5.58 into Eqn. 5.52 and Eqn. 5.53, it yields
0

(5.59)

0

(5.60)

Solving Eqn. 5.59 and Eqn. 5.60, it can be concluded that
(5.61)
0

(5.62)

If the initial condition of the beam is
0

1

10

sin

(5.63)

then
1

10

(5.64)
(5.65)

The solution of the problem should have the following forms
,
, ,
where

,

and

sin

cos

(5.66)

cos

cos

(5.67)

can be found from Eqn. 5.64, Eqn. 5.65 and Eqn. 5.61, respectively.

The same problem can be simulated by peridynamics. The displacement history of point A and
point B (Fig. 5.11) are recorded. Fig. 5.12 compares the peridynamic results with those from the

166

composite beam theory for the aspect ratio

5. Fig. 5.13 compares the peridynamic results

with those from the composite beam theory for the aspect ratio
methods show good match in vertical displacement

20. Results from the two

of point B. For horizontal displacement

of point A, peridynamic result is almost the same as the beam theory result when the aspect ratio
20. However, there is difference between the peridynamic result and the beam theory result
5. This is because the beam theory assumes no variation of vertical

when the aspect ratio

displacement when the beam is slender. With a small aspect ratio, such as

5, this variation is

not neglible and the beam theory does not provide a good approximation.

5.6 Comparisons of crack propagation velocity with experiments
Experiment studies on dynamic crack propagation in fiber-reinforced composite materials have
been conducted by Zheng[9], Rosakis[8], Stout[10] and Coker[11,12]. It has been shown that a
weak fracture plane usually occurs between fiber and matrix in unidirectional fiber-reinforced
composites. Due to material anisotropy, the wave speed along the fiber direction is very different
from that along the perpendicular direction. Dynamic crack propagation has been commonly
investigated by finite element method. A limited number of computational studies have been
reported by Huang [13], Hwang [14], Kumar [15], Stout [10], Lo [16], Sun [17] and Pandey[18].
The limit of computational works is likely due to the requirement of remeshing and the
complexity of handling elements once crack starts. In this section, dynamic crack propagation in
a unidirectional graphite/epoxy composite is studied with the use of peridynamics. The
computational results are validated by the experiment results given in refefrence [8].

Consider a 76	

152	

unidirectional graphite/epoxy fiber-reinforced composite plates

167

under three-point bending [8] as shown in Fig. 5.14. The fiber is in 90° direction. The material
properties are shown in Table 5.5 [8, 22]. There is a notch with a length of 15.2	

, i.e. 20% of

the plate width, at the left boundary of the plate. This crack length is used because it was used in
the past to produce reliable results in dynamic fracture experiements [19]. To minimize residual
stresses due to machining, a low-speed diamond saw was used to produce the initial notch with a
widtch of approximately 1.5	

.

Table 5.5 Material propteties of graphite/epoxy
Longitudinal Young’s modulus,

150

Transverse Young’s modulus,

11.6
0.36

Poisson’s ratio,

3.5

Shear modulus,
Mode I intralaminar fracture energy for longitudinal loading,
Mode I intralaminar fracture energy for transver loading,

77.9 /
5 /
1590

Mass density,

/

A drop-weight tower is used to introduce impact on the opposite side of the notch with an
impacting speed of

4	 / . After the impact, stress waves propagate to the interior of the

plate and then reflects from the boundaries. Because of the anisotropy of the material, stress
waves travel in different directions at different velocities. Experiments show that the crack starts
to propagate at about 25

after the impact so the effects of dispersion are not very important

since the applied stress pulse is very long (about 120	

) compared with the time of crack

initiation. Therefore, loading is continuously applied throughout the entire event.

168

The real-time visualization of dynamic fracture is produced by an optical method of coherent
gradient sensing (CGS) in reflection [20, 21]. Imaging is performed with a rotating-mirror highspeed camera. Details of the CGS system can be found in [8, 20, 21].

Fig. 5.15 shows the crack propgation velocity from the peridynamic simulation. The initial ime
0 is used to denote the beginning of the crack propagation. For negtive time,
The crack starts from about 700	 / and accelerates to 900	 / within about 10	
decelerates to less than 500	 / in 40	

0	 / .
. It then

after the initiation. This deceleration is believed to be

due to the fact that the growing crack tip enters a region of high compressive stress zone as it
appoaches the loading area. The peridynamic computational results are compared with the
experiment results from Ref. [8]. As shown in Fig. 5.15, they agree each other reasonably well.

5.7 Dynamic fracture mode in unidirectional composites
In order to investigate the behavior of cracks, Wu [24] conducted experiments with
unidirectional, fiberglass-reinforced Scotch composites with a centered precrack in the direction
of fibers. The composites were loaded with tension, pure shear and combined tension and shear.
In all three cases, it was observed that the crack propagated in the same direction as the fiber
direction. Finite element analysis were also used to study the damage path and failure initiation
of pre-notched composites by Boger [25] and Satyanarayana [26]. They prediected damge in
composite plates notched in the center for different layups under tension. Both experiment results
and simulation results showed that the crack path and failure initiation depends on fiber
orientation.

169

In this section, the crack propagation path and dynamic fracture mode of unidirectional fiberreinforced composites are studied using the proposed peridynamic model. The qualitive
comparison of the peridynamic results with those from experiments are of interest.

Consider the compact tension test on a 100	
composite plate with a 20	

200	

carbon/epoxy unidirectional

pre-notch at the center as shown in Fig. 5.16. The plate is loaded

at the top and the bottom boundary by a uniform stress . The fiber orientation is . The material
properties of the carbon/expoxy plate are shown in Table 5.6 [23].

Table 5.6 Material propteties of carbon/epoxy
329

Longitudinal Young’s modulus,
Transverse Young’s modulus,

6

Poisson’s ratio,

0.346
4.4

Shear modulus,
Mode I intralamina fracture energy for longitudinal loading,

15.49 	 /

Mode I intralamina fracture energy for transver loading,

0.168 	 /

Mass density,

1630

/

Fig. 5.17, Fig. 5.18, Fig. 5.19, Fig. 5.20 and Fig. 5.21 show the peridynamic simulation results
for

0°,	

30°,	

45°,	

60° and

90°, respectively. In each figure, there are three

contour plots with (a) vertical diaplacement of the plate, (b) strain energy density distribution and
(c) local damge of the plate. Different color bars are associated with different plots with red

170

indicating the highest value while blue the lowest value. Crack paths can be clearly seen from the
strain energy density contour plot as there are always strain energy concentraions at the crack
tips. The local damage is defined by Eqn. 2.12 and Enq. 2.13, which show different levels of
damage. A proper cuttoff value can be defined to judge if there is a crack. In all cases, the crack
propagates in the same direction as the fibers, which is consistent with the experiment
observations from Wu [24]. The damage is due to the seperation between matrix and fiber. There
is no fiber breakage.

As expected, in the

0° case, the crack propagates in the same direction as the pre-notch.

This also matches with the computational and experiment results in Section 5.6. In the smaller
angle case

30°, aside from the major crack, which propagates in 30°, there is matrix

shattering at the sides of the plate in 0° direction. The matrix shattering happens before the crack
starts to propagate as shown in Fig. 5.22. It starts at the lateral of the plate and propagates to the
interior of the plate. From Fig. 5.18 (c), the matrix shattering is not as severe as the major crack
since only 20% of the bonds are broken. However, in the major carck, more than 70% of the
bonds are broken. This is why matrix shattering is only a material softening and may not be seen
from experiment observation as reported in [24]. For the

90° case, the composite plate fails

due to splitting caused by shear stress in the matrix. It which matches with the findings of Boger
[25].

5.8 Conclusions
A peridynamic orthotropic material model baesd on the beam model is proposed in this chapter.
There are four independent mateiral parameters in this model and it matches with the four

171

material properties for two-dimensional orthotropic materials. The bond material properties
depend on these four mateiral parameters and the angle between the bond orientation and fiber
orientation. This results in the continuity of

the bond stiffness function with no need of

remeshing for different fiber orientations. This model is verified by a static tensile test and a
vibration problem of a laminated beam.

Dynamic damage propagation problems in composite materials can br greatly benefited from
peridynamics. The prediction of damage initiation and crack proagation of composite materials is
complex using traditional methods, such as finite element analysis, due to its anisotropy. As
investigated in peridynamic simulations, there is no need of tracking each crack propagation,
finding different damage modes and applying different damage rules. Damage happes
automatically. A single edge notch test is simulated in this chapter and the results match with the
experiment results. Crack path and failure initiation of a center notch plate is predicted
successfully by peridynamics when comparing with the experiment results.

172

Figure 5.1 A composite plate with fiber in ° direction and a bond in ° direction.

173 
 

Figure 5.2 Coordinate system for calculating strain energy density at point .

174 
 

Figure 5.3 Beam model for orthotropic materials.

175 
 

Figure 5.4 Pulling test in a composite laminate.

176 
 

Figure 5.5 Comparison of

of 0° laminate calculated from peridynamics (top) and composite

theory (bottom).

177 
 

Figure 5.6 Comparison of

of 45° laminate calculated from peridynamics (top) and the

composite theory (bottom).
178 
 

Figure 5.7 Comparison of

of 60° laminate calculated from peridynamics (top) and the

composite theory (bottom).
179 
 

Peridynamics

Figure 5.8 Comparison of

of 0° laminate calculated from peridynamics (top) and the

composite theory (bottom).
180 
 

Figure 5.9 Comparison of

of 45° laminate calculated from peridynamics (top) and the

composite theory (bottom).
181 
 

Figure 5.10 Comparison of

of 60° laminate calculated from peridynamics (top) and the

composite theory (bottom).
182 
 

Figure 5.11 Free vibration of a laminated beam with fibers in

183 
 

direction.

Figure 5.12 Peridynamic results compared with composite beam theory results with
Top: horizontal displacement

of point A. Bottom: vertical displacement

184 
 

of point B.

5.

Figure 5.13 Peridynamic results compared with composite beam theory results with
Top: horizontal displacement

of point A. Bottom: vertical displacement

185 
 

of point B.

20.

Figure 5.14 An unidirectional composite plate with single edge notch under three point bending.

186 
 

Figure 5.15 Comparison of crack propagation velocity between peridynamics and experiment.

187 
 

Figure 5.16 Compact tension test for a unidirectional composite plate.

188 
 

(a)

(b)

(c)

Figure 5.17 Simulation results for

50	

and

0°:

(a) vertical displacement, (b) strain energy density, (c) local damage.

189 
 

(a)

(b)

(c)

Figure 5.18 Simulation results for

70	

and

30°:

(a) vertical displacement, (b) strain energy density, (c) local damage.

190 
 

(a)

(b)

(c)

Figure 5.19 Simulation results for

70	

and

45°:

(a) vertical displacement, (b) strain energy density, (c) local damage.

191 
 

(a)

(b)

(c)

Figure 5.20 Simulation results for

90	

and

60°:

(a) vertical displacement, (b) strain energy density, (c) local damage.

192 
 

(a)

(b)

(c)

Figure 5.21 Simulation results for

100.5	

and

90°:

(a) vertical displacement, (b) strain energy density, (c) local damage.

193 
 

(a)

(b)

(c)

Figure 5.22 Simulation results for

40	

and

30°:

(a) vertical displacement, (b) strain energy density, (c) local damage.

194 
 

REFERENCES

195 
 

References

[1] Dwivedi, S. K., Espinosa, Modeling dynamic crack propagation in fiber reinforced
composites including frictional effects, Mechanics of Materials, 35, 481-509 (2003)
[2] Xu, J., Askari, A., Wechner, O., Damage and Failure Analysis of Composite Laminates under
Biaxial Loads, 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and
Materials Conference, April 23-26, 2007, Honolulu, Hawaii
[3] Hu, W., Ha, Y.D., Bobaru, F., Peridynamic model for dynamic fracture in unidirectional
fiber-reinforced composites, Comput. Methods Appl. Mech. Engrg., 217-220 (2012)
[4] Silling, S. A., Reformulation of elasticity theory for discontinuities and long-range forces, J.
Mech. Phys. Solids, 48, 175-209 (2000)
[5] http://www.wolfram.com/mathematica/
[6] Daniel, I. M., Ishai, O., Engineering mechanics of composite materials, the Second Edition,
Oxford University Press, Inc. (2006)
[7] Whitney, J. M., Structural Analysis Of Laminated Anisotropic Plates, Technomic Publishing
Company, Inc. (1987)
[8] Lambros, J., Rosakis, A. J., Dynamic crack initiation and growth in thick unidirectional
graphite/epoxy plates, Composite Science and Technology, 57, 55-65 (1997)
[9] Zheng, S., Sun, C.T., A double-plate finite-element model for the impact-induced
delamination problem, Composite Science and Technology, 53, 111-118 (1995)
[10] Stout, M.G., Liu, C., Ellis, R.W., Haberman, K.S., Bennet, J.G., Williams, Addesio, F.L.,
Rosakis, A.J., Mechanics analysis and simulation of a graphite/epoxy composite plate,
Experimental Mechanics, 187-192 (1998)
[11] Coker, D., Rosakis, A.J., Experimental observations of intersonic crack growth in
asymmetrically loaded unidirectional composite plates, GALCIT-SM Reports 98-16. CALTECH,
Pasadena, CA (1998)
[12] Coker, D., Rosakis, A.J., Experimental observations of intersonic crack growth in
asymmetrically loaded unidirectional composite plates, Philosophical Magazine, 81, 571-595
(2001)
[13] Huang, Y., Wang, W., Liu, C., Rosakis, A.J., Analysis of intersonic crack growth in
unidirectional fiber-reinforced composites, Journal of the Mechanics and Physics of Solids, 47,
1893-1916 (1999)

196 
 

[14] Hwang, C., Geubelle, P., A spectral scheme to simulate dynamic fracture problems in
composite, Computer Modeling in Engineering and Science, 4, 45-56 (2000)
[15] Kumar, P., Kishore, N.N., Initiation and propagation toughness of delamination crack under
an impact load, journal of the Mechanics and Physics of Solids, 46, 1773-1787 (1998)
[16] Lo, C. Y., Nakamura T., Kushner, A.S., Dynamic failure analysis along interfaces in
composite materials, Advanced Computational Methods for Material Modeling, AMD-vol.180,
ASME, 115-124 (1993)
[17] Sun, C. T., Qian, W., The use of finite extension strain energy release rates in fracture of
interfacial cracks, International Journal of Solids and Structures, 34, 2595-2609 (1997)
[18] Pandey, R.K., Sun, C.T., Calculating strain energy release rate in cracked orthotropic beams,
Journal of Thermoplastic Composite Materials, 9, 381-395 (1996)
[19] Tippur, H.K., Krishnaswamy, S., Rosakis, A. J., Optical mapping of crack tip deformations
using the method of transmission and reflection Coherent Gradient Sensing: A study of crack tip
K-dominance, Int. J. Fract., 52, 91-117 (1991)
[20] Tippur, H. V., Krishnaswamy, S., Rosakis, A. J., A coherent gradient sensor for crack tip
measurements: Analysis and experimental results, Int. J. Fract., 48, 193-204 (1991)
[21] Rosakis, A. J., Two optical techniques sensitive to gradients of optical path difference: The
method of caustics and the coherent gradient sensor (CGS), Experimental Techniques in Fracture,
J. Epstein, VCH, New York, 327-425 (1993)
[22] Donaldson, S. L., Fracture toughness testing of graphite/epoxy and graphite/PEEK
composites, Composites, 16, 103-112 (1985)
[23] Jose, S., Kumar, R. R., Jana, M. K., Venkateswara, R., Intralaminar fracture toughness of a
cross –ply laminate and its constituent sub-laminates, Composites Science and Technology, 61,
1115-1122 (2001)
[24] Wu, E. M., Fracture mechanics of anisotropic plates, Composite material workshop,
Lancaster, PA, Technomic Publishing Co. Inc, 20-43 (1968)
[25] Boger, P. B., Satyanarayana, A., Chunchu, P. B., Comparison of damage path prediction for
composite laminates by explicit and standard finite element analysis tools, 47th
AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference,
Newport, RI, May 1-4, 2006
[26] Satyanarayana, A., Bogert, P. B., Chunchu, P. B., The effect of delamination on damage
path and failure load prediction for notched composite laminates, 48th
AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference,
Honolulu, HI, April 23-26, 2007
197 
 

Chapter 6
Conclusions and Recommendations

6.1 Conclusions
The governing equations of commonly used continuum mechanics are of differential equations.
In simulating damage process, special techniques and updated boundary conditions must be
applied on the crack surfaces. Since the differential governing equations require displacement
continuity in the domain of study, they become troublesome when damage takes place. In
modeling dynamic crack propagation, knowledge of crack positions and propagation directions
are also required. It then is difficult and often times impossible for studying dynamic crack
propagation of fiber-reinforced composite materials based on the traditional continuum
mechanics. The following conclusions can be drawn from this dissertation research.

1. Instead of differential equations, peridynamics uses integral equations. Similar to molecular
dynamics, peridynamics assumes that the domain of interest is organized by points. Each point
interacts with every other point within the horizon through a bond. Damage at a point takes place
when a critical number of bonds associated with the point are broken. There is no need to impose
a separate damage theory such as fracture mechanics used in continuum mechanics. Chapter 2
introduced the governing equation of peridynamics and some basic properties of the
peridynamics.

2. Being a novel method, peridynamic theory is still in its infant stage. New models are needed
and require analytical verifications and experiment validations. In Chapter 3, the governing

198 
 

equation and a bond function for one-dimensional problems were proposed. The model was
verified by the analytical solution of a wave propagation problem. The analytical solution was
also used to study the convergence of numerical method and a balanced set of parameters were
identified. Split Hopkinson’s pressure bar (SHPB) tests were conducted and the experiment
results were used to validate the one-dimensional peridynamic model. Cross-interface bond
models were proposed to study the interface related issues and contact problems of peridynamics.
Both elastic model and plastic model were proposed so the peridynamic theory could be used to
simulate SHPB tests in validating both elastic and plastic deformations. Several issues involved
in the experiments were discussed and it was shown that peridynamic results matched well with
the experiment results. Finally, using the plastic model, the shaping effects in the SHPB tests
were simulated and the results suggested that the peridynamic model could be used for designing
the shaper. The one-dimensional model was validated in this chapter and it was concluded that
the impact process could be successfully simulated by peridynamics with its long range forces
without involving contact mechanics.

3. In Chapter 4, a bar model and a beam model were proposed for two-dimensional problems.
The two-parameter beam model was capable of presenting two independent material properties
required by isotropic materials. The proposed model was verified by a two-dimensional vibration
analysis. Compared with the solutions from classical mechanics, numerical parameters were
studied to show the convergence of peridynamic solutions to classical mechanics solutions. Two
computational techniques were also proposed and implemented. The computational efficiency
was greatly improved using these techniques. The ability for executing two-dimensional analysis
was thus proved feasible. Three dynamic damage propagation experiments from the literature

199 
 

were also used to validate the peridynamic model and its damage theory. In the first experiment,
the crack velocity in the single edge notch (SEN) test was simulated by peridynamics and the
results matched well with those from experiment given in the literature. In the second experiment,
peridynamics was shown to be capable of producing crack branching as reported in literature. In
the third experiment, an impulse loading was simulated and peridynamic results showed the
same crack propagation direction as given in the experiment. It was shown that the damage
propagated in peridynamic simulation without any additional damage theory. There was no need
of tracking individual crack paths or using special elements.

4. The orthotropic material models and their applications were studied in Chapter 5. A fourparameter peridynamic model for orthotropic materials was proposed to coincide with the four
independent material properties required for orthotropic materials. This four-parameter model
was different from the previous developed ones as it had four parameters and was mesh
independent. The model was verified by a static tensile test and a vibration excitation of a
laminated beam. An SEN (single edge notch) test of a 0° laminated plate was simulated by
peridynamics and the computational results matched with published experiment results. Fracture
initiation and crack path of laminated plates with different fiber orientations were also studied
using peridynamics. The mesh-free peridynamic model was convenient and efficient since there
was no need to have different meshes for different fiber orientations. Without using special
meshes and requiring prior knowledge of fracture paths, a general peridynamic code was able to
predict fracture velocity and crack path successfully.

200 
 

In conclusion, a peridynamic code was generalized in this work using Matlab. Two
computational techniques were proposed and implemented to improve the computational
efficiency dramatically. This rendered simulations of large domains feasible and paved a way for
peridynamics in simulating even larger structures in the future. Besides, material models for onedimensional problems, two-dimensional problems and orthotropic materials were proposed and
they matched well with the two material properties in isotropic materials and the four material
properties in orthotropic materials, respectively. Elastic models, plastic models and damage
theories of peridynamics were proposed and studied in this thesis work. These models were
compared with either analytical solutions or experiment results. All the comparisons showed
good matches. Finally, it is worth noting that the integral governing equation of peridynamics
covers the entire domain regardless the continuity or discontinuity of the displacement.
Therefore, as shown in this work, there is no need to apply special techniques or models on
displacement discontinuity. Cracks or material weakening can grow automatically whenever and
wherever they have to grow.

6.2 Recommendations
1. A three-dimensional peridynamic model can be developed in the future to simulate threedimensional problems. Silling’s bar model [1] is only able to model materials with Possion’s
ratio of ¼. The three-dimensional model can be developed from the beam model in Chapter 4 to
accommodate two material properties for isotropic materials. With the help of this threedimensional model, a three-dimensional axisymmetric peridynamic model can be developed
from the axisymmetric model in Chapter 4. The experiment data in Chapter 3 can be
conveniently used to validate this model.

201 
 

2. More complex constitutive functions can be proposed to accommodate material nonlinearity.
It has been shown that experiment based stress-strain data can be used as the computational input
for the peridynamic code. More experiments can be performed to study material properties and
the experiment data can help refining peridynamic models. Strain-rate dependent materials can
also be modeled by peridynamics. Besides the spring model, elastic or plastic, a dash pot can be
added to the constitutive bond function. Since a dash pot is velocity dependent, a material model
with a dash pot can be strain rate dependent.

3. Multi-physical loads can be modeled using peridynamics. Currently, peridynamics is only
used in solid mechanics. Gas and liquid can be modeled using similar governing equations as the
peridynamics for solid mechanics. The constitutive functions of gas and liquid, however, will be
greatly different from those for solids. Other than the mechanical loads, thermal loads, acoustic
loads and electromagnetic loads can also be modeled in peridynamics. In peridynamics, all
external loads are applied as body forces. The effects of all non-mechanical loads need to be
studied and can be applied as mechanical effects in the peridynamic domain.

4. Once a three-dimensional peridynamic model is developed, a three-dimensional orthotropic
materials can be modeled similar to the two-dimensional orthotropic material models in Chapter
5. This model can be conveniently used to model woven composite, which is a difficult subjects
in commonly used finite element analysis. When modeling woven composite, the weaving
structures of the woven composite need to be known to set up the material properties of bonds in

202 
 

different directions. Deformations and damages of each bond are independent and crack can
grow in any directions.

5. Currently, the peridynamics code is developed using Matlab. The domain grid is generalized
manually. It will be better to use commercial mesh generalizers to define problem domains and
generalize mesh in any size. Moreover, it will be better to embed the peridynamics code in a
commercial software because commercial software are robust and likely more efficient.

6. Validations of new developed theories and models are necessary. Problems requiring
peridynamic analysis are likely complex and it is almost impossible to find analytical solutions
for these problems. To develop more complex peridynamic models, more experiments
specifically designed for validating the peridynamic models are necessary.

203 
 

 

REFERENCES

204 
 

 

References

[1] Silling, S. A., Reformulation of elasticity theory for discontinuities and long-range forces, J.
Mech. Phys. Solids, 48, 175-209 (2000)

205