## Abstract

In this paper, we compare small deformations in an infinite linear elastic body due to the presence of point loads within the classical, local formulation to the corresponding deformations in the peridynamic, non-local formulation. Owing to the linearity of the problem, the response to a point load can be used to obtain the response to general body force loading functions by superposition. Using Laplace and Fourier transforms, we thus obtain an integral representation for the three-dimensional peridynamic solution with the help of Green’s functions. We illustrate this new theoretical result by dynamic and static examples in one and three dimensions. In addition to this main result, we also derive the non-local three-dimensional jump conditions, as well as the weak formulation of peridynamics together with the associated finite element discretization.

## 1. Introduction

The prediction of the spontaneous nucleation of cracks as well as the subsequent propagation in load-carrying structures such as the wing of an aeroplane presents a long-standing problem in continuum mechanics of solids. In a complex loading situation such as a bird strike, multiple interacting cracks can be present at the same time. The generally curvilinear path along which a crack propagates in a three-dimensional structure is not known *a priori* and must be determined as part of the solution. In addition, the crack path also depends on the *material*. Recently, anisotropic composite materials (such as carbon fibre-reinforced polymer) are replacing more traditional isotropic materials (such as aluminium) in part because of their higher specific strength and promising significant weight savings. For a given loading, the crack path in a composite structure depends on the details of the underlying microstructure (table 1).

The level of fidelity of simulations using traditional finite element (FE) codes in predicting inelastic material behaviour has lagged behind their capabilities in elastic stress analysis. This difficulty arises in part because the mathematical framework on which these methods are based assumes that the body remains continuous as it deforms. Hence, these methods break down at a crack tip and special techniques must be used, which typically require the path of the crack to be known in advance, among other difficulties.

As an attempt at improving this situation, a theory of continuum mechanics known as peridynamics has been recently proposed by Silling (2000). The objective of peridynamics is to reformulate the basic mathematical description in such a way that the same equations hold both at a crack tip as well as in the far field. In this approach, internal forces are expressed through pairwise interactions, the so-called bonds, between pairs of material points. The finite interaction distance introduces non-locality. The complete constitutive model, including damage, is determined at the bond level. Cracks grow when and where it is energetically favourable for them to do so.

Non-local theories in continuum mechanics have been known since the 1970s from articles by Kröner (1967), Eringen (1972, 1992), Eringen & Edelen (1972), Kunin (1983), Rogula (1982) and co-authors. These theories aim to describe certain effects that are not captured accurately in the corresponding local formulation, e.g. the physically unreasonable infinite stresses found at a crack tip in local linear elasticity. More recently, non-local approaches have been discussed, for example, in Altan (1989, 1991), Bazant & Jirásek (2002), Pisano & Fuschi (2003), Polizzotto (2001), Silling (2000), Wang & Dhaliwal (1993*a*,*b*), Chen *et al.* (2004) and Lei *et al.* (2005). Among these approaches, peridynamics falls into the category of strongly non-local methods owing to the presence of an integral operator as opposed to higher order gradient operator in the equation of motion. Recent theoretical developments in peridynamics can be found in Weckner & Abeyaratne (2005), Emmrich & Weckner (2005), Silling *et al.* (2007), Emmrich & Weckner (2007), Weckner & Emmrich (2007), Lehoucq & Silling (2008), Silling & Lehoucq (2008), Bobaru *et al.* (2009) and Warren *et al.* (2009). In Silling (2003), Bobaru & Silling (2004), Gerstle & Sau (2004), Bobaru *et al.* (2005), Gerstle *et al.* (2005), Silling & Askari (2005), Silling & Bobaru (2005), Askari *et al.* (2006), Dayal & Bhattacharya (2006), Bobaru (2007), Xu *et al.* (2007), Parks *et al.* (2008), Silling *et al.* (2008) and Xu *et al.* (2008) simulations based primarily on a numerical implementation of peridynamics called EMU (Silling 2006) cover a wide range of very interesting problems involving the spontaneous initiation of discontinuities followed by their unguided propagation.

However, the study of strictly elastic problems and their relationship with the classical local formulation has been somewhat neglected. This neglect stems from the following. Local elasticity is well understood and yields satisfactory results for a large class of important problems typically involving the determination of the stress and strain fields. While it is possible to solve the same set of elastic problems using the non-local peridynamic formulation, the computational costs would be considerably higher. However, the fidelity of the classical approach in determining the initiation and propagation of cracks is clearly lagging behind its ability to determine elastic stress and strain fields. Peridynamics has been developed in an effort to close this gap in predicting *inelastic* deformations.

However, in order to increase the fidelity in predicting inelastic material behaviour a better understanding of the *elastic* case is essential, which is the focus of this paper. On the one hand, we study the (expected) convergence of the non-local formulation to local elasticity in the large wavelength limit. On the other hand, several fundamental differences arise in peridynamics, which might help understand phenomena at length scales where local elasticity is not applicable.

The paper is organized as follows. In §2, previous research using Green’s functions for elastic problems in one dimension is summarized. In §3, this approach is extended to the three-dimensional case. This section also contains a careful comparison between local and non-local elasticity, including the previously unpublished derivation of the three-dimensional jump conditions as well as the weak formulation and the corresponding FE discretization. Examples are given in §4. Section 5 concludes and discusses open questions.

## 2. Local and non-local elastic deformations in one dimension

The equation of motion at time *t* for the material point *x* in an infinite, homogeneous body in one spatial dimension is
2.1
where the linear operator acting on the displacement field *u*(*x*,*t*) captures internal forces while *b*(*x*,*t*) captures external forces. *ρ* is the density.1

In local elasticity, the internal forces are represented by the differential operator
2.2
with Young’s modulus *E*. This is the well-known wave equation. On the other hand, the non-local peridynamic formulation for an infinite linear micro-elastic material leads to the integral operator
2.3
with the so-called micromodulus function *c*(*ξ*)=*c*(−*ξ*).2

The associated energy balance is obtained by multiplying the equation of motion (2.1) by the velocity and subsequent integration over the body:3 2.4 Given the initial data , it has been shown in Weckner & Abeyaratne (2005) that the solution of (2.1) has the following integral representation 2.5 where and

## 3. Local and non-local elastic deformations in three dimensions

### (a) Kinematics

The material particles *X* are addressed by their position in the reference configuration at say *t*=0, represented by their position vectors . At time *t*, the particle *X* has moved to its current position ** y**(

**,**

*x**t*)=

**+**

*x***(**

*u***,**

*x**t*), where

**is the displacement field. The velocity of particle**

*u**X*is defined as . The relative position of two particles

*X*and

*X*′ in the reference configuration is denoted by

**=**

*ξ***′−**

*x***and is called a peridynamic bond. The corresponding relative position in the current configuration is given by**

*x***(**

*y***′,**

*x**t*)−

**(**

*y***,**

*x**t*)=

**+**

*ξ***=(**

*η**ξ*+

*η*)

*n*_{ξ+η}with the relative displacement

**=**

*η***(**

*u***′,**

*x**t*)−

**(**

*u***,**

*x**t*) and the unit vector pointing from particle

*X*towards particle

*X*′,

*n*_{ξ+η}. For smooth deformation fields, we can introduce the deformation gradient

**(**

*F***,**

*x**t*)=(∇

**(**

*y***,**

*x**t*))

^{T}:

**+**

*ξ***=**

*η***·**

*F***+**

*ξ**O*(

*ξ*

^{2}) or d

**(**

*y***,**

*x**t*)|

_{t}=

**(**

*F***,**

*x**t*)·

*d*

**. Owing to the balance of linear and angular momentum, the non-local force that particle**

*x**X*′ exerts on

*X*must be a central force,

**f**(

**,**

*x***′,**

*x**t*)=

*f*(

**,**

*x***′,**

*x**t*)

*n*_{ξ+η}. These quantities are illustrated in figure 1.

### (b) Comparison of local and non-local elasticity

#### (i) Equation of motion

The equation of motion at time for the material point in an infinite, linear elastic, isotropic, homogeneous body as formulated within the framework of local continuum mechanics is given by the Navier equations (see Elmore & Heald 1985).
3.1
with the reference density *ρ*, the second Piola–Kirchhoff stress tensor ** S**, Cauchy’s infinitesimal strain tensor

**ϵ**and the external force field

**.**

*b**λ*and

*μ*are the Lamé constants which can alternatively be expressed in terms of the Young’s modulus

*E*=

*μ*(3

*λ*+2

*μ*)/(

*λ*+

*μ*) and the Poisson ratio

*ν*=

*λ*/2(

*λ*+

*μ*). At time

*t*=0, we have the initial conditions 3.2

In the strongly non-local peridynamic formulation of continuum mechanics (see Silling 2000), the equation of motion for an infinite, isotropic, homogeneous, linear microelastic, pairwise equilibrated material takes the form
3.3
and
3.4
where the interaction ‘horizon’ is taken to be the sphere with centre ** x** and radius , which is the collection of all the bonds

**with length**

*ξ**ξ*=∥

**∥<**

*ξ**δ*. The symmetric micromodulus tensor

**(**

*C***)=**

*ξ***(−**

*C***)=**

*ξ*

*C*^{T}(

**), more precisely, the micromodulus function**

*ξ**Λ*(

*ξ*), contains all constitutive information.

Formally, the differential operator of local elasticity has been replaced by the non-local peridynamic integral operator: 3.5 3.6 and 3.7

#### (ii) Jump conditions

One first important difference between the local and non-local formulation is the jump condition for linear momentum formulated across a moving discontinuity surface . The spatial (immaterial) point momentarily occupies the material point , so . At , has the surface normal and velocity4. At the fixed time *t*, the linear momentum jump condition reads as5
3.8
where is the mass per unit volume in the actual configuration which is related to the reference density by . The jump conditions for mass and continuity of displacement are identical in both formulations
3.9
and
3.10
Substituting equation (3.9) into the non-local equation (3.8) and assuming that the deformation is such that the density always remains positive (det(** F**)>0), it follows that
3.11
If we choose an arbitrary, immaterial discontinuity surface with , it follows from equation (3.11) that we cannot have a jump in the velocity field. In a continuous deformation, equation (3.10) then implies that the deformation gradient tensor, and therefore the strain tensor, must be continuous as well. This is the trivial case where all the fields are smooth.

On the other hand, if we assume an (initial) velocity field with a jump discontinuity as in the case of the well-known Riemann problem, equations (3.9) and (3.11)) imply or with . This means that in a non-local deformation, which respects both balance of mass and linear momentum, the normal component of the velocity field is always continuous and any discontinuity surface moves like a material interface in the normal direction. The only possible velocity jumps lie in the tangent plane that locally coincides with the generally curved discontinuity surface: with the projector . As , it follows that . For arbitrary and det(*F*^{+})>0, this implies . However, according to equation (3.10) this is no longer compatible with a discontinuous velocity field.

In summary, a jump in the velocity field can only occur in the components that lie in the plane tangent to the discontinuity surface . It is always accompanied by a jump in displacement field, the location of which is fixed at the Lagrangian point .6 This is an important difference between the non-local and local formulation, the latter allowing for shock-waves in which both velocity and strain field can simultaneously be discontinuous.

#### (iii) Energy balance

The energy balance associated with the equation of motion (3.5) is obtained by multiplication with the velocity field ** v**(

**,**

*x**t*) and integration over the domain 3.12 The definition of the kinetic energy and the power of the external forces is the same in both local and non-local elasticity while the elastic energy differs: the non-local elastic energy density is given by the integration of the pairwise potential

*w*(

**,**

*ξ***):∂**

*η*_{η}

*w*(

**,**

*ξ***)=**

*η***(**

*C***)·**

*ξ***=**

*η***(**

*f***,**

*x***′,**

*x**t*) over the horizon. For example, in a homogeneous deformation characterized by a constant deformation gradient

*F*_{0}, we have and with the first and second invariants of the strain tensor and

*I*

_{2}=ϵ

_{I}ϵ

_{II}+ϵ

_{II}ϵ

_{III}+ϵ

_{III}ϵ

_{I}(see Weckner

*et al.*2007). Requiring that the elastic energy density for any homogeneous deformation is identical in both local and non-local elasticity leads to the restriction7 3.13 It further relates the micromodulus function

*Λ*to the Lamé constant

*λ*3.14

#### (iv) Equation of motion in Fourier space

Applying the Fourier transform with respect to the spatial coordinate ** x**, we can equivalently characterize the equation of motion (3.5) by the acoustic tensor

**(**

*M***) as follows:8 3.15 3.16 3.17 3.18 3.19 3.20 3.21 and 3.22 The transformed initial conditions are 3.23 Using equation (3.14) together with equations (3.21) and (3.22), we see that the first non-zero term in the Taylor expansion of the non-local acoustic tensor**

*k*^{NL}

**(**

*M***) in the large wavelength limit**

*k**k*→0 coincides with the local acoustic tensor

^{L}

**(**

*M***) for materials with**

*k**λ*=

*μ*. Alternatively, the convergence of the non-local peridynamic equation (3.3) towards the local Navier equations (3.1) can be shown directly in (

**,**

*x**t*) space (see Emmrich & Weckner 2007 (linear bond-based formulation); and Silling & Lehoucq 2008 (non-linear state-based formulation)).

#### (v) Wave propagation

A physical interpretation of the acoustic tensor can be obtained by studying the propagation of plane waves , where *ω*(*k*) is the dispersion relation relating the angular frequency *ω* to the wavenumber *k*=∥** k**∥. Substituting the wave ansatz into the equation of motion (3.5) leads to the eigenvalue problem . From this, we can identify pressure and shear waves where and ,

*t*_{k}

**·**

*n*_{k}=0 travelling with the phase velocities 3.24 and 3.25 In local elasticity, the phase velocity does not depend on the wavelength for either pressure or shear waves: , (see equations (3.17) and (3.18)). In contrast, peridynamics always leads to wave dispersion. Note that wave dispersion is present in most augmented models of continuum mechanics such as the weakly non-local higher order gradient theories (see Mindlin 1965) or the strongly non-local Eringen-type models (see Eringen 1972). These models aim to describe certain effects that are not captured accurately in local linear elasticity. One such example is the non-linearity found in experimentally measured dispersion relations, reflecting the inability of real materials to sustain waves of arbitrarily small wavelength as described in Graff (1991). In this context, it remains an important open question whether it is possible to determine the micromodulus function

*Λ*(

*ξ*) from experimentally obtained dispersion data.

As an example, consider the case where all points interact: . The micromodulus function is assumed to be either exponential or trigonometric , where the length scale determines the degree of non-locality. As shown below, the exponential form of the micromodulus function leads to wave dispersion for any finite wavelength *χ*=2*π*/*k* while the trigonometric micromodulus function behaves like a low-pass filter: waves with a wavelength larger than *χ*_{c}=2*π**δ* travel with the same phase velocity as in the local formulation, waves smaller than the cut-off wavelength *χ*_{c} are dispersed. The exponential form of the micromodulus function results in the following phase velocities as a function of the normalized wave number :
3.26
and
3.27
while the trigonometric micromodulus leads to
3.28
and
3.29
The large wavelengths expansion confirms the convergence towards local elasticity for materials with *λ*=*μ* as illustrated in figure 2.

Furthermore, one can see that the components of the acoustic tensor become independent of the wave number in the small wavelength limit 3.30 so the phase velocity goes to zero as . In the example above, and .

#### (vi) Weak formulation

Another way to characterize the motion ** u** is given by the variational problem
3.31
where the Lagrangian density is given by

*L*(

**,**

*x**t*)=

*e*

_{kin}(

**,**

*x**t*)−

*e*

_{el}(

**,**

*x**t*)−

*e*

_{b}(

**,**

*x**t*) and the potential of the external force field is defined as

*e*

_{b}(

**,**

*x**t*)=−

**(**

*b***,**

*x**t*)·

**(**

*u***,**

*x**t*). The Euler–Lagrange equation associated with the variational problem (3.31) is the equation of motion (3.5) (see Weckner & Emmrich 2007). The only difference between the local and the non-local formulation is again the elastic energy density.

#### (vii) Finite-element discretization

Finally, one can also characterize the local and non-local formulations by their corresponding stiffness matrices. Introducing the Ritz ansatz
3.32
into the variational problem (3.31) leads to the discretized equation of motion for the 3*N* unknown displacements
3.33
The mass matrix and the inhomogeneity are identical in local and non-local elasticity
3.34
and
3.35
while the local and non-local stiffness matrices differ
3.36
3.37
with
where the limits of integration in the inner integral have been formally extended from to the whole domain since by definition ** C**(

**)=**

*ξ***0**∀∥

**∥≥**

*ξ**δ*. The bandwidth of the stiffness matrix of local elasticity depends on the support of the basis function

*g*

_{α}(

**). For computational efficiency, one typically introduces a numerical length scale**

*x**δ*

_{n}such that

*g*

_{α}(

**)=0 ∀**

*x**x*≥

*δ*

_{n}, ∀

*α*. Within the local formulation, this leads to a sparse stiffness matrix with a band structure. Within the non-local formulation, the situation is more complex as there are two length scales present: the numerical length scale

*δ*

_{n}and the peridynamic horizon

*δ*. To simplify the discussion, we assume that . Then the resulting non-local stiffness matrix will be dense, representing the physical interaction between any two particles in the body. Returning to the general case , the non-local stiffness matrix will generally have a larger bandwidth than the local stiffness matrix. The higher computational costs are justified, for example, in the presence of propagating cracks (corresponding to discontinuous displacements) since this class of solutions are not contained in the Sobolev space associated with the weak formulation of the Navier equations.

Using the concept of Dirac distributions Δ(*x*)
3.38
we recover the stiffness matrix of local elasticity for materials with *λ*=*μ* by formally setting9
3.39
in equation (3.37). This demonstrates the convergence of the discretized non-local formulation towards the discretized local formulation.

### (c) Integral representation of the solution

Using the notation in appendix A*b*, we apply the Laplace transformation with respect to time *t* to equation (3.15) and find the transformed solution
3.40
3.41
Knowing the individual Laplace transforms
3.42
and
3.43
we can use the convolution theorem of Laplace transforms to obtain the solution in Fourier space
3.44
Finally, we use the convolution theorem of Fourier transforms to obtain the following integral representation of the solution of equation (3.5) in (** x**,

*t*) space 3.45 with Green’s tensor10 3.46 with and where the notation has been used. Substituting equations (3.45) and (3.46) into equation (3.5) confirms that the equation of motion is identically satisfied. Since

**(**

*g***,0)=**

*x***0**and , the initial conditions are satisfied as well.

## 4. Examples

In this section, we consider the following examples:

(i) transient dynamics of a bar (one dimension),

(ii) initial value problem for initial data with Fourier-series representation (three dimensions), and

(iii) static solution for a single point load (three dimensions).

### (a) Transient dynamics of a bar (one dimension)

In the following, we consider a bar at rest for *t*≤0, subject to a pair of self-equilibrated point loads for *t*>0
4.1
This problem has been the topic of previous research. In Silling (2003) the *static* analytical solution for less smooth micromodulus functions is discussed. In Bobaru *et al.* (2009), the static numerical solution using dynamic relaxation is discussed in the context of adaptive grid refinement. In this paper, we discuss the fully *dynamic* response. We begin by introducing the normalization
where the displacement field has been normalized with the static elongation of a bar of length 2*L* subjected to a static force *F* in local elasticity. Then the solution for both local and non-local formulation is given by equations (2.5) and (4.1)
4.2
From equation (4.2), we see that the displacements are antisymmetric, .

#### (i) Classical local solution

Local elasticity corresponds to a linear dispersion relation, . In this case we can evaluate equation (4.2) explicitly. The solution can be written as the superposition of the transient solution and the static solution :
4.3
with
4.4
and
4.5
where sgn(*x*)=*x*/|*x*| is the sign function. Because the initial conditions are not in equilibrium with the applied static load, the resulting displacement field is time dependent. However, if, for example, only the spatial interval is considered, then the displacements become independent of time for . In general, we have , while is the static solution for equations (2.1), (2.2), and (4.1). This is illustrated in figure 3.

#### (ii) Peridynamic non-local solution

Consider the following peridynamic material, characterized by the micro-modulus function:11
4.6
In the following, we normalize the *material* length scale associated with the non-local material model, the horizon *δ*, with the *geometrical* length scale of the bar, *L*: .

The static solution defined in equations (2.1), (2.2), and (4.1) can be found analytically using Fourier transformations, 4.7 4.8 and 4.9 We can solve the integrals required for the inverse transformation analytically when identifying the presence of the one-dimensional Dirac distribution in the solution 4.10 where is the integral sine function. Returning to the dynamic problem, we find 4.11 with 4.12 Note that unlike in the local theory the ‘transient’ response does not vanish for large . Instead, we have undamped oscillations with the normalized angular frequency and amplitude around the static solution . 4.13

#### (iii) Comparison of classical and peridynamic solution

As previously pointed out by Silling *et al*. (2003), unlike in local elasticity, the presence of a Dirac distribution in the external loads in peridynamics leads to a Dirac distribution in the displacement field under the point load.12 In this sense, the local formulation seems more smooth. However, for all other points, the opposite is true: while the strain field has a jump discontinuity in local elasticity (dashed line) it remains continuous in peridynamics, as illustrated in figure 4.

In the local formulation, a static point load eventually (i.e. for ) leads to a static deformation field (see equation (4.5)). This is not true in peridynamics where in the long time (LT) limit any point continues to oscillate around the static solution with frequency and amplitude as shown in equation (4.13). In the limit as , the frequency of these oscillations becomes infinite while the amplitude goes to zero: . At the same time, the static non-local solution converges to the static local solution (4.5), . In figure 5, the local displacement under the point load is plotted together with the static non-local solution, , the dynamic non-local solution as well as the amplitude of the predicted oscillations.

According to equations (2.4) and (4.1), the displacement under the point load can also be interpreted as the total energy of the system: *E*_{tot}(*t*)=(2*F*/*A*) *u*(*x*=*L*,*t*). In local elasticity, the total energy becomes constant for and the system is conservative. In peridynamics, the point under the point load never comes to rest so the point forces continue to change the total energy.

### (b) Initial value problem for initial data with Fourier-series representation (three dimensions)

As an example, consider the following initial displacement:13 4.14 Then the solution follows from equation (3.44): 4.15 and 4.16 While the spatial distribution is identical in local and non-local elasticity, the temporal dependence is generally different due to the nonlinear dispersion relation.

### (c) Static solution for a single point load (three dimensions)

The non-local static solution of equation (3.15)
4.17
subjected to a point load at the origin, ** b**(

**)=**

*x*

*P**δ*(

**), is given by 4.18 with 4.19 where**

*x**x*=||

**||. Substituting the component of the local acoustic tensor equations (3.17) and (3.18), we obtain the well-known result (see Eason**

*x**et al*. (1956)) 4.20 In the non-local case, the integrals given in equations (4.19) do not converge: the last term in the integrands is unbounded for large

*k*since, unlike in the local case, the acoustic tensor is constant in this limit: . The reason for this divergence is the presence of a Dirac distribution in the solution. We rewrite equations (4.19) explicitly revealing the presence of a Dirac distribution in Green’s function in the non-local formulation 4.21 and 4.22 As an example, we consider the exponential micromodulus function discussed in §3

*b*(v). The results of the numerical integration is presented in figure 6 for different degrees of non-locality

*δ*>0 , together with the local solution

*δ*=0 with

*λ*=

*μ*. In figure 6

*a*, the normalized displacement component in the direction of

**is shown while in figure 6**

*x**b*the normalized displacement component orthogonal to

**is shown.**

*x*This numerical study indicates that in the limit of vanishing horizon *δ*→0 the non-local solution converges towards the local solution almost everywhere. However, for any finite horizon *δ*>0, the displacements under the point load remain bounded14 when ignoring the Dirac distribution.,15 In contrast, the displacements in local elasticity are unbounded due to the presence of the 1/*x* singularity in the solution. This is perhaps not surprising as the motivation in some of the earlier work on weakly non-local methods was to remove the presence of the unphysical singularities in the stress field surrounding a crack tip.

## 5. Conclusions

We derive an integral representation for the solution of the linear elastic three-dimensional bond-based peridynamic equation of motion using Green’s functions. Then we apply this theoretical result to a series of static and dynamic examples in one and three dimensions and compare the classical, local solution to the peridynamic, non-local solution. We found that when subjecting a peridynamic material to a point load represented by a Dirac distribution, the response will also include a Dirac distribution, independent of the spatial dimension of the problem. In general, it seems that the peridynamic solution is as smooth as the external forces/the initial data, as already pointed out by Silling *et al*. (2003) for the one-dimensional case. Applying a point load to a classical local material does not lead to a Dirac response, which suggest that in some sense the local response is more smooth. However, in the neighbourhood of the point load, the non-local deformations are actually more smooth: in one dimension, the strain field across a point load has a jump discontinuity within the local formulation but remains continuous in peridynamics. In three dimensions, the local displacements themselves have a 1/*x* singularity when approaching the point load while the non-local displacements remain finite.

Other results obtained in this paper include the derivation of the non-local three-dimensional jump condition, the weak formulation of peridynamics as well as the associated FE discretization.

There are several interesting and challenging directions for future research in non-local elasticity.

(i) More complex/realistic external forces such as time-dependent point loads (either at a fixed spatial position or moving in space) or spatially distributed loading.

(ii) Non-local deformation and stress field surrounding a crack tip.

(iii) More complex non-local material models, e.g. for anisotropic carbon fibre-reinforced polymer composites.

(iv) Non-local boundary conditions and their relationship to the local boundary conditions (Dirichlet, Neumann, Robin).

Studying these problems would provide important theoretical insight and increase the fidelity of peridynamic simulations.

## Acknowledgements

We would like to thank Prof. Florin Bobaru for a fruitful discussion at the ASME conference in Boston, Dr Rich Lehoucq, Dipl. Math. Dipl. Ing. Juliane Dunkel for a careful and constructive critical reading, Dr Michael Bieterman for initial finanical support and Dr Greg Shubin for loaning us his copy of Aki & Richards (2002). Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

## Appendix A

#### (a) Fourier transforms

##### (i) Definition

In the Cartesian basis {*e*_{i}}, the Fourier transform of the vector-valued function of the three independent variables *x*_{j}=*e*_{j}·** x**,

*j*=1, 2, 3, is defined as A 1 and A 2 or in component form and

##### (ii) Fourier transforms of derivatives

The Fourier transform of the gradient of a tensor ** T** is given by
A 3
where ⊗∈{·, ,×} is the dot product, dyadic product or the cross product. As an example, the Fourier transform of the derivative of a scalar function

*f*(

*x*,

*y*,

*z*) with respect to

*x*is given by .

##### (iii) Convolution theorem

The convolution of the two functions ** f**(

**) and**

*x***(**

*g***) is defined as follows: A 4 where ⊗∈{·, ,×} is the dot product, dyadic product or the cross product between the two vector-valued functions**

*x***and**

*f***. Then the Fourier transform of**

*g***(**

*h***) is given by .**

*x*#### (b) Laplace transforms

##### (i) Definition

In Cartesian basis {*e*_{i}}, the Laplace transform of the vector-valued function ** f**(

*t*) is defined as A 5 and A 6

##### (ii) Laplace transforms of derivatives

A 7

##### (iii) Convolution theorem

The convolution of the two functions ** f**(

*t*) and

**(**

*g**t*) is defined as follows A 8 where ⊗∈{·, ,×} is the dot product, dyadic product or the cross product between the two vector-valued functions

**and**

*f***. Then the Laplace transform of**

*g***(**

*h***) is given by .**

*x*## Footnotes

↵1 Equation (2.1) has the same form for both local and non-local formulation. However, the operator is different and so is the solution. This will be indicated by the superscripts

^{L}(·),^{NL}(·) in the following.↵2 It is physically intuitive that as the distance between a pair of particles gets very large, the interaction between them becomes negligible. In what follows, we shall assume that this happens fast enough to ensure the convergence of the various infinite integrals encountered. Convergence is trivial when the material has a ‘horizon’, i.e. if there is no interaction between particles that are more than some finite distance apart.

↵3 This can be verified by taking the time derivative inside the integral and substituting the corresponding equation of motion.

↵4 Here

*ϕ*^{+}and*ϕ*^{−}are the values offrom positive and negative sides of of and [[*ϕ*]]:=*ϕ**ϕ*^{+}−*ϕ*^{−}.↵5 See Eringen (1975) for a derivation.

↵6 In Weckner & Abeyaratne (2005), analogous results were obtained for the one-dimensional case.

↵7 This restriction is only present in the bond-based formulation and no longer present in the so-called state-based peridynamic formulation (see Silling 2000; Silling

*et al.*2007).↵8 The local case can be obtained by using the formula for the Fourier transforms of derivatives given in appendix A

*a*(ii). In the non-local case, one first uses the convolution theorem given in appendix A*a*(iii) to obtain the acoustic tensor . Carrying out the required integration for the inverse Fourier transform in spherical coordinates and using the result that , where is the unit sphere we obtain (3.15).↵9 Using the operator identity , we could alternatively set

*Λ*(*ξ*)=(15*λ*/2*π*)(Δ(*ξ*)′′/2*ξ*^{4}) or*Λ*(*ξ*)=−(15*λ*/2*π*)(Δ(*ξ*)′/*ξ*^{5}).↵10 Applying an external force which is localized in both time and space,

(*b*,*x**t*)=*ρ*Δ(*t*)Δ()*x*, to a body initially at rest leads to the solution*n*(*u*,*x**t*)=·*n*(*g*,*x**t*). Whenis non-zero over a finite spatial domain of size*b**r*_{b}, the ratio of the material length scale and the loading length scale*δ*/*r*_{b}determines the degree of non-locality in the solution. Additionally, the degree of non-locality decreases with increasing distance from relative to*δ*.↵11 This specific one-dimensional micromodulus function is the analogue of the three-dimensional trigonometric micromodulus function discussed in §3

*b*(v): both have a characteristic wavelength*χ*_{c}=2*π**δ*above which waves do not disperse. The mathematical simplicity of the resulting dispersion relation allows for analytical solutions.↵12 The presence of Dirac distributions in the deformation field is not unphysical because any physical loading function is always applied over a finite spatial domain.

↵13 Adding additional initial velocity is straight forward. Also this methodology can be easily extended to initial data with countably many wavenumbers, i.e. with Fourier-series representation.

↵14 In Kunin (1983), the analogous result is obtained for the static Green’s function within the context of the so-called quasicontinuum model of a perfect lattice of identical particles.

↵15 If we had applied the external force field to a finite region in space, the Dirac distribution would not be present in the solution.

- Received May 1, 2009.
- Accepted July 17, 2009.

- © 2009 The Royal Society