On the Computational Derivation of Bond-Based Peridynamic Stress Tensor

The concept of ‘contact stress’, as introduced by Cauchy, is a special case of a nonlocal stress tensor. In this work, the nonlocal stress tensor is derived through implementation of the bond-based formulation of peridynamics that uses an idealised model of interaction between points as bonds. The method is sufficiently general and can be implemented to study stress states in problems containing stress concentration, singularity, or discontinuities. Two case studies are presented, to study stress concentration around a circular hole in a square plate and conventionally singular stress fields in the vicinity of a sharp crack tip. The peridynamic stress tensor is compared with finite element approximations and available analytical solutions. It is shown that peridynamics is capable of capturing both shear and direct stresses and the results obtained correlate well with those obtained using analytical solutions and finite element approximations. A built-in MATLAB code is developed and used to construct a 2D peridynamic grid and subsequently approximate the solution of the peridynamic equation of motion. The stress tensor is then obtained using the tensorial product of bond force projections for bonds that geometrically pass through the point. To evaluate the accuracy of the predicted stresses near a crack tip, the J-integral value is computed using both a direct contour approximation and the equivalent domain integral method. In the formulation of the contour approximation, bond forces are used directly while the proposed peridynamic stress tensor is used for the domain method. The J-integral values computed are compared with those obtained by the commercial finite element package Abaqus 2018. The comparison provides an indication on the accurate prediction of the state of stress near the crack tip.


Introduction
Peridynamics is a rather recent reformulation of nonlocal continuum mechanics [1] which tries to unify, and incorporate within a single framework, the mathematical modelling of continuous media, discontinuities, cracks, and particle mechanics [2][3][4][5]. In this formulation of a solid mechanics problem, the spatial derivatives of related field quantities such as displacement, strain, or stress are not required. Thus, peridynamic theory replaces the governing partial differential equations of classical continuum mechanics by integro-differential equations, and in doing so removes singularities associated with classical continuum mechanics, e.g. at the tip of a sharp crack or right under a concentrated force. Peridynamics can be considered as the continuum version of molecular dynamics and the derivation of its associated parameters can be assimilated to that of molecular dynamics or atomistic simulations [6]. It has been shown in the literature that the peridynamic theory can be implemented in finite element as well as in molecular dynamics codes [4,7,8].
As a nonlocal continuum model, peridynamics has a different formulation and relevant quantities from classical continuum mechanics [9]. The link between peridynamics and classical continuum field theories of mechanics can, nonetheless; be established through the response variable fields, which provide, uniquely, the solution to the solid mechanics problem in classical theories viz. the stress, strain, and displacement fields. As the study of strain and displacement fields is essentially a geometric problem, the most important, yet non-trivial, link between the two models lies in the derivation of the stress tensor.
While the concept of a contact stress tensor was first introduced by Cauchy nearly two centuries ago, it can be viewed as a special case of the notion of stress in the more comprehensive paradigm of rational mechanics in which particles of matter can interact over a finite distance, a fact that was noted in the early twentieth century by Love [10] and elaborated later by researchers in the field [11][12][13]. It has been shown that all measures of stress, as are known in classical continuum mechanics, are interrelated and can be derived from a single set of quantities known as base forces [14]. As such, once the relation between a peridynamic stress tensor and one counterpart in classical continuum is established, the description of stress at a point within the nonlocal framework can be deemed as complete. The peridynamic stress tensor obtained as such must rely on one of the formulations of the nonlocal theory.
Peridynamics falls into two fundamental categories viz. bond-based and state-based. In bond-based peridynamics, the interaction between two particles (or material points [15]) is fully defined based on their relative position in the undeformed (base or reference Ω 0 ) and deformed (current Ω) configurations. The equation of motion at a point of a rate-independent, time-independent, and history-independent medium can in this case be re-written as follows: ρ x ð Þ :: Where ρ(x) is the mass density in the reference configuration, u(x, t) is the displacement vector field, f is a pairwise force function, which has dimensions of force per unit volume squared, that the particle at position x exerts on the particle at position x ′ , and b is the body force vector field. While it is intuitively expected that the integration should be over the entire volume of the medium, it is neither economic nor essential for this to be the case. As such, the concept of a horizon, i.e. H x = {x ′ ∈ Ω 0 , |x ′ − x| < δ} is introduced which signifies, as its radius, the finite distance over which two particles can interact and the long-range forces are non-negligible. The integration in equation above is thus performed over the finite volume of H x surrounding x, referred to as the neighbourhood or horizon of x. As the size of the horizon shrinks to zero nonlocal and classical formulations coalesce [13]. Bond-based peridynamics essentially suggest a Cauchy crystal with a reduced number of independent material constants (e.g. for a triclinic body 15 rather than 21). In state-based peridynamics, the interaction between two points is dependent upon the way the two interact with other points adjacent to them. As such, the limitations imposed on the cardinality of the set of elastic constants by the Cauchy crystal formulation is lifted, though the formulation renders calculations computationally more demanding. Drawing parallels with molecular dynamics, bond-based peridynamics is analogous to two-body interactions and state-based peridynamics to many-body interactions. As such, micro-potentials can be defined in peridynamics similar to two-body and many-body potentials of molecular dynamics [7,8].
Bond-based peridynamics is used as an analytical tool but also has been implemented in inhouse meshless and conventional finite element packages and is used to solve a wide range of problems. To mention but a few, Silling et al. [16] studied the deformation of a bar using peridynamics formulated using complex analytic functions. Bobaru et al. [17,18] studied adaptive mesh convergence in numerical simulation of 1D and 2D problems. Askari et al. [19] reported on multiscale simulation of materials using peridynamics. Mechanics of heterogeneous media, as such, can be formulated within the framework of this nonlocal method as denoted by Alali and Lipton [20]. Other researchers studied aspects of composite mechanics for heterogeneous or anisotropic media [21][22][23][24][25][26]. Macek and Silling [4] provided examples on damage due to high-velocity impact in plates and membranes. They looked into perforation and penetration mechanisms and concluded that the finite element package Abaqus was much faster in solving these problems compared to the in-house developed meshless code EMU. Cracks in 2D anisotropic media as well as trans-granular and inter-granular cracks in polycrystalline structures were also studied by researchers and good correlation with experiment was achieved [25]. Finally, ceramic-metal bi-materials with cracks at the interface as well as in the brittle ceramic part were studied by bond-based peridynamics implemented using truss elements in Abaqus by Beckmann et al. [27].
As stated before, constitutive modelling in bond-based peridynamics has limitations. The Poisson's ratio is not independent and cannot be assigned any permissible value. For 3D and 2D plane strain problems, the Poisson's ratio is ν = 1/4 while for plane stress problems ν = 1/3. Statebased peridynamics avoids this limitation by allowing the pairwise force function f for a single interaction at a point to be a function of relative positions and displacements of all points adjacent to that point. Constitutive modelling in state-based has been discussed in detail by Silling et al. [28].
Bonds have been implemented through truss elements in commercial finite element software [4,27] as well as through beam elements [29]. Every point in such a model is connected to all points within its horizon with a single element of the relevant type. Once the elastic stiffness and fracture-related parameters between the continuum model and the truss mesh are established, the truss mesh can be used as the nonlocal peridynamic mesh to study the problem. The quantisation of truss element properties is not a unique process and there are several rational ways to achieve this.
As for response quantities, stress, strain, and displacement fields draw a parallel between the peridynamic and continuum models. While being considered as essential response parameters in solid mechanics, stress and strain are not intrinsic to peridynamics. In a bond-based peridynamic pseudo-mesh built through embedded trusses, the displacement fields of the two models can be related directly, i.e. the values of displacement field variables at discrete points of the peridynamic pseudolattice can be compared to their counterpart in the continuum model or finite element mesh. Strain at a point is attributed to the symmetric part of the displacement gradient field and can either be obtained through numerical differentiation of the displacement field or through solving for strain components the stretches (direct strains in inclined fibres) in the truss elements passing through the point. As for the stress tensor, which relates peridynamic bond forces to surface tractions in the continuum model, direct force projection and force averaging through strain energy-internal energy equivalence are the two methods that render it obtainable.
Fracture and damage have also been studied in peridynamics. Indeed the strength of the peridynamic method lies in the accurate prediction of fracture [4]. Based on Griffith's fracture criterion [30], as a necessary condition for crack growth, Aidun and Silling [31] studied dynamic fracture in quasi-brittle materials using peridynamics. Agwai et al. in a recent study compared existing crack propagation models with peridynamics [32]. There are several issues when it comes to the prediction of crack growth, crack velocity, and path based on classical continuum theory. It is wellestablished that due to the singularity of stress and strain fields at the crack tip in the classical theory, one must resort to energy-based propagation criteria. J-integral [33] could thus be used to predict the propagation of fracture in brittle materials. It has, in fact, been shown that J-integral can also be obtained based on peridynamics [34]. As for the crack velocity and propagation direction, it has been shown, time and again, that methods such as the Extended Finite Element Method (XFEM) and Continuum Damage Mechanics (CDM) can lead to erroneous crack velocities and mesh-sensitive crack paths [27] to some extent. There is another issue with crack paths in classical continuum theory which should be addressed. The direction of a crack in a 2D problem depends solely on the ratio of modal stress intensity factors K II /K I . Nonetheless, the stability of the predicted path is contingent upon the degree of stress biaxiality. This requires introduction of an extraneous parameter, T-stress [35], which, if taken into account, can determine the stability or instability of the predicted path but does not provide any further information as to what the actual crack path will be. Branching and coalescing cracks can also pose issues when dealing with curved or kinked cracks. None of these phenomena poses an issue in peridynamics and no special consideration is needed to deal with them. Based on fracture mechanics, a damage parameter can be introduced consistent with its continuum counterpart [4].
The present study deals with the numerical derivation of the peridynamic stress tensor. A built-in MATLAB code is developed and used to construct a 2D peridynamic grid and subsequently approximate the solution of the peridynamic equation of motion. The stress tensor is then obtained using the tensorial product of bond force projections for bonds that geometrically pass through the point. The main contributions of the present work are: A straightforward computational approach to extract stress fields from a peridynamic code is presented.
The stress field approximated is compared with analytical and finite element solutions for two problems: a plate with a hole and a plate with a pre-existing crack. The results are in agreement with the expected values; however, the stresses suffer from the skin effect of peridynamics near the boundaries.
An approximation of the J-integral using the equivalent domain method is also presented for the estimation of the stress state near the crack tip. Compared to employing the direct integration method presented in [34], better convergence to the finite element J-integral value is obtained.
The proposed method can easily be extended to 3D problems and combined with the inherent ability of peridynamics to accommodate discontinuous displacement fields, the stress state near the crack tip during propagation can be monitored.
The present work is structured as: a brief overview of bond-based peridynamics and its formulation is given in Section 2. In Section 3, the theory of peridynamic stress tensor is covered and it is shown how the different components of stress may be obtained using an integration scheme across a surface. In Section 4 two examples are used to evaluate the performance of the proposed expression, a plate with a hole and a plate with a pre-existing edge crack. The examples are used to compare the stresses predicted using the peridynamic model with those computed from analytical expressions and finite element approximations, either directly and indirectly. Both problems are solved assuming linear elastic behaviour under plane stress conditions. Finally, concluding remarks are included in Section 5. The present work does not deal with a detailed study of displacement and strain fields; neither does it include geometric, material, or boundary nonlinearities. As such, plasticity, visco-elasticity, and contact force-induced stresses are not considered here.

Bond-Based Peridynamics
As mentioned above, the peridynamic model is a framework for nonlocal continuum mechanics based on the idea that pairs of particles exert forces on each other across a finite distance. Let us suppose two points within the horizon of each other are placed on two sides of an imaginary plane in a stress-free and undeformed medium. As the body undergoes deformation, the bond between the two points stretches (elongates or contracts) and a force in the bond is generated. When the bond is cut through by the plane, the force has several projections at the point on the plane through which it passes (where it intersects the plane). This simple fact is the basis of the bond-based peridynamic stress tensor. A peridynamic stress tensor was first derived in [12], despite the fact that the basis for such derivation existed for some time [10]. From the continuum mechanics perspective, the peridynamic stress tensor is analogous to the first Piola-Kirchhoff stress tensor of the classical theory. This fact is, nonetheless, of little significance when small deflections are concerned as all measure of stress in that case coalesce. The peridynamic stress tensor provides the force per unit area across any imaginary internal surface through which the bonds pass. Conceptually, the nonlocality of the peridynamic stress tensor is due to nonlocal forces in the bonds that cross from one side of the surface to the other. Since the peridynamic operator for the internal force density is expressed exactly as the divergence of the peridynamic stress tensor, the peridynamic equation of motion turns formally into the classical Eq. (12).
Let us now proceed to derive the bond-based stress tensor. We shall denote, as is customary, by ξ = x ′ − x the relative position of two points in the undeformed configuration, and by η = u(x ′ , t) − u(x, t) = u ′ − u the relative displacement of the two points at time t ≥ 0.
It is assumed that a scalar potential field function W(η, ξ) exists which contains all necessary information for the derivation. The bond force can then be obtained as follows, where w(η, ξ) in Eq. (1) is the micro-potential to which some works of literature have alluded (see e.g. [4,24]).
This model may be implemented using simple truss elements; however, in the implementation of this model using linear truss elements there is a restriction on the form W(η, ξ) may assume as long as its dependence on ξ is concerned. The energy stored in a single bond represented as a linear truss element is obtained in the sequel. Let the magnitude of the bond force be given by Eq. (2).
where E and A are the Young's modulus and cross-sectional area of the truss element, respectively, and n is the unit vector parallel to the deformed bond, i.e. n¼ ξþη j ξþη j jj . This implies the micro-potential w(η, ξ) is of the quadratic form with respect to η.n which renders W(η, ξ) to be of the form presented by Eq. (3).
Other features of f(η, ξ) would be direct consequences of the definition, the specific type of formulation or master balance laws, and are summarised as follows: Equation (4a) arises from the definition of the 'horizon'. Equation (4b) is the implication of Newton's third law and Eq. (4c) originates from the balance of angular momentum. The deformation of a single bond in the truss system is defined using the stretch parameter of Eq. (5) as follows: This definition yields the same value as that of the direct engineering strain in the bond and has a similar function. To formulate the problem, a constitutive relation is required. The magnitude of the pairwise force function f(η, ξ) is hence related to the stretch in the bond which for the linear elastic case is as follows: where c in Eq. (6) is the micro-modulus and is related to the elastic constants of the medium. The parameter c can be obtained by equating strain energies under isotropic expansion of the classical continuum model and that of the peridynamic model for a single horizon size. It has been shown that if c assumes a conical variation with respect to the relative position ξ of two particles, better convergence to classical elasticity is achieved [3]. Assuming plane stress conditions, the micro-modulus is given as: Ignoring the inertia effects, the peridynamic equation of motion can be approximated numerically using a midpoint collocation method, as described in [36]. Then, the integral in Eq. (1) can be written as a finite summation as: where N is the number of particles such that |x j − x i | < δ and V j is the material volume associated with particle j. Although the expression in Eq. (6) describes a linear relationship between bond stretches and bond forces, the final system of equations is nonlinear. Here, the Newton-Raphson method is implemented to solve the final system of Eq. [37].

Derivation of the Bond-Based Peridynamic Stress Tensor
Let us now move on to derive the peridynamic stress tensor based on the bond-based interactions between the pairs of adjacent points. As Macek and Silling [4] demonstrated, standard truss elements available in Abaqus can be used to represent peridynamic bonds. The stress tensor is essentially a dyad relating the bond forces and radius vectors of adjacent points for bonds passing through a point to the state of stress at that point. Figure 1 shows the schematic of a truss-like interaction in which every point is connected to its adjacent ones within its horizon using rod/spring/truss elements. A node incident to such an element is thus an end point of that element. This figure also highlights the topology of the system. In the simplest case where the peridynamic model is isotropic and homogeneous, the node spacing in the mesh is uniform and a lattice constant Δx can be used to represent the entire nodal set. The horizon is always an integral multiplier of the lattice constant. Two types of convergence are then conceivable: an m-convergence and a δ-convergence. Given δ = mΔx which relates the horizon size to the lattice constant as an integer multiplier of the former, one can conduct an m-convergence study in which the size of the horizon enlarges (or shrinks) while the lattice constant remains the same, or a δ-convergence study in which both the lattice constant and horizon size increase (or decrease) proportionally. In each case, the mechanical properties of the truss system must be related to the bulk matter through the equations presented in the literature (see e.g. [1,13]).
Taylor expansion of the pairwise force function about an arbitrary points leads to the following equation: where o( * ) denotes the lowercase Landau order symbol. Ignoring the higher order terms, the expansion about the point of zero deformation leads to the following equation: where ξ 0 signifies the initial bond vector. If the original configuration corresponds to the forcefree (hence stress-free) state, f(0, ξ 0 ) is identically zero. The constitutive micro-modulus tensor is then defined as C ¼ ∂ f ∂η 0; ξ 0 ð Þ. In a linearized theory, where the bonds are defined as entities of fixed length ab initio and then deformed (stretched or contracted), the last term in Eq. (10) disappears as it related the variations in pairwise force function to perturbations in initial geometry. Hence, the pairwise force function takes the following familiar form: It can be shown (see e.g. [4]) that C(ξ) is obtained as follows: where c is given by Eq. (7). The basis for numerical calculation of stress in this work is the bond force in the form presented by Eq. (11), however, the formulation is general. As it was shown by Lehoucq and Silling [12], the peridynamic stress tensor can be defined as follows: Here, S denotes a unit sphere and dΩ m denotes a differential solid angle on S in the direction of any unit vector m. They provided classical, mechanical, and variational interpretations of this quantity and showed its divergence-free nature (∇.v(x) = 0) which renders it analogous to the classical Cauchy stress tensor. They also showed that for a finite number of discontinuities, the definition of stress remains the same and established its uniqueness conditions. The details of the relevant theorems and their proofs can be found in [12]. We present here merely the most related interpretation, namely mechanical interpretation of the force flux. Referring to Fig. 2, let the definition of Eq. (13) hold. The associated traction vector corresponding to a plane with normal n is thus obtained as follows: Following [12] and as depicted in Fig. 2, let y = x + ym and z = x − zm. The differential area dA y is defined on a sphere centred at z and containing y that subtends a differential solid angle shown as dΩ. Thus: In a similar fashion, the analogous quantity on a sphere centred at y and containing z is indeed identical: Now if a plane with normal n cuts through the cylinder of cross-sectional area dA y = dA z at point x, the element of area is obtained as follows: Thus, the total force that element of volume dA y dy exerts on the element of volume dA z dz is as follows: Fig. 2 Interpretation of stress at the point x across a plane with unit normal n (adopted from [12]) Thus the differential force per unit area on the plane through x with unit normal n is: Which constitutes the integrand of Eq. (14). The factor ½ appears before this equation due to the fact that bond forces are counted for twice in their enumeration (handshaking lemma). This definition provides the mechanical interpretation of stress and is analogous to the one of elasticity as discussed by researchers [2,12]. With regard to computational calculation of stress, a discretisation scheme will render the system essentially a discrete one. When truss elements are used to implement the peridynamic mesh the end nodes of each element are the discrete particles and the elements represent the interaction between them. Just as in a discrete atomic system with a particular form of potential energy, as discussed before, the system comprises interacting nodes and interactions are of classical type.
Following Eq. (13), the components of the nonlocal stress tensor are as follows: Now we consider this equation applied to a discrete system. Consider a point centred at R and let C(R) be the Voronoi cell of R. We define the function Δ(r − R) as follows: It is obvious that where v R is the volume of the Voronoi cell situated at point R. For simplicity, it is assumed that the volume of each Voronoi cell is constant. Furthermore, the numerical examples presented in this study are discretized using a structured grid of particles of equal volume. For any piecewise continuous function g(r), the average is defined as follows: This average value can be assumed to be the value of g(r) associated with the Voronoi cell C(R). Re-writing Eq. (18) leads to the following definition for stress (designated by σ(R) = σ ij (R)e i ⨂ e j as is conventional) at the point R: To make the transition from a continuum to a discrete system, the pairwise force density function, derived as follows, is used: where F(R 1 , R 2 ) is the force that particle at R 1 exerts on the particle at R 2 . It is always possible to find this force. Analogous to what transpires in a molecular structure, once the system is in equilibrium, the force that atom at R 1 exerts on the atom at R 2 is equal in magnitude and opposite in direction to the force that the rest of the system exerts on the atom at R 2 . Following Eq. (22) and by replacing r with R + αr and r ′ by R − (1 − α)r one obtains: Substituting from Eq. (23) into Eq. (21) and after some algebraic manipulation, the stress at R is obtained as follows: Equation (24) always applies to a discrete atomic system as well as a system of interacting particles. The force between pairs of particles is the gradient of a higher order potential the existence of which is assumed and the form of which is contingent upon the type of interaction. This equation is the basis for the calculations of this work and is analogous to the BDT atomistic stress tensor defined by Basinski, Duesbery, and Taylor [6]. The BDT stress tensor is derived as follows: Which, for the static case, coalesces with the one derived in the present work. Here, similar to the 'horizon' in peridynamics, a 'cut-off radius' defines the extent of spatial influence of a point upon another. Thus, Ω α is a sphere of radius equal to this cut-off radius. Interacting force between two particles, as mentioned before, is obtained as follows: The scalar field function ψ represnts the pairwise potential function. As the size of Ω α shrinks to zero fluctuations in stress emerge. This is a natural result of the definition of stress as an averaged quantity.
Dealing with the problem of fracture, the critical fracture energy (G c ) must be compared to the path independent J-integral defined as follows: where W is the elastic energy density, K the kinetic energy density, and T i denotes the traction force. For a sharp crack parallel to x 1 , J = G c is the necessary condition for the crack to propagate. For a linear elastic solid under static loading, J can be re-defined as follows: Comparison of J E and G c provides a robust and universal method, within the framework of linear elastic fracture mechanics (LEFM), for determining whether a crack will be propagating or not. Cracks can move or arrest during a loading process based on this comparison. Another important and relevant concept is that of T-stress which determines the stability of a crack path. It arises from the asymptotic expansion of the 2D stress field around the crack tip in a linear elastic body as follows: where r and θ designates the distance from the crack tip and the angle at which the point is situated in the polar coordinate system. The relation as presented here is valid when crack faces are stress-free. As Cotterell and Rice [35] discussed in detail for small amounts of Mode I crack growth, a straight crack path is stable when T-stress is negative (T < 0), whereas the path will be unstable and, hence, could deviate from being straight when T-stress assumes positive values (T > 0). As such, the point at which the crack length renders a sign shift in this quantity from negative to positive can be reckoned of as the point of possible initiation of kinking.

Case Studies
In this section, the peridynamic stress tensor derived earlier is used in two case studies: (i) a plate containing a hole and (ii) a plate containing a notch. The first example is used to illustrate the convergence of the stresses computed using the peridynamic model to the stresses computed from an analytical solution using Airy stress functions and those predicted using the commercial finite element software Abaqus. In the second example, the peridynamic stresses are used to compute the path-independent contour J-integral using the Equivalent Domain method [38].

Stress Concentration Around a Hole
The familiar problem of a plate containing a hole is selected as it exhibits both normal and shear stresses, with an analytical solution for an infinite plate being available that renders asymptotic comparison possible. Furthermore, the accuracy of peridynamic stresses near geometrical boundaries that cause stress concentration is evaluated through comparison with the finite element and analytical solutions. The plate is assumed to behave elastically under plane stress conditions. The problem geometry is illustrated in Fig. 3 while the problem parameters are summarized in Table 1. Use of bond-based peridynamics restricts the value of the Poisson ratio to v = 1/3. Using Airy stress functions, the analytical relations for stress components around a hole in an elastic infinite medium are given as: where r and θ are the polar coordinates measured from the center of the circle. Transformation to a Cartesian system of coordinates is carried out with the familiar tensor transformation laws as defined in Eq. (31), using both indicial and matrix notation: where A is the affine transformation matrix, which for a 2D case is given as: The stresses computed though the analytical solution of Eq. (30a-c) is compared with the stresses computed using the finite element software Abaqus and the stresses predicted using an in-house peridynamic code written in MATLAB. The stresses are computed as a postprocessing step through Eq. (24). For the finite element, solution 8-node biquadratic plane stress quadrilateral elements were used for the discretization of the computational domain. The mesh was selected after convergence tests to ensure accuracy. To reduce the computational cost, the mesh was refined near the hole. In total, 58,239 nodes where defined leading to a system with 116,478 degrees of freedom.
The peridynamic model is set-up by assuming a uniform discretization of the problem domain. The computational efficiency could be improved by locally refining the grid. For  simplicity, however, uniform grids are used throughout this study. The interested reader can find further information on the refinement of peridynamic grids in other works of literature (e.g. [17,18]). The boundary conditions on the left edge of the plate were imposed through the addition of a fictitious material layer of thickness equal to the peridynamic horizon δ, as suggested in [39]. The value of the peridynamic horizon was set equal to δ = 3Δx PD i.e. m = 3, which is common for macroscale problems [40]. The formulation of peridynamics assumes particles completely embedded within the material. This assumption is violated near the geometric boundaries where the so-called skin effect of peridynamics appears. The skin effect is reduced by employing the volume correction method defined in [41]. The discretization Fig. 4 Comparison of the normal and shear stress fields obtained by the plane stress finite element and peridynamic analyses of plate with a hole length used for the peridynamic model is Δx PD = 0.0833 mm leading to approximately 36,0000 particles and 72,0000 degrees of freedom. The normal and shear stress fields from the finite element and peridynamic analyses are compared in Fig. 4. The peridynamic model can capture the stress concentration near the hole and a general agreement between the stress fields can be observed in all directions. During the derivation of the peridynamic stress tensor, however, the assumption that the particle is surrounded by other particles has been made again. This assumption is violated not only for the particles near the edges of the plate but also for the particles near the hole.
To better illustrate the accuracy of the peridynamic stress tensor, the normal stresses are plotted along a horizontal and a vertical section (Fig. 5 and Fig. 6 respectively) while the shear stresses along a circular path with radius r 1 = 1.4r are shown in Fig. 7. Small deviations from the analytical solution are anticipated due to the finite dimensions of the problem considered. As expected, a stress concentration factor of approximately 3 is predicted near the hole from the finite element and the analytical solution. The stresses predicted using the peridynamic model, however, are affected in the vicinity of the boundaries and the hole. The peridynamic solution predicts a higher stress concentration factor of approximately 3.4. Moving towards the interior of the plate, the stress values predicted converge to those from the other solutions. Indicative of this are the shear stresses plotted away from the hole.
To better illustrate the accuracy of the proposed expression (Eq. 24), the relative error of the stress component σ xx between the peridynamic and the finite element simulation, defined as: is plotted in Fig. 8. Near the hole, where the stress concentration takes place, the error between the two solutions is significant (approximately 20%). This plot, however; exemplifies how moving away from the geometrical boundary defined by the hole reduces the erroneous phenomenon. Particularly, by moving a distance equal to the peridynamic horizon, where particles are completely surrounded by other particles, the error between the solutions reduces to 1.19%. It is noted that the comparison takes place on interpolated values as positions of the peridynamic particles do not coincide exactly with the positions of the finite element nodes.

Computation of the Nonlocal J-Integral Using the Peridynamic Stress Tensor
It is a well-established fact of classical continuum mechanics that stress components at the tip of sharp cracks in elastic quasi-brittle materials are singular. In linear elastic materials, the order of singularity is of square root, i.e. σ ij r; θ ð Þ ¼ r − 1 2 F ij θ ð Þ, where r is the radial distance from the crack tip. As such, the value of a stress component cannot determine whether or not a crack is propagating. Extraneous parameters such as stress intensity factors (K-values) or Jintegral are required to study the propagation of cracks. Not only can the J-integral or stress intensity factors be used as flag parameters to dictate the status of a crack with regard to arrest or propagation but also the modal stress intensity factors are used to determine the crack path. For instance, in the case of maximum tangential stress criterion for crack propagation, the angle of crack rotation respective to initial crack plane is given as follows: Fig. 6 Comparison of the normal stresses along the horizontal section indicated with an orange line Fig. 7 Comparison of the shear stresses along a circular path with radius r 1 = 1.4r, indicated with an orange circle where χ ¼ K I = K II is the ratio of stress intensity factors. As it is clear from Eq. (34), the path of a crack depends merely upon this ratio and is independent of other geometric ratios of the problem. This is only correct when a notch size is small compared to other meaningful geometric dimensions involved. As crack grows larger, the non-singular term in the expansion of stress is required to determine the stability of the predicted path.
Peridynamics on the other hand simulates crack propagation and initiation naturally without the need of any other external criteria [42]. Furthermore, the nonlocal formulation of peridynamics avoids spatial derivatives and can thus accommodate naturally discontinuous displacement fields. The exact position of the crack tip does not need to be tracked explicitly at each iteration nor is a special treatment required for crack branching or coalescence, avoiding the need for complex algorithms to define the crack path. These attributes make peridynamics very suitable for fracture problems.  Consider the example of a plate with a pre-existing crack at one of its edges under plane stress conditions (Fig. 9). The plate is fixed on the other edge while a uniformly distributed load is applied on the top and bottom faces. The dimensions of the plate are H = 50 mm and L = 100 mm and the length of the crack is a = 10 mm. The material properties are the same as the example in the previous case study for the round hole in a plate.
The plate with a pre-existing crack again solved with both peridynamics and finite elements to compare the results. For the finite element mesh, special crack tip (collapsed) elements are used. Since elastic material behaviour is assumed, the singular fields are expected to exhibit 1= ffiffi r p singularity. This type of singularity can be approximated accurately by placing the middle node at a distance 1/3 from the tip position [43]. The domain was discretized by 14,975 biquadratic plane stress quadrilateral elements leading to a total 45,450 nodes. The mesh refinement near the crack tip region is used to improve computational efficiency.
A uniform discretization was implemented for the peridynamic model. Additionally, a fictitious material layer equal to the peridynamic horizon, δ, is used to apply the boundary conditions. In total, 721,800 peridynamic particles are used for the discretization leading to a total of 1,443,600 degrees of freedom. The discretization length is Δx = 0.0833 mm and the peridynamic horizon was set equal to δ = 3Δx. The pre-existing crack was introduced by breaking the bonds that cross the crack position prior to the simulation. An illustration of the grid employed is presented in Fig. 10.
The simulated stress fields from the finite element and the peridynamic models are compared in Fig. 11. To enhance clarity, the stress distribution near the crack tip is presented. There is a very close agreement in all stress components between the two solutions. Exactly at the crack tip, the stresses obtained by the two solutions deviate. Two factors contribute to this effect: the first has to do with the influence of removing bonds from particles that was discussed previously and the second has to do with the nonlocal nature of peridynamics. The internal length scale, introduced through the peridynamic horizon, leads to the stresses at the crack tip to remain finite due to the nonlocal particle interactions across a finite distance. Analogous to the example of the plate with the hole, stresses converge to the same solution away from the tip. Referring to the coordinate system indicated in Fig. 11, we define two sections at the tip of the crack, one horizontal with x = x(x, H/2), 0 ≤ x ≤ a and one vertical with x = x(0, y), 0 ≤ y ≤ a. The normal stresses σ xx and σ yy are plotted in Fig. 12, along the horizontal and the vertical section, respectively. This comparison is purely qualitative near the crack tip Since a direct comparison between the stresses at the crack tip is not possible, the parameter J-integral was used as it can uniquely describe the stress state in problems where linear elasticity is employed. The concept of the J-contour integral was introduced by Rice [33] and has met wide applicability in linear and nonlinear elastic fracture mechanics. Assuming that the faces of a crack are stress free and static conditions apply, the J-integral on an arbitrary path that contains the crack tip is given as [38]: where W is the strain energy density, T i are components of the traction vector and Γ is the integration contour. Hu et al. proposed in [34] a numerical procedure for the computation of J-integral using the peridynamic theory. The integration is performed on a closed contour ∂R that contains the crack tip as illustrated in Fig. 13. Two layers of thickness equal to the peridynamic horizon δ, called R 1 and R 2 , respectively, are defined on each side of the contour. Then the contour integral is computed as: where W i is the strain energy density of particle i, n c , n 1 , and n 2 are the number of particles on ∂R and in R 1 and R 2 , respectively, A k is the area associated with particle k, t is the material thickness, and n i, x is the horizontal component of the outward unit vector. The spatial derivatives in Eq. (36) can be approximated using a central difference scheme as suggested in [34] as: where Δx is the peridynamic discretization length. The numerical efficiency of the J-integral computation can be improved by converting the contour integral into a domain integral. Shih et al. [44] introduced the equivalent domain integral (EDI) method. This method is robust and applicable to problems with elastic, plastic, or viscoelastic material properties, with thermal loading, under static, or dynamic conditions [38]. Furthermore, numerical integration of the area or volume integral in 2D and 3D cases, respectively, is easier than contour or surface integration. As such, many commercial finite element packages (e.g. Abaqus) implement this method.
Using EDI, the contour integral of Eq. (36) is converted to an area integral using the divergence theorem by introducing an arbitrary function q(x). For a problem under static conditions, elastic material properties and stress-free crack face, the J-integral is computed as [38]: Fig. 12 Comparison of the normal stresses σ xx and σ yy along a horizontal and a vertical path from the crack tip respectively where A is the integration area and δ 1, j is the Kronecker delta. Typically, a pyramid or a trapezoidal function is defined for q. The function is constructed so its value is 0 at the edge of the domain and 1 at the crack tip. The domain integral in Eq. (38) requires the computation of the stresses using Eq. (24). A pre-existing crack in a peridynamic grid, however, is defined by removing the bonds it crosses. As such, near the crack faces inaccuracies in the computation of stresses will appear due to the same skin effect that was discussed previously. These inaccuracies will affect the value computed using Eq. (38). In their work, Hu et al. [34] also point out that this skin effect phenomenon near the crack faces is one of the reasons why the J-integral value computed using Eq. (36) deviates from the value obtained using classical elasticity. For this reason, the function q is assumed to be of trapezoidal shape, as illustrated in Fig. 14. The derivatives ∂q ∂x i vanish at the plateau and as such the J-integral from Eq. (38) will be insensitive to any inaccuracies near the tip of the crack.
The integral described in Eq. (38) can be approximated numerically as a finite summation using the particle positions as midpoint collocation points. It is then evaluated as: 2 Fig. 13 Definition of the integration contour ∂R and areas R 1 and R 2 near the crack tip for the computation of the nonlocal J-integral according to Hu et al. [34] where k tot is the total number of particles in the integration domain, σ ij, k is the stress component σ ij at the k th particle, and A k is the area associated with particle k. The spatial derivatives appearing in Eq. (39) can be approximated using the central difference scheme from Eq. (37). The J-integral value computed from Abaqus was J = 6.794 MPa • mm . This value was used to compare the accuracy of the value obtained by the peridynamic model. For comparison, both the J-integral approximations from Eqs. (36) and (39) will be used in the subsequent analyses, denoted as J cont and J dom , respectively. The motivation behind this comparison is two-fold. The first lies in the different formulation of the numerical approximations. The contour integral in Eq. (36) uses bond forces directly while for the one in the domain integral approximation (Eq. 39), the bond forces have been converted to stresses. Thus, it can provide an indication on the accurate prediction of the stresses. The second motivation is the easy extendibility of domain formulation to 3D cases. It is noted though that in neither formulation the bond forces or stresses near the crack tip are considered. In Fig. 17 in the Appendix, a flowchart is included that describes the computational procedure for calculation of the J-integral using either the contour or the domain integration formulation. It can be seen that the steps required for the evaluation of each expression are comparable.
The convergence of the two approximations is first evaluated by computing the relative error between the J-integral value computed using finite elements and the one from the peridynamic model as: where in J PD the values J cont and J dom are used. In Fig. 15, the relative error is plotted for different values of peridynamic discretization. The ratio of the peridynamic horizon to the discretization length is kept constant and equal to δ/Δx = 3. This type of analysis is termed δ-convergence as the peridynamic horizon becomes smaller as the grid is refined [3,36]. The path used for the integration was defined as a square centered at the crack tip with edge 2l c = 10, as illustrated in Fig. 13 and in Fig. 14 for the two numerical formulations. Both solutions approach the finite element solution with similar rates as the peridynamic discretization is refined. Additionally, the difference between the two expressions is less than 1%, with J dom exhibiting slightly better results.
In their work, Hu et al. [34] establish the path independence of J cont , even for contours very close to the crack tip, as long as the contours are feasible, i.e. it is possible to define the areas R 1 and R 2 . In practice, use of contours very close to the crack tip is avoided for the computation of the J-integral as they can lead to numerical inaccuracies. Discretizing the grid with Δx = 0.0833 mm, the analysis is repeated for different values of l c /a ∈ [0.05, 0.95]. The interval is selected so that the feasibility of the contour is satisfied both near the crack tip and the boundary of the plate.
The results for the different contour paths are plotted in Fig. 16. The maximum variation in the J-integral for the different paths are 0.7% and 0.3% for J cont and J dom , respectively. Although the improvement is small, J dom leads to approximations closer to the finite element solution compared to J cont . Additionally, assuming that function Fig. 15 Convergence of J cont and J dom to the J-integral value computed using finite elements, for different discretization values of the peridynamic domain and l c = 15. Logarithmic scale has been used for both the horizontal and vertical axes Fig. 16 Absolute error between the J-integral value computed using finite elements and J cont and J dom for different paths q(x) has a trapezoidal shape, better accuracy is achieved as the contour approaches the tip of the crack.

Conclusions
Using the work of Lehoucq and Silling [12], a numerical procedure for the calculation of stresses within the framework of bond-based peridynamics was proposed. The stress tensor is obtained using the tensorial product of bond force projections for bonds that geometrically pass through the point of interest. The peridynamic stress tensor was shown to be able to replicate situations with non-zero direct and shear stresses.
Two familiar yet important examples of plate with a hole and a cracked plate were presented to evaluate the accuracy of the stresses predicted by the peridynamic model. The results were compared with those obtained using finite element analysis and available analytical solutions. It was shown that the proposed numerical procedure, described in Eq. (24), is prone to inaccuracies near the geometrical boundaries of the domain, including the tip and faces of a crack. This effect is similar to the so-called peridynamic 'skin effect'. Nevertheless, it was shown that when the distance to the boundary is greater than the peridynamic horizon δ, stresses are captured accurately. In [45], Madenci et al. present a least square minimization methodology for the approximation of the temporal and special derivatives of a field variable. This approach is approach is able to solve ordinary and partial differential equations without being prone to errors near the boundaries. Employing this formulation could lead to a more accurate approximation of the stress field. In the present work, however, a similar 'approximate stress tensor' to that obtained by Silling et al. [46] is derived rather than the 'full' peridynamic stress tensor.
In the case of the cracked plate problem, the stress state predicted by the finite element method and the peridynamic solutions were compared using the J-integral. Using the peridynamic model, the J-integral value was approximated using both the direct contour integration expression suggested in [34] and the equivalent domain integral method. In both cases, the results were very closely correlated, with similar convergence rates. This indicates that the peridynamic and the finite element solutions converge to the same stress state around the crack tip. The domain integration method exhibits a small improvement in terms of convergence compared to the contour integration. The domain integration method is generally considered more robust and easier to extend to 3D problems that can be a valuable tool for future studies. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.