University of Nebraska - Lincoln

DigitalCommons@University of Nebraska - Lincoln

Mechanical & Materials Engineering Faculty Mechanical & Materials Engineering,

Publications Department of

2012

The formulation and computation of the nonlocal J-integral in

bond-based peridynamics
Wenke Hu
University of Nebraska-Lincoln

Youn Doh Ha
Kunsan National University, Korea

Florin Bobaru
University of Nebraska-Lincoln, [email protected]

Stewart A. Silling
Sandia National Laboratories, NM

Follow this and additional works at: https://digitalcommons.unl.edu/mechengfacpub

Hu, Wenke; Ha, Youn Doh; Bobaru, Florin; and Silling, Stewart A., "The formulation and computation of the
nonlocal J-integral in bond-based peridynamics" (2012). Mechanical & Materials Engineering Faculty
Publications. 80.
https://digitalcommons.unl.edu/mechengfacpub/80

This Article is brought to you for free and open access by the Mechanical & Materials Engineering, Department of at
DigitalCommons@University of Nebraska - Lincoln. It has been accepted for inclusion in Mechanical & Materials
Engineering Faculty Publications by an authorized administrator of DigitalCommons@University of Nebraska -
Lincoln.
Int J Fract (2012) 176:195–206
DOI 10.1007/s10704-012-9745-8 This article is a U.S. government work, and is not subject to copyright in the United States.

ORIGINAL PAPER

The formulation and computation of the nonlocal J-integral

in bond-based peridynamics
Wenke Hu · Youn Doh Ha · Florin Bobaru ·
Stewart A. Silling

Received: 31 January 2012 / Accepted: 14 June 2012 / Published online: 11 July 2012
© Springer Science+Business Media B.V. 2012

Abstract This work presents a rigorous derivation Keywords J-integral · Peridynamics · Nonlocal
for the formulation of the J-integral in bond-based peri- methods · Path-independence · Fracture
dynamics using the crack infinitesimal virtual exten-
sion approach. We give a detailed description of an 1 Introduction
algorithm for computing this nonlocal version of the
J-integral. We present convergence studies (m-conver- The J-integral formulation for a linear elastic body in
gence and δ-convergence) for two different geome- 2D was introduced by Rice (1968). The J-integral value
tries: a single edge-notch configuration and a double was shown to be path independent and equivalent with
edge-notch sample. We compare the results with results the energy release rate. The J-integral can be under-
based on the classical J-integral and obtained from stood both as a fracture energy parameter and as a stress
FEM calculations that employ special elements near intensity parameter because J uniquely characterizes
the crack tip. We identify the size of the nonlocal region crack-tip stresses and strains (Hutchinson 1968; Rice
for which the peridynamic J-integral value is near the and Rosengren 1968). The J-integral has been exten-
classical FEM solutions. We discuss how the boundary sively used to compute energy flow to the crack tip, to
conditions and the peridynamic “skin effect” may influ- estimate crack opening and as part of failure criteria for
ence the peridynamic J-integral value. We also observe, ductile materials. Recently, the peridynamic theory was
computationally, the path-independence of the peridy- introduced in order to handle problems with discontin-
namic J-integral. uous fields (Silling 2000). Peridynamics is a nonlocal
formulation of classical mechanics which allows for a
natural treatment of crack initiation, growth and prop-
agation. The peridynamic theory has been successfully
W. Hu · Y. D. Ha · F. Bobaru (B)
Department of Mechanical & Materials Engineering, applied to damage analysis of viscoplastic materials
University of Nebraska-Lincoln, Lincoln, NE 68588, USA (Foster et al. 2010, 2011) dynamic fracture and crack
e-mail: [email protected] branching in glass (Ha and Bobaru 2010, 2011), dam-
age in composite materials from impact or shock load-
S. A. Silling
Sandia National Laboratories, Multiphysics Simulation ing (Xu et al. 2008; Kilic et al. 2009; Hu et al. 2011,
Technology Department, Albuquerque, NM 87185, USA 2012) and nano-scale structures (Silling and Bobaru
2005; Bobaru 2007). The main motivation for study-
Present Address:
ing the J-integral in peridynamics is that it allows com-
Y. D. Ha
Department of Naval Architecture, Kunsan National University, putation of the energy consumed by a growing crack
Kunsan, Korea for any bond failure criterion, and for any dissipative

123
196 W. Hu et al.


mechanisms and nonlinearities that may be occurring, ρ ü(x, t) = f(u(x̂, t)−u(x, t), x̂−x)d Vx̂ +b(x, t)
provided these occur sufficiently close to the crack tip.
H
Also, the J-integral itself could be used as a failure cri-
(1)
terion in computations: it is possible to use a critical
value of J as a criterion for breaking bonds. where f is the pairwise force function in the peridynam-
While for problems in elastic domains without ic bond that connects point x̂ to x, and u is the displace-
cracks, peridynamics converges to classical elasticity ment vector field. ρ is the density and b is the body force
in the limit of the nonlocal region (the horizon) going to intensity. The integral is defined over a region H called
zero, it is important to determine the size of this region the “horizon region”. The horizon region is the com-
for which the peridynamic J-integral results are close pact supported domain of the pairwise force function
to those obtained from the classical J-integral. Silling around a point x, and is taken most often as a disk in
and Lehoucq (2010) presented a state-based peridy- 2D applications (and sphere in 3D). We will abuse the
namic J-integral formulation based on energy balance language and call the radius δ of such a disk the hori-
approach and related it to the Griffith criterion. In the zon. Then, the meaning of the word “horizon” should
present study, we derive the nonlocal J-integral for be clear from the context. The interpretation, selection,
the bond-based peridynamic formulation by means of and use of the peridynamic horizon and its relation to
an infinitesimal virtual crack extension. We describe dynamic crack propagation is discussed in Bobaru and
an algorithm for computing this nonlocal version of Hu (2012).
the J-integral and perform several convergence stud- A micro-elastic material (Silling 2000) is defined as
ies in order to find the size of the nonlocal region (the one for which the pairwise force derives from a poten-
peridynamic horizon size) for which the nonlocal J- tial ω:
integral value is nearing the classical J-integral value. ∂ω(η, ξ )
f(η, ξ ) = (2)
We also investigate the dependence of the peridynamic ∂η
J-integral value on: (1) the boundary conditions; (2) the
where ξ = x̂ − x is the relative position in the reference
peridynamic “skin effect”, and (3) the location of the
configuration and η = û − u is the relative displace-
integral contour (a path-independence study). We com-
ment.
pare our results with results obtained from FEM solu-
A linear micro-elastic potential, which leads to a lin-
tions of the classical Linear Fracture Mechanics using
ear relationship between the bond force and the relative
special elements. The special elements use the exact
elongation of the bond, is obtained if we take
analytical solution of linear elastic fracture mechanics.
The paper is organized as follows: in Sect. 2, we c(ξ )s 2 ξ 
ω(η, ξ ) = (3)
review the basic formulation for bond-based peridy- 2
namics; in Sect. 3 we give the derivation for the peridy- where s the bond relative elongation
namic J-integral based on the infinitesimal virtual crack
ξ + η − ξ 
extension; in Sect. 4 we present the algorithm used for s= (4)
computing the peridynamic J-integral; in Sect. 5, we ξ 
perform two types of convergence (m-convergence and The corresponding pairwise force becomes
δ-convergence) studies for a thin plate with single edge ∂ω(η, ξ ) ∂ ξ + η
notch and double edge notch and we compare the peri- f(η, ξ ) = = c(ξ )s (5)
∂ η ∂η
dynamic results with those obtained from finite element
with
calculations for the classical J-integral using Abaqus
6.10. Calculations related to the path-independence of ∂ ξ + η
=e
the nonlocal J-integral are also shown. Conclusions are ∂η
given in Sect. 6. where e is the unit vector along the direction of the
deformed bond, ξ + η, between x̂ and x in the deformed
configuration.
2 Introduction to the peridynamic theory The function c(||ξ ||) is called micromodulus func-
tion and it represents the bond elastic stiffness. Some
The peridynamic equations of motion are given as possible choices for the micromodulus function are

123
The formulation and computation of the nonlocal J-integral 197

Fig. 1 The skin effect in peridynamics for points close to the

boundary

Fig. 2 The conical micromodulus function

given in 1D (Bobaru et al. 2009), 2D (Ha and Bobaru
2010) and 3D (Silling and Askari 2005). These are
obtained by matching [see Eq. (6)] the elastic strain
energy density at a material point in the bulk (see point
A in Fig. 1) of a microelastic peridynamic material
to the elastic strain energy density from the classical
theory, when both materials are under the same homo-
geneous deformation.

1 c(ξ )s 2 ξ  Fig. 3 Two-dimensional body containing a straight crack
Wclassical = dAx̂ (6)
2 2
H
For simplicity, we use the form of the micromodulus where E and ν are Young’s modulus and Poisson’s
function obtained for a point in the bulk for all points, ratio, respectively. δ is horizon size.
including those that are within δ from the surface (like
point B in Fig. 1). This leads to an effectively softer mate-
rial close to the boundary, since now the integral in Eq. 3 Derivation of the peridynamic J-integral
(6) is over a smaller region and for the match to hold one
would have to increase the micromodulus value. There- In this section we derive the peridynamic J-integral for
fore, the strains will be larger for the “skin” of the domain a mode I crack that grows in a self similar manner in
than in the bulk of the material, for a deformation that a 2D microelastic peridynamic body. We consider a
classically would be homogeneous (see Bobaru and Ha certain domain , which contains a straight crack of
2011). We call this behavior in peridynamics the “skin length a, in equilibrium and zero body force density.
effect”. The skin effect can be compensated for by apply- The peridynamic equations are:

ing a correction factor to the micromodulus of bonds
near a surface (Macek and Silling 2007; Kilic 2008). A f(u(x̂, t) − u(x, t), x̂ − x)dAx̂ = 0 (8)
similar observation can be made about the effect of sur- 
faces on the critical relative elongation parameter that is where Ax̂ represents the nodal area in 2D instead nodal
correlated with the fracture energy of a material, as dis- volume Vx̂ in 3D.
cussed in Ha and Bobaru (2011). The elastic energy density at any x in  for the
In this paper, we use the “conical” micromodulus microelastic material is
functions given below (see Fig. 2) since it gives better 
1
convergence rates to the classical (local) solutions in W (x; a) = ω(η, ξ )dAx̂ (9)
2
elasticity problems compared to the simpler, constant (a)
micromodulus (see Bobaru et al. 2009; Ha and Bobaru where a is the crack length (see Fig. 3).
2010). Consider a fixed coordinate system (X1 , X2 ) as in
   
ξ  24E ξ  Fig. 3. The local coordinate system (x1 , x2 ) is attached
)
c (ξ = C1 1 − = 1−
δ π δ 3 (1 − ν) δ to the crack tip and moves as the crack grows. Assume
(7) that the crack grows by an infinitesimal amount in the x1

123
198 W. Hu et al.

direction. Considering the change in total elastic strain Defining the domain Q(a) = (a)\R(a), the right hand
energy resulting from an infinitesimal virtual extension side of Eq. (13) can be rewritten as
 
of the crack we get: 1   ∂ û
f(η, ξ )dAx dAx̂
 2 ∂a
dU d  (a) R(a)
=− W (x; a)dAx (10)  
da da 1   ∂ û
= f(η, ξ )dAx dAx̂
R(a) 2 ∂a
Q(a) R(a)
 
The right hand side of Eq. (10), by adding and subtract- 1   ∂ û
+ f(η, ξ )dAx dAx̂ (14)
ing R(a) W (x; a + a)dAx , can be written as: 2 ∂a
R(a) R(a)
 Also, from Eq. (8), we have
d  1  
W (x; a + a)dAx = lim
da a→0 a  f(η, ξ )dAx̂ = −  f(η, ξ )dAx̂ (15)
R(a)
⎡ R(a) Q(a)

⎢  By the change of variables x̂ → x and using the linear
×⎣ W (x; a+ a)dAx
admissibility condition (see Silling 2000)
R(a+ a)−R(a)
⎤ f(−η, −ξ ) = −f(η, ξ ) (16)
 and employing Eqs. (15), (14) can be written as
⎥
+  W (x; a+ a) − W (x; a)dAx ⎦  
1   ∂ û
f(η, ξ )dAx dAx̂

R(a)
 2 ∂a
∂ W (x; a) (a) R(a)
=  W (x; a)n1 d S +  dAx (11)  
∂a 1   ∂ û
∂R R(a) = f(η, ξ )dAx dAx̂
2 ∂a
Q(a) R(a)
 
Consider now the second term on the right hand side 1 ∂u
+   f(−η, −ξ )dAx̂ dAx
of Eq. (11) and by using Eqs. (9) and (2), we have 2 ∂a
R(a) R(a)
  
1   ∂ û
 ∂ W (x; a) dAx = f(η, ξ )dAx dAx̂
∂a 2 ∂a
Q(a) R(a)
R(a)
     
∂ û ∂u 1 ∂u
=
1  
f(η, ξ ) · − dAx −   f(η, ξ )dAx̂ dAx
2 ∂a ∂a 2 ∂a
R(a) (a) R(a) R(a)
   
∂ û 1 ∂ û
1
=   f(η, ξ ) · dAx̂ dAx =   f(η, ξ ). dAx̂ dAx
2 ∂a 2 ∂a
R(a) (a) R(a) Q(a)
   
∂u 1 ∂u
−
1  
f(η, ξ )dAx̂ dAx̂ dAx (12) +   f(η, ξ ). dAx̂ dAx
2 ∂a 2 ∂a
R(a) (a) R(a) Q(a)
(17)
The second term on the right hand side of Eq. (12) is Substituting Eq. (17) into Eq. (12), we get

 ∂ W (x; a) dAx
equal to zero because of Eq. (8). Thus, Eq. (12) reduces
to ∂a
R(a)
      
 ∂ W (x; a) dAx = 1   f(η, ξ )dAx ∂ û dAx̂ 1   ∂ û ∂u
= f(η, ξ ). + dAx̂ dAx
∂a 2 ∂a 2 ∂a ∂a
R(a) (a) R(a) R(a) (a)\R(a)
(13) (18)

123
The formulation and computation of the nonlocal J-integral 199


Jperi =  W (x; a)n1 d S
∂R
   
1  ∂ û ∂u
− η
f( , ξ )· + dAx̂ dAx
2 ∂x1 ∂x1
R2 R3

Clearly, when the horizon δ goes to zero, the domain

integral becomes a contour integral and the regions R2
Fig. 4 Integration domain for the peridynamic J-integral. The and R3 reduce to the contour ∂ R. In this case, the for-
red curve is the contour of integration mulation of peridynamic J integral coincides with the
classical J-integral formula. Indeed, the peridynamic
J-integral formulation when horizon goes to zero is
Since the coordinates (x1 , x2 ) are attached at the crack  
∂u
tip, and are defined by Jperi = W (x; a)n1 d S − 
 · τ (x, n)d S
∂x1
x1 = X1 − a, x2 = X2 ∂R ∂R
where τ (x,n) is the force flux vector at any x in the
we have that ∂∂xa1 = −1 direction of unit vector n normal to the tangent of the
By using the chain rule, we obtain contour (see Lehoucq and Silling 2008).
∂ û ∂ û ∂u ∂u
=− , =− (19) Remark Notice that for a fixed horizon size, the con-
∂a ∂x1 ∂a ∂x1 tour needs to be large enough so that the regions R2
Then, we place Eq. (19), (18), and (11) back into Eq. and R3 exist as shown in Fig. 4. Examples of feasible
(10) to get the peridynamic J-integral formula as fol- and unfeasible contours are shown in Sect. 5.
low:
  
 1  
Jperi = W (x; a)n1 d S − f(η, ξ ) 4 An algorithm for calculating the peridynamic
2 J-integral
∂R R(a) (a)\R(a)
 
∂ û ∂u
· + dAx̂ dAx (20) We use a uniform discretization with grid spacing x
∂x1 ∂x1 over the domain to compute the peridynamic J-inte-
where W (x; a) is the strain energy density ∂ R is the gral. Non-uniform discretizations (and even variable
integral contour. R(a) is the region inside the integral horizon size) are possible in peridynamics, see Bobaru
contour, and (a)\R(a) is the region outside the inte- and Ha 2011, but for simplicity here we use uniform
gral contour (see Fig. 3). Silling and Lehoucq (2010) grids. Using the mid-point integration scheme for the
obtained the state-based peridynamic J-integral formu- domain integral in Eq. (20) and the trapezoidal rule for
lation based on an energy approach. We observe that the contour integral in Eq. (20) for the case when the
in the particular case of the bond-based theory, the for- contour is piecewise linear, such as the rectangular box
mula in Silling and Lehoucq (2010) coincides with the in Fig. 6, leads to the following approximation in 2D:
one we obtained here in Eq. (20). n boundary
Jperi ≈ W n1 x
n=1
Discussion: The first term of peridynamic J-integral   
n inner n outer
[see Eq. (20)] is a contour integral, along ∂ R, where 1 j ∂u 1 ∂ û 1
j
− f (u 1 −u i1 , x1 −x1i )
+
R = R1 UR2 (see red curve in Fig. 4). The set R2 is a 2 ∂ x1 ∂ x1
i=1 j=1
“band” of thickness δ inside the contour ∂ R, and R3 is  
∂u 2 ∂ û 2
a band of thickness δ outside of the contour ∂ R. The j j
+ f (u 2 −u i2 , x2 −x2i ) + A j Ai (21)
second term of the peridynamic J-integral is a double ∂ x1 ∂ x1
domain integral. This domain integral is zero unless the where n contour and n inner are the number of nodes along
points in R are in R2 , and the points in \R are in R3 . the integral contour boundary ∂ R and the number of
Thus, Eq. (20) can be written as follows nodes in the inner region (the R2 region in Fig. 4),

123
200 W. Hu et al.

(a) m-convergence (b) δ-convergence

Fig. 5 Graphical description for the m-convergence and δ—convergence in peridynamics. m = δ/ x

respectively. n outer is the number of nodes in the outer peridynamic nonlocal solutions to the classical solu-
region (the R3 region in Fig. 4). Ai is the nodal area tion is expected in the limit of the horizon going to
or node i. The central difference scheme is used for zero. In this paper we perform convergence studies to
∂u/∂x1 and ∂ û/∂x1 as follows investigate for what horizon size, relative to the size of
the sample, and which m-values does the peridynam-
∂u 1 (x1i , x2i ) u 1 (x1i + x, x2i ) − u 1 (x1i − x, x2i )
≈ , ic J-integral get close (with a relative difference of a
∂ x1i 2 x few percentages) to the classical value of the J-inte-
∂u 2 (x1i , x2i ) u 2 (x1i + x, x2i ) − u 2 (x1i − x, x2i ) gral. We then analyze the influence of boundary con-
≈
∂ x1i 2 x ditions (when symmetry conditions are used) and of
j j j j j j the peridynamic “skin effect” (see Sect. 2) on the non-
∂ û 1 (x1 , x2 ) û 1 (x1 + x, x2 ) − û 1 (x1 − x, x2 )
j
≈ , local J-integral. We compare the peridynamic results
∂ x1 2 x with the classical J-integral value as approximately
j
∂ û 2 (x1 , x2 )
j j
û 2 (x1 +
j j
x, x2 ) − û 2 (x1 −
j
x, x2 ) given by Finite Element results using Abaqus with spe-
j
≈ cial crack-tip elements for a plate with an edge notch
∂ x1 2 x
under tension. In order to investigate how symmetric
We compute the peridynamic J-integral with the fol- boundary conditions influence the peridynamic J-inte-
lowing algorithm based on Eq. (21). gral value we analyze a double-edge notch plate under
uniform tension, for which we use symmetry bound-
ary conditions. We also perform calculations that show
5 Numerical examples: convergence studies, the path-independence of the peridynamic J-integral
path-independence, and effects (Table 1).
from the boundaries In following examples, we use the same material
parameters: Young’s modulus is 72 GPa and Poisson’s
In peridynamics, two types of convergence studies ratio is 1/3. Along top and bottom boundaries of the two
are typically used: the δ-convergence (fix the num- different samples, a uniform tensile stress σ = 1 MPa
ber of nodes covered by a horizon, which is propor- is applied. Uniform discretization is used for all com-
tional to m = δ / x, and decrease the horizon size), putations.
and m-convergence (keep the horizon size fixed and In order to test the correctness of the implementa-
increase m) (See Fig. 5 and Bobaru et al. 2009). It is tion, we compute the peridynamic J-integral for a plate
important to note that in m-convergence studies, the (see Fig. 6) without a notch, using the integral contour
peridynamic approximate solutions are not supposed in Fig. 6. In this case, peridynamic J-integral values are
to converge to the classical solution, but to the non- in the range of 10−15 for any m-values and horizon size
local solution for the particular horizon size for which we tried. For instance, the peridynamic J-integral value
the m-convergence study is executed. In this instance, is 3.3 × 10−15 Pa m for m = 9 and δ = 1.5 mm. The
comparing the results with the classical solution is only J-integral on a closed contour in an elastic material is
made for illustration purposes. Convergence of the zero.

123
The formulation and computation of the nonlocal J-integral 201

Table 1 Algorithm to compute the peridynamic J-integral

1: Get the nodal displacements and elastic strain energy density
at all nodes
2: Define the geometry of the integral contour
3: Compute the first term (contour integral) of the peridynamic
J-integral
(a) Find nodes along the contour ∂ R and the corresponding
outer normal vector to the integral contour at such nodes
(b) Evaluate the line integral using trapezoidal rule
4: Compute the second term (domain double integral) of the
peridynamic J integral
(a) Define the inner region (R2 ): one horizon size inside of the
integral contour ∂ R
(b) Define the outer region (R3 ): one horizon size outside of
the integral contour ∂ R
(c) Search all viable source nodes xi in inner region
(d) Search all viable nodes xj in outer region Fig. 7 A sample Abaqus mesh with special elements around the
(e) Compute the second term based on the discretization in Eq. crack tip, used to compute the classical J-integral value. The dark
(21) solid line is the crack
5: J peri = First_term–Second_term

gives us a converged value for the classical J-integral

of about 19.76 Pa m.
Since the J-integral involves the components of the
displacement field, it is of interest to compare the dis-
placements obtained from peridynamics with a “large”
and a “small” horizon size with the classical displace-
ments obtained from a sufficiently fine FEM mesh. As
shown in Figs. 8 and 9, the x- and y-displacements
from peridynamic computation are almost identical to
the FEM calculations when the horizon size, relative to
the sample size and crack size is below 1/20. In these
figures the same legend is used for both the peridynamic
results and the finite element results. The legend levels
are produced automatically by Abaqus, and we set the
Fig. 6 Geometry configuration for single edge-notched plate same values in Tecplot to post-process the peridynamic
(red square is the J-integral contour) results.
Two types of convergence studies (m-convergence
and δ-convergence) are shown in Table 2 and Fig. 10.
5.1 Single-notched specimen For a fixed m, the peridynamic J-integral approaches
the classical solution when the horizon size decreases,
Consider a single edge-notch square plate with dimen- so we do observe δ-convergence to the classical J-inte-
sion of 10 cm × 10 cm with a length of the notch of gral value. We also note that the m-convergence curves
5 cm, as shown in Fig. 6. We use the square contour (for a fixed δ) start to level off when m becomes larger
shown in Fig. 6 to compute the peridynamic J-inte- than 6. Increasing m for a fixed horizon size delivers
gral and the corresponding classical value with Abaqus. solutions that approach the exact nonlocal solution. The
In Abaqus, we use the “seam crack” option to create results of Fig.10 indicate good candidates for both δ and
the crack and use special crack-tip quadratic elements m to use in order to obtain the peridynamic J-integral
around the crack tip as shown in Fig. 7. In the actual values within a few percentages from the FEM results.
FEM computation, the total number of nodes is about We stress that, in general, the FEM results may also
30,000 (many more nodes than shown in Fig. 7), which be a few percentages away from an analytical solution

123
202 W. Hu et al.

Fig. 8 Comparison of vertical displacements obtained with: a δ = 24 mm, m = 3 (about 144 nodes) (geometry and loading as
the FEM (about 30,000 nodes) , b peridynamics using δ = given in Fig. 6). The legend is in meters
1.5 mm, m = 9 (about 360,000 nodes), and c peridynamics using

Fig. 9 Comparison for horizontal displacements obtained with: a the FEM, b peridynamics using δ = 1.5 mm, m = 9, and c peridy-
namics using δ = 24 mm, m = 3 (geometry and loading as given in Fig. 6). The legend is in meters

of the classical model. Thus, when the intention is to 5.2 Double-notched specimen
obtain peridynamic J-integral values close to the clas-
sical ones, an m of about 6 and a horizon size smaller In this section, we consider a double-notch plate with
than 6 mm (a ratio of at least 1/10 to the crack length, dimensions 20 cm × 10 cm. The length of each notch
and 1/20 to the sample dimensions) are good choices. is 5 cm. By making use of symmetry, we can reduce
the problem to analyzing the same configuration as
before (see Fig. 6) except that now we impose sym-
Remark the peridynamic “skin effect” mentioned in metry conditions on the displacements along the right-
Sect. 2 exists on the crack surfaces. Hence, we expect hand boundary (see Fig. 11). We use the same integral
peridynamic J-integral results to be higher than the contour as in Sect. 5.1 and three different horizon sizes
FEM results since we effectively have a softer mate- δ = 6 mm, δ = 3 mm, and δ = 1.5 mm and m = 3, 6,
rial (Young’s modulus value half of that in the bulk) and 9 to observe how convergence is influenced by the
around the crack tip than the material in the bulk. We presence of the symmetry boundary condition. While
assigned a softer material to a thin region, of thickness in a classical model discretized with finite elements
equal to 1.5 mm, the same as the size of the smallest the symmetry boundary condition is simply imposed
peridynamic horizon used here, around the crack line on the boundary nodes where zero horizontal displace-
in the FEM model and the J-integral value from the ments are enforced, in a nonlocal peridynamic model
Abaqus calculation increased by about 2 %. this type of condition is more delicate. In principle, in

123
The formulation and computation of the nonlocal J-integral 203

Table 2 Peridynamic
δ = 24 mm δ = 12 mm δ = 6 mm δ = 3 mm δ = 1.5 mm
J-integral values for the
single edge notched m=3 30.15 Pa m 23.64 Pa m 21.74 Pa m 20.83 Pa m 20.41 Pa m
specimen
m=6 26.52 Pa m 23.05 Pa m 21.54 Pa m 20.85 Pa m 20.51 Pa m
m=9 26.17 Pa m 22.95 Pa m 21.53 Pa m 20.86 Pa m 20.53 Pa m

Fig. 11 The double edge notched specimen with symmetric

Fig. 10 The relative difference between the peridynamic results boundary conditions and the J-integral contour
for the single-notch sample and the classical J-integral obtained
with the FEM (Abaqus)

the right-hand boundary, where the symmetric bound-

the discrete peridynamic model we should enforce the ary condition is applied (see Fig. 12). This influences,
boundary nodes to have zero horizontal displacements, to a certain extent, the results for the J-integral, as we
but in addition to that, for a sliver of thickness δ on will see below. We also notice that the difference in
the corresponding boundary we should also prescribe brightness between the FEM and peridynamics results
the horizontal displacement field (see Silling 2000 for a are not due to differences in the values, but are caused
discussion on nonlocal Dirichlet boundary conditions). by differences in the colors used by Abaqus (for the
The problem is that those values are not known. Ways FEM results) and Tecplot (which we used to post-pro-
around this are possible: for example, one could set zero cess the peridynamic results).
horizontal displacements for all nodes in the region of The m-convergence trends are the same for all three
thickness δ, and in the limit of the horizon going to cases as shown in Fig. 13. We stress again that in
zero this willapproach the actual symmetry condition m-convergence studies, the peridynamic approximate
desired; alternatively, one could impose zero horizon- solutions are not supposed to converge to the classical
tal displacements only on the boundary nodes, just like solution, but to the nonlocal solution for the particu-
in a finite element model. We choose the latter strategy lar horizon size for which the m-convergence study is
because of simplicity and similarity to the way bound- performed. Since we do not have the exact solution for
ary conditions are imposed for the classical model. the nonlocal problem, we compare the results with the
The finite element solution from Abaqus with classical solution for illustration purposes. See also the
special crack-tip elements is 12.90 Pa m. comments in Bobaru and Duangpanya (2010, 2012),
As shown in Fig. 12, the strain energy density on a peridynamic model for diffusion processes. The
obtained with peridynamics for a sufficiently small data in Table 3 shows that the peridynamic solutions are
horizon size is very close to that given by a converged all larger than the FEM result, which is consistent with
FEM solution. However, some high strain energy den- the results obtained in section 5.1. The effectively softer
sity values are observed from the peridynamic result on material around the crack tip (induced by the peridy-

123
204 W. Hu et al.

Fig. 12 Strain energy density results with: a FEM (30,000 nodes), b peridynamics with δ = 1.5 mm, m = 9 (about 360,000 nodes)
c peridynamics with δ = 6 mm, m = 3 (about 2700 nodes)

Table 3 Convergence study with double edge notched specimen

δ = 6 mm δ = 3 mm δ = 1.5 mm

m=3 14.83 Pa m 14.01 Pa m 13.58 Pa m

m=6 15.49 Pa m 14.50 Pa m 13.99 Pa m
m=9 15.90 Pa m 14.77 Pa m 14.17 Pa m

section by choosing different integration paths around

the crack tip.

5.3 The path-independence of the peridynamic

J-integral

As the peridynamic J-integral over a closed contour in

Fig. 13 The relative difference between peridynamic results for a micro-elastic material is zero, analytically, the peri-
the double edge-notch (with symmetric boundary conditions) dynamic J-integral is path independent. We test here
and the classical J-integral value obtained from a FEM (Abaqus) different contours to observe if, and by how much, the
analysis value of the peridynamic J-integral (computed using the
algorithm in Table 1) changes with changing the inte-
gration contour. For these tests we use the smallest hori-
namic “skin effect”) is responsible for this. In addition, zon and finest grid from the computations above (m =
for a given m-value, the smaller the horizon size the 9 and δ = 1.5 mm) and perform calculations for the
closer the peridynamic result is to the classical value, single edge-notch problem and the double edge-notch
confirming that δ-convergence takes place. With a suf- problem with symmetry boundary conditions. We use
ficiently small horizon, the relative difference between the three different contours shown in Fig. 14: the clos-
the classical J-integral value and the peridynamic J- est contour to the crack tip that remains feasible from
integral becomes smaller than a few percentages. The the point of view of evaluating the integrals in Eq. (20)
differences between the peridynamic results and the (see also Fig. 4), another away from the crack tip and
classical J-integral value are, however, larger now than the boundaries, and finally one near the boundaries, but
before, for a similar horizon size. The likely motive for at least one horizon size away from the boundaries to
this change is the imposition of the symmetric boundary remain feasible. Thus, contour (a) in Fig. 14 is one hori-
conditions. The influence of the imposition of the sym- zon size away from the crack tip, contour (b) is in the
metric boundary conditions on the calculation of the middle of the structure, and contour (c) is one horizon
peridynamic J-integral is further investigated in the next away from the boundary of the sample (Table 4).

123
The formulation and computation of the nonlocal J-integral 205

Fig. 14 Three different

integral contours to
compute peridynamic J
integral (top). Feasible and
unfeasible contours (bottom
row) for the computation of
the nonlocal J-integral (see
the remark at the end of
Sect. 3) near the crack tip
and near a domain boundary

Table 4 Values of the peridynamic J-integral on the three dif- Table 5 Relative difference between the peridynamic J-integral
ferent contours from Fig. 14 on contours (a) and (c) in Fig. 14 and the corresponding values
obtained on contour (b)
Contour (a) Contour (b) Contour (c)
Contour (a) Contour (c)
Single edge-notch 20.11 Pa m 20.54 Pa m 19.84 Pa m
Double 13.88 Pa m 14.17 Pa m 13.37 Pa m Single edge-notch 2.09 % 3.41 %
edge-notch Double edge-notch 2.05 % 5.65 %
(using (using symmetric
symmetric boundary conditions)
b.c.’s)

approximation is supposed to approach the exact solu-

tion of the nonlocal problem and not the exact solution
We measure the relative difference between the val- of the classical problem. It is known that, in certain
ues obtained based on contour (b) and those on contours cases, nonlocal results may be closer to the classical
(a) and (c). As shown in Table 5, the relative differences results for some δ1 and m 1 , and higher for a δ2 < δ1
between contour (a) and contour (b) for both the single and m 2 > m1. This means that we may match the exact
edge-notch and double edge-notch specimens are very classical solution for some finite δ3 in [δ2 , δ1 ] and m 3
similar, equal to about 2 %. However, the relative dif- in [m 1 , m 2 ] (see also Bobaru and Duangpanya 2010,
ference between contour (c) and contour (b) is larger in 2012).
the double-edge case (with symmetric boundary con- In conclusion, a contour for computing the peridy-
ditions) than in the single edge-notch case. The reason namic J-integral needs to avoid the boundary region
is likely the way the nonlocal Dirichlet boundary con- because of the skin effect and imposition of nonlocal
dition is enforced and the peridynamic “skin-effect”. Dirichlet boundary conditions. On the other hand, the
Notice that the results on contour (c) are actually closer contour may be selected near the crack tip, with the only
to the classical results obtained with the FEM, but these condition being the feasibility requirement mentioned
trends should not be misinterpreted. The horizon size, in the remark in Sect. 3 and shown in Fig. 14.
the number of nodes in the horizon, and the type of
boundary conditions, all “influence” the approximate
solutions of nonlocal models, and these results may 6 Conclusions
be lower than a classical value for some m-values,
and higher for some other m-values. We recall that in In this paper we derived the 2D peridynamic J-integral
m-convergence, as m tends to infinity, the numerical formulation based on the infinitesimal virtual extension

123
206 W. Hu et al.

of the mode-I crack. We discussed the computation of Bobaru F, Hu W (2012) The meaning, selection, and use
the peridynamic J-integral value which involves a dou- of the peridynamic horizon and its relation to crack
branching in brittle materials. Int J Fract. doi:10.1007/
ble integral in contrast with the classical J-integral. We
s10704-012-9725-z
showed that the peridynamic J-integral converges to Bobaru F, Yang M, Alves LF, Silling SA, Askari E, Xu J (2009)
the classical J-integral formulation when the horizon Convergence, adaptive refinement, and scaling in 1D peri-
size goes to zero. Using a simple discretization, we dynamics. Int J Numer Meth Eng 77(6):852–877
Foster JT, Silling SA, Chen WW (2010) Viscoplasticity using
performed convergence studies (m-convergence and
peridynamics. Int J Numer Meth Eng 81:1242–1258
δ-convergence) for two specimens: a single edge-notch Foster J, Silling SA, Chen WW (2011) An energy based failure
and a double edge-notch sample, for which symmetry criterion for use with peridynamic states. Int J Mult Comp
boundary conditions were enforced. We observed that Eng 9:675–688
Ha YD, Bobaru F (2010) Studies of dynamic crack propaga-
the peridynamic results approach within a few percent-
tion and crack branching with peridynamics. Int J Fract
ages the classical J-integral values obtained using the 162:229–244
FEM with special element around the crack tip, when Ha YD, Bobaru F (2011) Characteristics of dynamic brittle frac-
the horizon size is less than about 1/20 of the crack size ture captured with peridynamics. Eng Fract Mech 78:1156–
1168
and sample size. In particular, we discussed the influ-
Hu W, Ha YD, Bobaru F (2012) Peridynamic simulations of
ence of the peridynamic “skin-effect” around the crack dynamic fracture in unidirectional fiber-reinforced compos-
tip and along the boundaries, on the value of the peri- ites. Comput Methods Appl Mech Eng. 217–220: 247–261
dynamic J-integral. We computed the peridynamic J- Hu W, Ha YD, Bobaru F (2011) Modeling dynamic fracture and
damage in fiber-reinforced composites with peridynamics.
integral along different paths and very small variations
Int J Mult Comput Eng 9:707–726
were seen (attributable to the approximation error) Hutchinson JW (1968) Singular behavior at the end of a ten-
between the contours, except when the integration path sile crack tip in a hardening material. J Mech Phys Solids
was within a horizon distance from the boundaries. Due 16:13–31
Kilic B (2008) Peridynamic theory for progressive failure pre-
to nonlocality, special care also needs to be paid when
diction in homogeneous and heterogeneous materials. Ph.D.
symmetric boundary conditions are imposed. thesis, Department of Aerospace and Mechanical Engineer-
ing, University of Arizona, pp 68
Acknowledgments The WH, YDH, and FB are thankful for Kilic B, Agwai A, Madenci E (2009) Peridynamics theory for
the financial support offered through research contracts between progressive damage prediction in center-cracked compos-
UNL and the ARO (Dr. Larry Russell), ARL (project coordi- ite laminates. Compos Struct 90:141–151
nators Dr. C.F. Yen and Dr. C. Randow), ARO award number Lehoucq RB, Silling SA (2008) Force flux and the peridynamic
58450EG, and Sandia National Laboratories (contract number stress tensor. J Mech Phys Solids 56:1566–1577
568428). Sandia National Laboratories is a multi-program lab- Macek RW, Silling SA (2007) Peridynamics via finite element
oratory operated by Sandia Corporation, a wholly owned sub- analysis. Finite Elem Eng Des 43:1169–1178
sidiary of Lockheed Martin company, for the U.S. Department Rice JR (1968) A path independent integral and the approxi-
of Energy’s National Nuclear Security Administration under con- mate analysis of strain concentration by notches and cracks.
tract DE-AC04-94AL85000. The computations in this work were J Appl Mech 9:379–386
completed utilizing the Holland Computing Center of the Uni- Rice JR, Rosengren GF (1968) Plane strain deformation near a
versity of Nebraska. crack tip in a power-law hardening material. J Mech Phys
Solids 16:1–12
Silling SA (2000) Reformulation of elasticity theory for dis-
continuities and long-range forces. J Mech Phys Solids
References 48:175–209
Silling SA, Askari E (2005) A meshfree method based on
Bobaru F (2007) Influence of van der Waals forces on increasing peridynamic model of solid mechanics. Comput Struct
the strength and toughness in dynamic fracture of nanofiber 83:1526–1535
networks: a peridynamic approach. Model Simul Mater Sci Silling SA, Bobaru F (2005) Peridynamic modeling of mem-
Eng 15:397–417 branes and fibers. Int J Nonlinear Mech 40:395–409
Bobaru F, Duangpanya M (2010) The peridynamic formulation Silling SA, Lehoucq RB (2010) Peridynamic theory of solid
for transient heat conduction. Int J Heat Mass Transfer mechanics. Adv Appl Mech 44:73–168
53:4047–4059 Xu J, Askari E, Weckner O, Silling SA (2008) Peridynam-
Bobaru F, Duangpanya M (2012) A peridynamic formulation for ic analysis of impact damage in composite laminates.
transient heat conduction in bodies with evolving disconti- J Aerosp Eng 21:187–194
nuities. J Comput Phys 231(7):2764–2785
Bobaru F, Ha YD (2011) Adaptive refinement and multiscale
modeling in 2D peridynamics. Int J Mult Comput Eng
9(6):635–660

123