ARL-RP-0531 ● AUG 2015

US Army Research Laboratory

A Finite Difference-Augmented Peridynamics 

Method for Wave Dispersion

by Raymond A Wildman and George A Gazonas

A reprint from International Journal of Fracture. 2014; 190:39-52.

Approved public release; distribution unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the 

Army position unless so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official 

endorsement or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
ARL-RP-0531 ● AUG 2015

US Army Research Laboratory

A Finite Difference-Augmented Peridynamics 

Method for Wave Dispersion

by Raymond A Wildman and George A Gazonas

Weapons and Materials Research Directorate, ARL

A reprint from International Journal of Fracture. 2014; 190:39-52.

Approved for public release; distribution unlimited.
ii

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing the 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid 

OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

Aug 2015

2. REPORT TYPE

Reprint

3. DATES COVERED (From - To)

January 2014 –October 2014

4. TITLE AND SUBTITLE

A Finite Difference-Augmented Peridynamics Method for Wave Dispersion

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Raymond A Wildman and George A Gazonas

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

US Army Research Laboratory

ATTN: RDRL-WMM-B

Aberdeen Proving Ground, MD 21005-5069

8. PERFORMING ORGANIZATION REPORT NUMBER

ARL-RP-0531

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

Approved for public release; distribution unlimited.

13. SUPPLEMENTARY NOTES

A reprint from International Journal of Fracture. 2014; 190:39-52.

14. ABSTRACT

A method is presented for the modeling of brittle elastic fracture which combines peridynamics and a finite difference method

to mitigate the wave dispersion properties of peridynamics. Essentially, a finite difference method is used in the bulk for 

wave propagation modeling, while peridynamics is automatically inserted in high strain areas to model crack initiation and 

growth. The dispersion properties of finite difference methods and discretized peridynamics are reviewed and the interface 

reflection properties between the two regions are investigated. Results show that the augmented method can improve the 

modeling of wave propagation and boundary conditions. In addition, the numerical stress intensity factor computed at a crack 

tip shows reduced oscillations in the augmented method, likely due to the improved dispersion properties of the bulk. 

Dynamic fracture simulations show a difference in crack paths between the methods.

15. SUBJECT TERMS

Finite difference; peridynamics; fracture; wave propagation

16. SECURITY CLASSIFICATION OF:

17. LIMITATION

 OF

 ABSTRACT

UU

18. NUMBER

 OF 

 PAGES

20

19a. NAME OF RESPONSIBLE PERSON

Raymond A Wildman

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

410-306-2232

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
Int J Fract (2014) 190:39–52

DOI 10.1007/s10704-014-9973-1

OR IG INAL PAPER

A finite difference-augmented peridynamics method

for reducing wave dispersion

Raymond A. Wildman · George A. Gazonas

Received: 9 July 2014 / Accepted: 7 October 2014 / Published online: 21 October 2014

© Springer Science+Business Media Dordrecht (outside the USA) 2014

Abstract A method is presented for the modeling of

brittle elastic fracture which combines peridynamics

and a finite difference method to mitigate the wave

dispersion properties of peridynamics. Essentially, a

finite difference method is used in the bulk for wave

propagation modeling, while peridynamics is automatically inserted in high strain areas to model crack

initiation and growth. The dispersion properties of

finite difference methods and discretized peridynamics are reviewed and the interface reflection properties

between the two regions are investigated. Results show

that the augmented method can improve the modeling

of wave propagation and boundary conditions. In addition, the numerical stress intensity factor computed at a

crack tip shows reduced oscillations in the augmented

method, likely due to the improved dispersion properties of the bulk. Dynamic fracture simulations show a

difference in crack paths between the methods.

Keywords Finite difference · Peridynamics ·

Fracture · Wave propagation

1 Introduction

Computational mechanics problems involving the fracture and failure of materials are difficult to solve using

R. A. Wildman (B) · G. A. Gazonas

U.S. Army Research Laboratory, Attn: RDRL-WMM-B,

US Army Research Laboratory, Aberdeen, MD 21005,

USA

e-mail: raymond.a.wildman.civ@mail.mil

partial differential equation formulations as any discontinuity introduced by such failure causes difficulties in both the computation and formulation. Peridynamics proposes to address this issue by replacing the

spatial differential operators in the standard continuum

mechanics equations with integral operators that converge to those differential operators under certain conditions (Silling 2000; Silling et al. 2007; Emmrich and

Weckner 2006). When discretized and coupled with

a damage model, peridynamics has been shown to be

capable of matching a number of experimental results

in both brittle and ductile fracture problems (Silling

and Askari 2005; Ha and Bobaru 2010, 2011; Wildman and Gazonas 2013; Foster et al. 2010). In addition, peridynamics has been applied to other fields of

physics, such as heat transfer (Bobaru and Duangpanya 2010) and electromigration (Gerstle et al. 2008).

While peridynamics has seen success in matching fracture patterns observed in experiments, issues remain,

such as its behavior at material boundaries and its dispersion characteristics (Weckner and Abeyaratne 2005;

Mikata 2012). Here, we attempt to address dispersion

in peridynamic simulations by coupling with a finite

difference solver to compute internal forces.

Finite difference methods have a long history, beginning with a proof of their stability properties for hyperbolic partial differential equations (Courant et al. 1928)

and similar methods exist in many fields; a short

list includes electromagnetics (Alterman and Karal

1968; Zahradnik and Hron 1992), elasticity (Yee 1966;

Taflove and Brodwin 1975), and diffusion (Brian 1961).

123
40 R. A. Wildman, G. A. Gazonas

Their first implementations relied on rectangular, uniform grids, though finite difference methods have been

developed for unstructured grids and curved regions

(Hesthaven and Warburton 2002).

A finite difference method (FDM) will be augmented with peridynamics in this work, as the most

basic implementations of FD methods are point-based

on uniform grids, much like standard peridynamic discretizations. The compatibility of the spatial discretization allows for straightforward switching between

the two methodologies. Essentially, all that changes

between the two formulations is the calculation of internal forces: in FD methods, an approximation of the

divergence of the stress tensor is used, and in peridynamics, the integral of a micro-force function is used

over some finite horizon. In fact, peridynamics can be

viewed as a sum of continuous differences, and when

discretized appears much like a weighted sum of finite

differences. In addition, peridynamics is known to be

less accurate on boundaries so the finite difference discretization can be used to implement boundary conditions in undamaged regions.

Adaptive methods for peridynamics have been introduced to improve its accuracy (Bobaru et al. 2009;

Bobaru and Ha 2011). These methods use an adaptive grid to refine the solution around crack tips and

other stress concentrations. It has been shown that peridynamics is more accurate and converges to the continuum mechanics formulation as the horizon shrinks.

Of course, a shrinking horizon must be coupled with

a refined grid, hence the adaptive refinement method.

Here, we attempt to gain more accurate wave propagation using a different formulation in the bulk, rather

than adaptive refinement. As will be shown, a more

accurate dispersion relation leads to more accurate

stress concentrations around crack tips than peridynamics alone, without the need for adaptive refinement.

Currently, hybrid peridynamics/finite element methods have been developed (Lubineau et al. 2012; Azdoud

et al. 2013, 2014). The combination of peridynamics with a local finite element method has the same

goal as this work, though current literature is restricted

to quasi-static problems and the method requires the

local and non-local regions be set a priori. In addition, Seleson et al. (2013) describe a hybrid finite difference/peridynamics model using a blending function

in 1D, though again, the focus is on preset, unchanging local/nonlocal regions. In contrast, this work will

focus on the dynamic insertion of peridynamic nodes

into a finite difference domain. Previous approaches

have used blending functions to ease the reflections

from these disparate computational domains, though

here the focus will be on dynamic insertion.

The remainder of the paper proceeds as follows:

In Sect. 2, peridynamics is summarized and a specific

finite difference method is formulated. In addition,

the dispersion properties of both discretizations are

reviewed, and finally a hybrid approach is detailed. In

Sect. 3, numerical results are given, including a demonstration of the dispersion properties of both methods,

interface reflections between the two computational

regions, a comparison of stress intensity factors at a

stationary crack tip, and finally a dynamic crack propagation problem. Section 4 concludes the paper.

2 Formulation

In this section, first peridynamics and FDM are summarized, followed by a discussion of their numerical

dispersion properties. Subsequently, a hybrid method

is presented, which essentially uses finite difference

approximations on discretization nodes undergoing

small strain, and peridynamics in areas of high

strain.

2.1 Peridynamics

To facilitate comparisons of dispersion characteristics,

the peridynamics formulation is the linearization of

the standard bond-based method described in (Silling

2000). Bond-based peridynamics can be stated as the

integro-differential relation

ρ

∂2

∂t2 u =



Hx

f



u− u, x− x



dx+ b, (1)

where bold type indicates a vector, u = u (x, t), u=

u



x

, t



, ρ is the density (assumed to be homogeneous), Hx is the horizon or region of influence at point

x and is typically a sphere of radius δ, b is a body force,

and f is a peridynamic force function. The linearized

version of Eq. (1) can then be stated as

ρ

∂2

∂t2 u =



Hx

c

ξ ⊗ ξ

|ξ |

3 ηdx+ b, (2)

where ξ = x

−x, η = u−u and c is a material constant

that can be related to classical elasticity constitutive

123
A finite difference-augmented peridynamics method 41

parameters. For example, in 3D, c is related to Young’s

modulus E and Poisson’s ratio ν as

c = 6E

πδ4 (1 − 2ν) (3)

where ν = 1/4, in 2D plane stress we have

c = 6E

πδ3 (1 − ν)

, (4)

where ν = 1/3 and in 2D plane strain we have

c = 6λ

πδ3ν

, (5)

where ν = 1/4 and λ is Lamé’s first parameter, all of

which are derived via matching the elastic strain energy

of the classical and peridynamic formulations (Ha and

Bobaru 2010).

A standard, point-based discretization of Eq. (2), is

used and is given by

ρ

∂2

∂t2 um = 

Nm

n=1

c

ξ mn ⊗ ξ mn 

ξ mn





3 ηmnVmn + bm, (6)

where Nm is the number of nodes in the horizon Hx of

node xm, ξ mn = xn − xm, ηmn = u (xn, t) − u (xm, t),

and Vmn is the volume (or area in 2D) of the nth node

associated with node m (Silling and Askari 2005). This

distinction on the volume is necessary because exact

integration weights that account for the edge of the

horizon are used. In other words, the regions or cells

associated with nodes in the reference configuration are

assumed to be rectangular (with sides given by the discretization sizes Δx and Δy), and if the entire rectangle

resides within a node’s horizon then its area is simply

the area of the square (in 2D). If the region intersects

the horizon, then its partial area is computed exactly

by using the divergence theorem to convert the area

integral into a line integral that encloses the integration

region. The line integral will then be a collection of

zero or more straight line segments and a circular arc,

all of which can be computed exactly.

The temporal discretization is a velocity Verlet

method, stated as

vk+1/2 = vk +

Δt

2

ak ,

uk+1 = uk + Δtvk+1/2

,

ak+1 = (F + b)

ρ ,

vk+1 = vk+1/2 +

Δt

2

ak+1,

(7)

where a = ∂2u/∂t2 is the acceleration, v = ∂u/∂t is the

velocity, F is the internal force given as the first term

of the right hand side of Eq. (6), b is a given body force

(that can also be used to apply boundary loads), and Δt

is the time step size.

Most importantly, peridynamics implements material failure by separating “bonds,” i.e. removing nodes

from a given node’s horizon so that they do not interact. The simplest damage model in the literature is to

remove a bond based on a bond strain measure, the

linearized version of which is given as

smn = ηmn · ξmn 

ξmn





2 , (8)

with the definitions of ξmn and ηmn given above. A

maximum strain can be given as a constant, giving the

failure criterion

smn > smax, (9)

where smax is a maximal bond strain, typically given in

terms of a fracture energy as

smax =

4πG0

9Eδ , (10)

where G0 is an energy release rate (Silling and Askari

2005). If Eq. (9) is met for a given bond pair {m, n},

node n is removed from node m’s family (and vice

versa), and not included in the summation of Eq. (6).

After a bond breaks, its supported load is transferred to

other bonds in the nodes’ families, encouraging more

bond breakage and leading to progressive crack formation.

While peridynamics provides a description of the

elastic wave equation that is accurate under a converging horizon and its discretized form gives a useful

damage model, in general its dispersion characteristics are not consistent with a linear elastic material.

Other numerical methods can provide more accurate

wave propagation modeling and augment peridynamics’ superior handling of discontinuities.

2.2 Finite difference method

A finite difference formulation follows standard elasticity theory and discretizes its differential operators

with finite difference stencils. The linear elastic wave

equation can be stated as

ρ

∂2

∂t2 u =∇· σ + b, (11)

123
42 R. A. Wildman, G. A. Gazonas

where bold type with an overbar represents a secondrank tensor, b is a body force, the stress, σ is given

by

σ = C :  = C :

1

2 (∇u + u∇), (12)

two overbars represent a fourth-rank tensor, and C is

a fourth-rank constitutive tensor. In 3D, the isotropic,

linear elastic constitutive tensor is given by

C =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎣

2μ + λλ λ 000

λ 2μ + λ λ 000

λ λ 2μ + λ 000

0 0 02μ 0 0

0 0 0 02μ 0

0 0 0 0 02μ

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎦

,

(13)

where the Lamé parameters λ and μ are

λ = Eν

(1 + ν) (1 − 2ν)

, (14)

and

μ = E

2 (1 + ν)

. (15)

Now, the above set of equations can be specialized

with the restrictions on Poisson’s ratio inherent in peridynamics. For example, for 2D plane strain the constitutive tensor becomes

C = 2

5

E

⎡

⎣

310

130

002

⎤

⎦ , (16)

resulting in the coupled equations

ρ

∂2

∂t2 ux = 25

E



3

∂2

∂x2 ux+2

∂2

∂x∂y

uy+

∂2

∂y2 ux



+bx,

ρ

∂2

∂t2 uy = 25

E



3

∂2

∂y2 uy+2

∂2

∂x∂y

ux+

∂2

∂x2 uy



+by.

(17)

The remaining two formulations follow similarly, for

example, the plane stress formulation is identical but

with the factor of 2/5E replaced with 3/8E.

Often, finite difference discretizations for elasticity

will use two or four staggered grids, solving for the displacement on one grid and the stress on the second, or

different components of each on four grids (Virieux

1986). (This approach is known as the Yee cell in

electromagnetics (Yee 1966).) Though staggered grid

approaches can be more accurate, here we use a single grid to ensure compatibility with the peridynamics

grid. (In an extension of this method, state-based peridynamics could be implemented in a dual, staggeredgrid discretization with the more standard finite difference formulation.) Rather, the second order spatial

derivatives will be discretized with a standard central

difference approximation given as:

∂2

∂x2 f (x, y)

≈ f (x + Δx, y) − 2 f (x, y) + f (x − Δx, y)

Δx2 ,

(18)

∂2

∂y2 f (x, y)

≈ f (x, y + Δy) − 2 f (x, y) + f (x, y − Δy)

Δy2 ,

(19)

and

∂2

∂x∂y

f (x, y) ≈ f (x + Δx, y + Δy)

ΔxΔy

− f (x + Δx, y − Δy)

ΔxΔy

+ f (x − Δx, y − Δy)

ΔxΔy

− f (x − Δx, y + Δy)

ΔxΔy ,

(20)

where Δx and Δy are the discretization sizes in x and

y respectively. The temporal discretization is again the

velocity Verlet method shown in Eq. (7).

2.3 Numerical dispersion

Next, we discuss the numerical dispersion inherent in

both peridynamics and FDM, as any interface between

the two methods with different dispersion characteristics will result in wave reflections. Previous work has

shown that the peridynamics formulation displays dispersion, even before discretization (Weckner and Abeyaratne 2005) and it is well known that FDM suffers

from numerical dispersion. In 1D, the linearized form

of peridynamics behaves as a continuous sum of differences, much like discretized finite difference methods.

For example, specializing Eqs. 2 and 6 to 1D gives

ρ

∂2

∂t2 u =



Hx

c u− u

|x− x|

dx

, (21)

123
A finite difference-augmented peridynamics method 43

which can be discretized as

ρ

∂2

∂t2 um = 

Nm

n=1

c un − um

|xn − xm|

Vmn. (22)

A central difference scheme can be recovered by setting

the horizon to be equal to the spatial discretization size

δ = Δx, resulting in Vmn = Δx, Nm = 2 (the two

nodes adjacent to node m), c = E/δ2, and giving



1

n=−1 n=0

E

δ2

un − um

|xn − xm|

Δx = E um−1 − 2um + um+1

Δx2 .

(23)

A standard discretization of 1D peridynamics with

a horizon size equal to the spatial discretization size

is then equivalent to a central difference scheme, thus

motivating an augmented method. With this discretization, the best that can be done is to match the dispersion

of a finite difference technique, though it will be clear

that including a larger horizon size impacts dispersion.

First, the dispersion relation for a second order accurate

in space and time FDM can be given by

cos(ωΔt) = γ 2 cos(βΔx) − γ 2 + 1, (24)

where γ = vΔt/Δx is the Courant number, v = √E/ρ is the wave speed, ω is the angular frequency

and β is the wavenumber.

In 1D, the discretized peridynamic operator becomes

a sum of central differences in all of the nodes within

the horizon. The addition of these points results in a

different dispersion relation for each value of the horizon size δ. For example, with a horizon size of 2Δx

and a central difference discretization in time, we have

cos(ωΔt) = γ 2

2



cos(βΔx) +

1

4

cos(2βΔx)



− 5

8

γ 2 + 1, (25)

and for a horizon size of 3Δx we have

cos(ωΔt) = 2γ 2

9

cos(βΔx)

+ γ 2

9

cos(2βΔx) + γ 2

27

cos(3βΔx)

− 10

27γ 2 + 1, (26)

As an example, Fig. 1 gives a comparison of the dispersion curves of FDM and peridynamics with two horizon

sizes, δ = 2Δx and δ = 3Δx for Δx = 1 mm, Δt =

10 ns, E = 65 GPa, and ρ = 2,235 kg/m3.

0 500 1000 1500 2000 2500 3000

0

0.5

1

1.5

2

x 107

β (rad/m)

ω (rad/s)

Exact

FDM

PD2

PD3

Fig. 1 Dispersion curves for FDM and discretized peridynamics

with two horizon sizes

Again, as we wish to apply FDM to some areas of

the computational domain, and peridynamics to others,

these interfaces will exhibit frequency-dependent wave

reflections because of the differing dispersion relations.

The magnitude of these reflections will be investigated

in Sect. 2.4.1.

2.4 Augmented method

The augmented method essentially implements a

straight-forward FD discretization while simultaneously computing peridynamics-based bond strains as a

fracture criterion. Unless fracture occurs, peridynamics

is not used to compute internal forces, it is only used as

a fracture criterion. After bond breakage occurs, peridynamics is then used to compute the internal forces.

2.4.1 Wave reflection at PD/FDM interfaces

As mentioned in Sect. 2.3 FDM exhibits numerical dispersion and peridynamics has inherent dispersion properties, which can be exacerbated when discretized. Our

computational grid will contain areas that use both

FDM and peridynamics to compute wave propagation,

and so any interface between the two will result in

reflections. The reflection coefficients can be computed

based on the dispersion relations given in Eqs. (24),

(25), and (26) using the numerical phase velocity. Here,

we will simply give the relation in 1D, which would be

123
44 R. A. Wildman, G. A. Gazonas

representative of a normally incident plane wave in an

infinite medium. The reflection coefficient between two

materials with different impedances is given by

Γ = Z2 − Z1

Z2 + Z1

, (27)

where Z1 and Z2 are impedances of the two media. In

1D, the acoustic impedance is given by

Z = 

Eρ = ρv. (28)

If we use the phase velocity vp = ω/β and the dispersion relations from Eqs. 24, 25, and 26, we can substitute a numerical impedance into Eq. (27), resulting

in

ΓPD2 = ωPD2 − ωFDM

ωPD2 + ωFDM

, (29)

where

ωFDM = cos−1



Δt

−1



γ 2 cos(βΔx) − γ 2 + 1



(30)

and

ωPD2 = cos−1



Δt

−1



γ 2

2



cos(βΔx) +

1

4 cos(2βΔx)



− 5

8

γ 2 + 1

 . (31)

Similar results follow for other combinations of dispersion relations, three of which are plotted in Fig. 2 using

the same parameters as those in Fig. 1. As expected, the

wave reflection between computational grids becomes

worse with increasing wavenumber and horizon size;

0 500 1000 1500 2000 2500 3000

−0.35

−0.3

−0.25

−0.2

−0.15

−0.1

−0.05

0

β (rad/m)

Γ

FD−PD2

FD−PD3

PD2−PD3

Fig. 2 Reflection coefficients for normal incidence between

FDM and discretized peridynamic computational grids

however, this is only an illustrative example as in 2D

and 3D the reflection coefficient will differ for angle of

incidence and wave type (longitudinal vs. shear). We

do of course expect similar behavior in general: Higher

reflection coefficients for higher wavenumbers (or frequencies). Finally, while larger horizon sizes result in

larger reflections between the FD/peridynamics interfaces, it has been observed that larger horizons yield

more accurate crack features (Ha and Bobaru 2011).

Ultimately, a large enough horizon should be used to

model fracture, though continually increasing the horizon will result in less accurate behavior at computational interfaces. Consequently, in the dynamic fracture

results (Sect. 3.4) a horizon of 3Δx will be used.

2.4.2 Methodology

The augmented method proceeds as follows: The set

of nodes is split into two sets, nodes at which internal forces are computed using finite differences NFD

and those at which internal forces are computed using

peridynamics NPD. These node sets do not overlap

(NFD ∩ NPD = ∅) and together they comprise the

set of all nodes N(NFD ∪ NPD = N). At each node,

rm ∈ NFD, both the finite difference approximation

of the internal force (divergence of the stress) and

the peridynamic bond strains are computed, restated

here as:

Fm =

⎡

⎣

2

5 E



3 ∂2

∂x2 ux + 2 ∂2

∂x∂y uy + ∂2

∂y2 ux



2

5 E



3 ∂2

∂y2 uy + 2 ∂2

∂x∂y ux + ∂2

∂x2 uy



⎤

⎦

r=rm

,

(32)

where ux and uy are functions of r and t, and with

the appropriate finite differences applied to each differential operator. Simultaneously, at each node in NFD

the peridynamic bond strain is computed as given in

Eq. (8), though it is not included in the internal force

estimation until the failure criterion is met (see Eq. 9).

After some bond strain smn exceeds the failure criterion,

the internal force estimation switches to peridynamics,

again given as

Fm = 

Nm

n=1

c

ξ mn ⊗ ξ mn 

ξ mn





3 ηmnVmn, (33)

thus moving the node from the set NFD to the set NPD.

In this way, peridynamics is used similar to cohesive zones in finite element methods: After some field

123
A finite difference-augmented peridynamics method 45

measure exceeds a critical value, a new constitutive

response is used to model the internal force computation giving progressive failure/fracture modeling. In

this way, the finite difference force computation can

more accurately model the wave propagation in the

bulk, and peridynamics can be used in damaged areas.

This approach is summarized in Algorithm 1.

Data: E, ρ, δ,Δx, Δy,smax

Result: u

compute integration weights at each node for peridynamic

horizon, δ;

compute finite difference coefficients;

for each time step do

for each node m in NFD do

compute force Fm using finite differences;

compute bond strains, smn;

if any smn > smax then

break bond {m, n};

switch node m to peridynamics force

computation;

end

end

for each node m in NPD do

compute force Fm using peridynamics;

compute bond strains, smn;

if any smn > smax then

break bond {m, n};

end

end

velocity Verlet temporal integration to update u;

apply boundary conditions to u;

end

Algorithm 1: FDM augmented peridynamics method

2.4.3 Boundary conditions

Boundary conditions can be handled using either

method, peridynamics or finite differences. Dirichlet boundary conditions are straightforward in either

method. A traction-free, Neumann-like boundary condition in peridynamics is straightforward, and is implemented by excluding bonds in a horizon. In other

words, a traction-free boundary conditions is set up

when bonds are broken along a fracture. In the finite

difference domain, the Neumann condition can be

set by expanding the traction as a finite difference

of boundary and internal nodes, and solving for the

unknown boundary displacements in terms of the

known internal displacements. For example, on a rectangular domain aligned with the Cartesian axis, the

traction-free boundary condition on the top edge can be

stated as

yˆ · σ = 0, (34)

giving the coupled equations

∂

∂y

ux +

∂

∂x

uy = 0,

∂

∂x

ux + 3

∂

∂y

uy = 0,

(35)

where the factor of 3 in the second equation is a result

of the constitutive relation of Eq. (16). Applying finite

differences to the above, the unknown displacements on

the boundary can be written in terms of known internal

displacements as

Δyuy (x + Δx, y) − Δyuy (x − Δx, y)

+ 2Δxux (x, y) = 2Δxux (x, y − Δy),

Δyux (x + Δx, y) − Δyux (x − Δx, y)

+ 6Δxuy (x, y) = 6Δxuy (x, y − Δy), (36)

where a central difference has been used for the x derivative and a backward difference has been used for the

y derivative. A sparse matrix can be formed to solve

for the displacement at all boundary nodes after the

displacement update in Eq. (7) is used to compute the

internal displacements. Of course, in Eq. (35), the zeros

on the right hand sides can be replaced with a loading

function, taking care to include Young’s modulus and

the correct scaling factor for either plane strain or plane

stress.

Unfortunately, it has been found that any interface between the finite difference method and peridynamics on a boundary can lead to high enough

stress concentrations that spurious fracture can form.

For now, the boundary of the material will be set

to not accumulate damage other than pre-cracked

regions.

3 Numerical results

In this section, several numerical results will be presented. The first will focus on comparing the dispersion

effects of the finite difference discretization and discretized peridynamics. Next, reflections at the interface

between the two computational regions will be demonstrated. In Sect. 3.3, the stress intensity factor of a semiinfinite line crack in an infinite medium will be compared in both the hybrid and peridynamic approaches.

Finally, a dynamic fracture problem will be illustrated

in Sect. 3.4.

123
46 R. A. Wildman, G. A. Gazonas

3.1 Dispersion comparison

In this subsection, the dispersion effects of a finite

difference method and discretized peridynamics with

two horizons (2Δx and 3Δx) in a 2D plane strain

setting are compared. The constitutive parameters are

given by Young’s modulus E = 65 GPa and density

ρ = 2,235 kg/m3 and the discretization parameters are

given by x-step of Δx = 1 mm, y-step of Δy = 1 mm,

and time step size of Δt = 10 ns. Poisson’s ratio is

fixed at ν = 1/4 as discussed above. The computational region is a 200 mm-by-200 mm2 (spanning the

region 0 ≤ x, y ≤ 200 mm), chosen to be large enough

so boundary effects do not influence the results. A

plane wave with a Gaussian profile is used as an initial

condition—so as to compare with an exact solution—

given by

u (r, t = 0) = exp −7.5 × 104 

r · xˆ − 0.1

2



xˆ.

(37)

The simulation is run for 10 μs and ux is observed at

the points r = (70, 100) mm and r = (80, 100) mm.

First, Fig. 3 shows the x-directed displacement ux

vs. time at the two points listed above. The displacement measured at 80 mm is shown as the solid lines

and the displacement at 70 mm is shown as the dashed

lines. An “exact” solution (determined assuming a 1D

problem, though the numerical solutions were computed in 2D) is shown as the red curves, the finite difference approximated solution is shown as the black

0 200 400 600 800 1000

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

Time step

PD2 (x=70)

PD2 (x=80)

PD3

PD3

FD

FD

Exact

Exact

Fig. 3 The x-directed displacement versus time step measured

at two points, demonstrating dispersion effects

105 106 107 108 10−20

10−15

10−10

10−5

100

Frequency (Hz)

PD2 (x=70)

PD2 (x=80)

PD3

PD3

FD

FD

Exact

Exact

Fig. 4 Fourier transform of Fig. 3

105 106 107 108 10−12

10−10

10−8

10−6

10−4

10−2

Frequency (Hz)

Absolute error (m)

FD (x=70)

FD (x=80)

PD3

PD3

PD2

PD2

Fig. 5 Absolute error in the Fourier transform of Fig. 3measured

against the exact solution

curves, and two peridynamics approximated solutions

are shown in green and blue, using horizons of 3Δx and

2Δx respectively. As can be seen, peridynamics with a

higher horizon displays errors characteristic of dispersion effects, namely a lowered amplitude and oscillations. Next, the Fourier transform of the x-directed displacement was computed and is shown in Fig. 4, along

with the absolute error in each solution compared with

the “exact” solution (Fig. 5). These figures demonstrate

that, compared to the finite difference approximation,

higher frequencies in the peridynamics methods display larger errors, consistent with dispersion errors.

123
A finite difference-augmented peridynamics method 47

3.2 Interface reflections

Wave reflections between computational regions are

examined here. This was discussed in an analytical

form for 1D in Sect. 2.4.1, though here we perform

numerical tests of the 2D algorithms. In the first test,

the parameters of the problem are identical to those

of the previous subsection, though now the left half of

the domain will be computed with the finite difference

method and the right half will use peridynamics with a

horizon size of either 3Δx or 2Δx. The initial condition is again a Gaussian pulse, now centered at the line

x = 75 mm and the x-directed displacement is now

measured at the point r = (90, 100) mm.

As shown in Fig. 6, the reflected wave can be examined by also computing the solution to this problem

with a fully finite difference domain, and subtracting

the results, thus leaving only the reflected wave in the

left half domain. In Fig. 6, the blue curve is the result

using peridynamics with a horizon of 3Δx in the righthalf and the red curve shows the simulation with a

horizon of 2Δx. Given that the magnitude of the incident pulse was 0.5 m, the reflections at the interface

can be significant for the 3Δx horizon size. By maximum magnitude, the reflected wave is approximately

12%; however, peridynamics will only ever be used

in regions that have experienced fracture or failure, so

some amount of wave reflection would be expected

regardless.

Given the above, it may however also be possible

that the anomalies seen in Fig. 6 may be due to the

0 200 400 600 800 1000

−0.06

−0.04

−0.02

0

0.02

0.04

0.06

0.08

Time step

Absolute error (m)

PD3

PD2

Fig. 6 Absolute error between a fully FD solution and a solution

with half peridynamics, half FD

0 200 400 600 800 1000

−0.015

−0.01

−0.005

0

0.005

0.01

0.015

Time step

Absolute error (m)

PD3

PD2

Fig. 7 Absolute error between a fully FD solution and a solution

with a vertical line of peridynamic nodes

dispersion properties of the peridynamic region. To test

this, we again perform the above simulation in a FD

domain with peridynamics “window” consisting of a

single vertical line of nodes at the line x = 100 mm.

The smaller window will reduce the dispersion effects

of a wave propagating through a peridynamic region,

but still demonstrate the effects of interface reflection.

Figure 7 gives the results of this test, with the blue curve

showing peridynamics with a horizon of 3Δx and the

red curve having a horizon of 2Δx.

For a fully fractured region, a free surface is formed

and 100% of a wave is reflected, so as a final test,

we introduce a free surface with boundary conditions

implemented using finite differences to compare a fracture surrounded by peridynamic nodes. The same set up

is used as described above, though here the free surface

is placed at x = 100 mm. The initial pulse is centered

at the line x = 70 mm and the result shown in Fig. 8 is

measured at the point r = (90, 100) mm. Again, Fig. 8

shows the difference between the fully FD simulation

and the simulation using a peridynamic boundary with

a thickness of one horizon.

3.3 Stress intensity comparison

Next, the augmented method and peridynamics are

compared using the dynamic mode I stress intensity

factor at a semi-infinite, stationary crack in an infinite medium. This problem has an analytical solution

for the stress intensity factor for both stationary and

123
48 R. A. Wildman, G. A. Gazonas

0 200 400 600 800 1000

−0.1

−0.05

0

0.05

0.1

Time step

Absolute error (m)

PD3

PD2

Fig. 8 Absolute error between a fully FD solution and a solution

with a peridynamic border

moving cracks (Liu et al. 2011), though here we focus

on the stationary case. As given by Liu et al. (2011),

the stress intensity factor for this geometry can be

stated as

KI (t) = H (t − tc)

2σ∗

1 − ν

vl (t − tc) (1 − 2ν)

π ,

(38)

where H (·) is the Heaviside step function, vl is the

longitudinal wave speed, and σ∗ is the applied stress.

To simulate an infinite problem, a large region is used

with a horizontal pre-crack through the center and the

simulation is stopped before the reflected wave reaches

the crack tip. The material properties used give a longitudinal wave speed of vl = 7,332 m/s with Poisson’s ratio 1/4. Finally, the applied tensile stress is

40 MPa.

To estimate the stress intensity factor in the simulation, the method outlined by Chen (1975) will be

used. To summarize, the stress at the crack tip can be

expressed

σyy ≈ K∗

I (t)

√2πr



1 + sin

θ

2

sin

3θ

2



cos

θ

2

, (39)

where K∗

I (t)is a low-order estimate of the stress intensity,r is the distance to the crack tip, and θ is the angle.

The stress intensity factor KI (t) can then be estimated

by choosing nodes along a line with constant angle with

respect to the crack tip, computing the stress at those

points using finite differences, and constructing a least

squares approximation of the stress near the crack tip.

0 0.5 1 1.5 2 2.5

x 10−5

0

0.2

0.4

0.6

0.8

1

Time (s)

Normalized stress intensity

FD

PD3

Exact

Fig. 9 Normalized stress intensity factor comparison between

hybrid and peridynamics

More precisely, we choose nodes along the line θ = 0

(typically five nodes, starting with the fourth closest to

the crack tip) and use the higher order approximation

of the stress intensity

K∗

I ≈ KI + α (θ)r

1/2 + β (θ)r + γ (θ)r3/2

, (40)

and, according to Chen (1975), the coefficients of noninteger powers ofr, α (θ) and γ (θ), are zero for θ = 0.

A least squares approximation is then used to estimate

KI (t) (Chen 1975). Figure 9 shows the normalized

stress intensity factor (normalized by σ

√πa where

a = 150 mm is the crack length) with the exact solution and both peridynamics and the hybridized method.

The oscillations seen in peridynamics (most likely due

to wave dispersion) are reduced in the hybrid method.

3.4 Dynamic fracture

Finally, a dynamic crack propagation problem is simulated. The problem setup is that described in (Kalthoff

1974): A 150 mm horizontal crack in a 300 mm-by100 mm-by-9 mm glass plate is propagated with a tensile load of 3,040 N. No material properties are given

in (Kalthoff 1974), so we use properties given in (Ha

and Bobaru 2010) for soda-lime glass: Young’s modulus of 72 GPa, density of 2,440 kg/m3, and Poisson

ratio of 1/3 (restricted for plane stress). We also have

chosen a fracture energy value of 3 J/m2, which gives

a crack pattern that most closely matches the experiment in (Kalthoff 1974). While this fracture energy

123
A finite difference-augmented peridynamics method 49

Fig. 10 Comparison of crack paths between the hybrid method

(red) and peridynamics (blue)

Fig. 11 Magnification of the crack branching area of Fig. 10

appears somewhat low, it is close to the value of approximately 8 J/m2 found in the literature for crack initiation

(Döll 1975). Finally, in the peridynamics simulations,

the load is applied as a body force as described in Ha

and Bobaru (2009), and as a traction boundary condition in the hybrid method, with a traction of 1.13MPa.

First, discretization sizes of Δx = 1 mm, Δy =

1 mm, and Δt = 10 ns were used in both methods and

the simulations were run for 100 μs (10,000 time steps)

with a horizon size of δ = 3Δx. The final crack patterns are shown in Fig. 10, where a crack is taken to

be any node with damage >0.4 (as measured by the

ratio of the number of broken bonds to the number of

initial bonds) and the hybrid result is shown in red and

the peridynamics is shown in blue. Figure 11 shows a

magnified view of the crack branch: As can be seen,

both methods branch at approximately the same location, though they differ somewhat in their subsequent

crack paths.

The same problem was repeated, though now with

discretization sizes of Δx =0.5 mm and Δy =0.5 mm,

and with a horizon again of δ = 3Δx. Figure 12 shows

the results of the simulations, again with the hybrid

method in red and the peridynamics in blue. Here, the

Fig. 12 Magnification of the crack branching area with a finer

discretization

peridynamic method alone appears to have an earlier

crack initiation time. This is investigated further by

tracking the crack tip position and speed versus time.

The crack speeds were estimated (for the finer discretization) by tracking the damage as it propagated

across the domain and are shown in Fig. 13. The location of the crack was recorded at every 5μs (500 time

steps) (shown in Fig. 14) and the crack speed was estimated by computing the central difference of the result.

The hybrid method has a crack initiation time slightly

later than the peridynamics results and a crack speed

slightly slower than that of the peridynamics. These differences are most likely due to the lower stress intensity

factor observed (see Fig. 9) at the crack tip of the hybrid

method. This can be confirmed by comparing the end

points of the cracks in Fig. 12. Two features can be seen

in the crack speed data in the form of the two dips in

crack speed around 35 and 65μs. The first drop in speed

occurs as the crack branches and the second when the

branch angles change slightly at 225 mm.

It is also illustrative to compare both methods

between the two discretizations given above. To that

end, Fig. 15 shows the peridynamic method at both

of the above discretizations. It can be seen in the figure that the crack initiation occurs slightly earlier in

the finer discretization (shown in red). In contrast, the

hybrid method shows a more stable fracture path, with

both following a very similar path and initiation time

(Fig. 16). In both figures, the coarse discretization is

plotted in blue and the fine discretization is shown in

red. These differences are most likely caused by the

oscillations seen in the stress intensity factor, due to the

higher dispersion in the peridynamic method alone.

123
50 R. A. Wildman, G. A. Gazonas

Fig. 13 Comparison of crack speeds of the hybrid method (red)

and peridynamics (blue)

0 20 40 60 80 100

0.16

0.18

0.2

0.22

0.24

0.26

0.28

Time (μs)

Crack tip position (m)

FD

PD

Fig. 14 Comparison of crack tip locations of the hybrid method

(red) and peridynamics (blue)

Finally, the dynamic fracture problem was repeated

with both discretizations used above, but with a higher

load of 6080 N. In this case, a more complicated crack

path is observed, with a larger difference between the

methods. Figure 17 shows the result using discretization sizes of Δx = 1 mm, Δy = 1 mm. A major

difference between the two methods is the shallower

crack branch angle at the branch onset and the additional branch close to the center in the hybrid method.

Otherwise, they both feature multiple branches, consistent with crack patterns seen at higher loads (Bowden et al. 1967; Ha and Bobaru 2011). Figure 18 shows

Fig. 15 Comparison of crack paths at different discretizations

(coarse in blue and fine in red) using peridynamics

Fig. 16 Comparison of crack paths at different discretizations

(coarse in blue and fine in red) using the hybrid method

Fig. 17 Comparison of crack paths at a higher load and a coarse

discretization

the result using discretization sizes of Δx = 0.5 mm,

Δy = 0.5 mm. At this discretization, the peridynamics solution features an additional branch near the far

edge and two arrested cracks not seen at the coarse

discretization of the peridynamics method.

123
A finite difference-augmented peridynamics method 51

Fig. 18 Comparison of crack paths at a higher load and a fine

discretization

4 Conclusions

A hybrid finite difference/peridynamics method was

presented to improve the dispersion properties of standalone peridynamics. The hybrid method essentially

switches the spatial internal force computation between

finite differences (divergence of the stress) and peridynamics (integral of a micro-force function) after a

damage criterion is met. Linearized peridynamics can

be shown to model linear elastic materials with specific

Poisson ratios, and this fact is used to develop a matching finite difference scheme. Analysis was presented

showing the dispersion properties of discretized peridynamics in 1D in comparison to a straight-forward

finite difference discretization. Numerical results were

also shown to demonstrate the dispersion properties of

the two methods in 2D. Other results show that the

numerically-computed stress intensity factor for peridynamics can oscillate, which may lead to differing

crack paths.

This method has a number of possible extensions,

such as 3D, state-based peridynamics, and non-linear

materials. Additional work may be done as well to ease

the transition region between peridynamics and finite

differences, as reflections can occur due to the differing

material representations.

References

Alterman Z, Karal F (1968) Propagation of elastic waves in layered media by finite difference methods. Bull Seismol Soc

Am 58(1):367–398

Azdoud Y, Han F, Lubineau G (2013) A morphing framework to

couple non-local and local anisotropic continua. Int J Solids

Struct 50(9):1332–1341

Azdoud Y, Han F, Lubineau G (2014) The morphing method as a

flexible tool for adaptive local/non-local simulation of static

fracture. Comput Mech 54:1–12

Bobaru F, Duangpanya M (2010) The peridynamic formulation

for transient heat conduction. Int J Heat Mass Transf 53(19–

20):4047–4059

Bobaru F, Ha YD (2011) Adaptive refinement and multiscale

modeling in 2D peridynamics. Int J Multiscale Comput Eng

9(6):635–659

Bobaru F, Yang M, Alves LF, Silling SA, Askari E, Xu J (2009)

Convergence, adaptive refinement, and scaling in 1D peridynamics. Int J Numer Method Biomed Eng 77(6):852–877

Bowden FP, Brunton JH, Field JE, Heyes AD (1967) Controlled

fracture of brittle solids and interruption of electrical current. Nature 216:38–42

Brian PLT (1961) A finite-difference method of high-order accuracy for the solution of three-dimensional transient heat conduction problems. AIChE J 7(3):367–370

Chen YM (1975) Numerical computation of dynamic stress

intensity factors by a lagrangian finite-difference method

(the hemp code). Eng Fract Mech 7(4):653–660

Courant R, Friedrichs K, Lewy H (1928) Über die partiellen

Differenzengleichungen der mathematischen Physik. Math

Ann 100:32–74

Döll W (1975) Investigations of the crack branching energy. Int

J Fract 11:184–186

Emmrich E, Weckner O (2006) The peridynamic model in nonlocal elasticity theory. Proc Appl Math Mech 6:155–156

Foster JT, Silling SA, Chen WW (2010) Viscoplasticity using

peridynamics. Int J Numer Methods Eng 81(10):1242–

1258

Gerstle W, Silling S, Read D, Tewary V, Lehoucq R (2008) Peridynamic simulation of electromigration. Science 8(2):75–

92

Ha YD, Bobaru F (2009) Traction boundary conditions in peridynamics: a convergence study. Tech rep, University of

Nebraska-Lincoln

Ha YD, Bobaru F (2010) Studies of dynamic crack propagation

and crack branching with peridynamics. Int J Fract 162(1–

2):229–244

Ha YD, Bobaru F (2011) Characteristics of dynamic brittle fracture captured with peridynamics. Eng Fract Mech

78(6):1156–1168

Hesthaven JS, Warburton T (2002) Nodal high-order methods on

unstructured grids: I. Time-domain solution of Maxwell’s

equations. J Comput Phys 181(1):186–221

Kalthoff JF (1974) On the propagation direction of bifurcated cracks. In: Sih GC (ed) Dynamic crack propagation.

Springer, Dordrecht, pp 449–458

Liu ZL, Menouillard T, Belytschko T (2011) An XFEM/Spectral

element method for dynamic crack propagation. Int J Fract

169(2):183–198

Lubineau G, Azdoud Y, Han F, Rey C, Askari A (2012) A morphing strategy to couple non-local to local continuum mechanics. J Mech Phys Solids 60(6):1088–1102

Mikata Y (2012) Analytical solutions of peristatic and peridynamic problems for a 1d infinite rod. Int J Solids Struct

49(21):2887–2897

Seleson P, Beneddine S, Prudhomme S (2013) A force-based

coupling scheme for peridynamics and classical elasticity.

Comput Mater Sci 66:34–49

123
52 R. A. Wildman, G. A. Gazonas

Silling S (2000) Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids

48(1):175–209

Silling S, Askari E (2005) A meshfree method based on the

peridynamic model of solid mechanics. Comput Struct

83(17–18):1526–1535

Silling SA, Epton M, Weckner O, Xu J, Askari E (2007) Peridynamic states and constitutive modeling. J Elast 88(2):

151–184

Taflove A, Brodwin ME (1975) Numerical solution of steady–

state electromagnetic scattering problems using the timedependent Maxwell’s equations. IEEE Trans Microw

Theory Tech 23:623–630

Virieux J (1986) P-sv wave propagation in heterogeneous

media: velocitystress finite difference method. Geophysics

51(4):889–901

Weckner O, Abeyaratne R (2005) The effect of long-range forces

on the dynamics of a bar. J Mech Phys Solids 53(3):705–728

Wildman RA, Gazonas GA (2013) A perfectly matched layer

for peridynamics in two dimensions. J Mech Mater Struct

7(8):765–781

Yee K (1966) Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media.

IEEE Trans Antennas Propag 14(3):302–307

Zahradnik J, Hron F (1992) Robust finite-difference scheme for

elastic waves on coarse grids. Studia Geophysica et Geodaetica 36(1):1–19

123
15

1 DEFENSE TECH INFO CTR

(PDF) DTIC OCA

2 US ARMY RSRCH LAB

(PDF) IMAL HRA MAIL & RECORDS

MGMT

RDRL CIO LL TECHL LIB

1 GOVT PRNTG OFC

(PDF) A MALHOTRA

52 US ARMY RSRCH LAB

(PDF) RDRL CIH 

 J KNAP

RDRL WM

 B FORCH

 S KARNA

 J MCCAULEY

RDRL WML B

I BATYREV

 J BRENNAN

 E BYRD

 S IZVYEKOV

 W D MATTSON

 B RICE

 D TAYLOR

 N WEINGARTEN

RDRL WML H

 C MEYER

 B SCHUSTER

RDRL WMM

 J BEATTY

 R DOWDING

 J ZABINSKI

RDRL WMM B

 T BOGETTI

 C FOUNTZOULAS

 G GAZONAS

 D HOPKINS

 B LOVE

 B M POWERS

 T WALTER

 R WILDMAN

 C YEN

 J YU

RDRL WMM D

 S WALSH

RDRL WMM E

 J LASALVIA 

RDRL WMM F

 M TSCHOPP

RDRL WMM G

 J ANDZELM

 T CHANTAWANSRI

 C RINDERSPACHER

 T SIRK

 Y SLIOZBERG

RDRL WMP

 S E SCHOENFELD

RDRL WMP B

 S SATAPATHY

 A SOKOLOW

 T WEERASOORIYA

RDRL WMP C

 D CASEM

 A TONGE

 B LEAVY

 C WILLIAMS

 J CLAYTON

 J LLOYD

 M GREENFIELD

 S BILYK

 S SEGLETES

 T W BJERKE

 D DANDEKAR

RDRL WMP D

 R DONEY

 C RANDOW
16

INTENTIONALLY LEFT BLANK.