The J-area integral applied in peridynamics

The J -integral is in its original formulation expressed as a contour integral. The contour formulationwas, however, found cumbersome early on to apply in the finite element analysis, for which method the more directly applicable J -area integral formulation was later developed. In a previous study, we expressed the J -contour integral as a function of displacements only, to make the integral directly applicable in peridynamics (Stenström and Eriksson in Int J Fract 216:173– 183, 2019). In this article we extend thework to include the J -area integral by deriving it as a function of displacements only, to obtain the alternative method of calculating the J -integral in peridynamics as well. The properties of the area formulation are then compared with those of the contour formulation, using an exact analytical solution for an infinite plate with a central crack in Mode I loading. The results show that the J -area integral is less sensitive to local disturbances compared to the contour counterpart. However, peridynamic implementation is straightforward and of similar The original online version of this article was revised because of errors made in the production process regarding the layout. C. Stenström (B) · K. Eriksson Mechanics of Solid Materials, Luleå University of Technology, Luleå, Sweden e-mail: christer.stenstrom@ltu.se K. Eriksson e-mail: kjell.eriksson@ltu.se C. Stenström Luossavaara-Kiirunavaara Aktiebolag (LKAB), Luleå, Sweden scope for both formulations. In addition, discretization, effects of boundaries, both crack surfaces and other boundaries, and integration contour corners in peridynamics are considered.


Introduction
The J -integral is an expression for calculating the strain energy release rate in a cracked body, or the energy available at the tip of a crack to form new crack surfaces as the crack extends (Rice 1968). It is a firmly established parameter that can advantageously be applied both to linear and nonlinear elastic materials. In its original formulation, the J -integral is expressed as a contour integral. However, in finite element (FE) analysis, the J -integral may vary depending on the fracture, mesh design and integration method. One point of concern is that coordinates and displacements are calculated at element nodes, while stresses and strains are obtained at Gaussian integration points. Other points of concern are crack tip element specification, domain of integration and numerical errors. Therefore, the more objective J -area integral is commonly used in FE software, as its integration area dependence is less than the path dependence of a corresponding J -contour integral in integration regions close to a crack tip (Li et al. 1985;Banks-Sills and Sherman 1992;Kuna 2013). As an example of a FE implementation, using quadrilateral elements, the contour integral may be evaluated along a line of three Gaussian points in each element along the integration path, while the area integral includes the area of an element and all of its nine Gaussian points (Banks-Sills andSherman 1986, 1992;Khoei 2015). The J -area integral is generally also simpler to implement in FE codes (Anderson 2005). Nevertheless, the first ring-shaped areas close to and around the crack tip of the J -area integral usually still display path dependency (Khoei 2015, Fig. 7.23;Simulia 2014;Ansys 2019).
The peridynamic theory is a nonlocal formulation of solid mechanics, introduced for handling crack initiation, extension and final failure of a body, without the need of supplementary methods (Silling 2000;Silling and Askari 2005). Peridynamics is based upon integral equations, thereby avoiding spatial derivatives, which are not defined at discontinuities, such as crack surfaces.
A peridynamic nonlocal J -integral has been derived by Silling and Lehoucq (2010) for state-based peridynamics, based on an energy balance approach. Later, Hu et al. (2012) presented a bond-based peridynamic J -integral, using an infinitesimal virtual crack extension method. These nonlocal J -integral formulations include displacement derivatives and force interactions over inner and outer regions associated with the contour integral. The width of these regions depends on the degree of nonlocality of the peridynamic model.
The displacements are the principal unknowns in peridynamics, from which other quantities subsequently can be obtained. A general expression of the J -contour integral as a function of displacement derivatives was derived by Bruck (1989). The expression can also be found in Sutton et al. (1992), Hamam et al. (2007) and Hart and Bruck (2021). Since the J -area integral is common in classical solid mechanics, we implement in this paper the J -area integral in peridynamics as an alternative to the J -contour integral. Our treatment of J is simpler than Bruck's, due to rectangular integration paths.
For assessing the accuracy of the proposed J -area integral, we apply it to the exact analytical stress-straindisplacement solution of an infinite plate with a central crack loaded in mode I and compare the result to the corresponding analytical solution J 0 = K 2 I /E. The stresses or displacements of the exact solution on a contour away from the crack tip can be used as input boundary conditions in a peridynamic model of the same crack problem. To allow comparison with previous work (Hu et al. 2012;Stenström and Eriksson 2019), stresses are used as boundary conditions in the peridynamic model in this work. The J -area formulation will then be evaluated for the peridynamic model, to test for discrepancy between the analytical J -area integral and the peridynamic solution.
Our approach can be seen as a two-step process. In the first step, the displacement of a body may be obtained with any suitable method. In the second step, the J -integral is calculated from displacements obtained in the first step. The method is in principle thus not limited to be used jointly with just peridynamics, but allows also other methods that obtain displacements of a body, analytical, numerical or experimental. In particular, the method is equally applicable to bondbased as well as state-based peridynamics. At present, however, peridynamic J -integral works are foremost bond-based. Therefore, in this work only comparisons with bond-based peridynamic J -integral applications are made.
In the next section, we briefly introduce the bondbased peridynamic theory; for the state-based peridynamic theory, the reader is referred to Silling et al. (2007). The sections thereafter will, in order, deal with: derivation of the J -area integral on displacement formulation, derivation of the exact analytical solution of the crack problem, description of the discretization for implementation, and algorithm for calculating the J -area integral. The accuracy of the method is then investigated by comparing the peridynamic solution with the exact analytical solution of the crack problem, followed by studying conformity when the peridynamic discretization method is changed. Thereafter, the results are discussed and conclusions are given in the last sections.

Bond-based peridynamic theory
The peridynamic equation of motion of the material point at position x at time t is given as (Silling 2000;Silling and Askari 2005) where is the domain of the body, u is the displacement vector field, ρ is the mass density and b is a prescribed A material is called microelastic if the pairwise force between material points is derivable from a micropotential ω (Silling 2000): where ξ = x − x is the relative position of two material points in the reference configuration, and η = u x , t −u (x, t) = u −u is the corresponding relative displacement in the deformed configuration. A linear microelastic material results in the micropotential where s is the relative elongation of a bond: Differentiation of (3) according to (2) gives where (ξ + η)/||ξ + η|| = e, is a unit vector along a line through the two points of a bond in the deformed configuration. As we assume that a material point x does not interact with material points outside its horizon, f = 0 for ||ξ || > δ. The particular kernel of the integrand in Eq. (1), here the ratio c(||ξ ||)/||ξ ||, is common in peridynamic mechanical problems. Other kernels are possible and the selection influences the nonlocality, convergence, and thus the the discretization applied (Chen et al. 2016).

Micromodulus
The elastic stiffness of a bond is determined by the micromodulus function c(||ξ ||), which is found by calibrating the peridynamic strain energy density against the classical strain energy density, for a homogeneous body under a) isotropic (dilatational) deformation and b) pure shear (distortional) deformation. The micropotential ω is the energy in a single bond with dimension 'energy per unit volume squared'. The strain energy density of a single point is therefore The factor 1 /2 appears as the points in a bond shares the bond energy between them equally. For a 2D body where t is the thickness of the body, and c 1 comes from assuming a constant micromodulus c(||ξ ||) = c 1 . Isotropic deformation and plane stress give the classical strain energy Setting W = W 0 yields the corresponding micromodulus: The isotropic and pure shear deformations must result in the same c 1 , which in 2D plane stress restricts Poisson's ratio to 1 /3, and to 1 /4 in 2D plane strain and 3D models (Silling 2000;Gerstle et al. 2005). This is because the forces within a bond depend only on the two material points of a bond (and no other points). This restriction was overcome in state-based peridynamics by letting each bond depend on the collective deformation within the horizon (Silling et al. 2007). The forces do not necessarily have to be pairwise equal in magnitude and opposite to each other as in bond-based peridynamics. However, an advantage of bond-based peridynamics is that it is computationally less expensive. Equation (8) is derived under the assumption of a constant micromodulus for plane stress in 2D (in which case ν = 1 /3), and it holds also for plane strain (ν = 1 /4) (Gerstle et al. 2005). Other types of micromoduli (triangular, conical, higher degree polynominal) are derived in a similar fashion and are available for 1D (Bobaru et al. 2009), 2D (Ha and Bobaru 2010) and 3D (Silling and Askari 2005). Note that these micromoduli are calculated assuming a continuous body, i.e., an infinite number of points inside the horizon. For a finite number of points in this region, the micromodulus is, depending upon the type, more or less dependent upon the number of points there Stenström 2020, 2021). Depending on the desired accuracy, and with respect to the type of micromodulus and the relative grid density factor (m), the dependence of the micromodulus upon the number of points inside the horizon might be necessary to take into account.
In this study, both c derived from assuming an infinite and finite number of material points within the horizon will be applied, i.e. c ∞ (Eq. (8)) and c(m). Approximate values of c(m = 3) and c(m = 5) have been derived in Appendix A and are c m=3 = 1.1507c ∞ and c m=5 = 1.0697c ∞ . It is seen here that the approximation error decreases as the grid density is increased from 3 to 5.

Skin effect
Each material point x interacts with its neighboring points, the family of x, within a radius δ. Therefore, material points at a distance smaller than δ from a free surface will experience a truncated horizon, resulting in a smaller domain of integration. This skin effect (Ha and Bobaru 2011) leads to a softer material and larger strains near boundaries. Various correction methods exist, and have conveniently been reviewed and compared by Le and Bobaru (2018).
In this study, 'the energy method' method by Madenci and Oterkus (2014) will be used for skin effect correction; see Appendix B.

Constraint conditions and stresses
The peridynamic equation of motion is a nonlinear integro-differential equation in time and space, without spatial derivatives, and therefore does not lead to any classical boundary constraints. Instead volume constraints can be applied through a nonzero volume rather than on a surface (Silling and Askari 2005;Du et al. 2012), commonly introduced within a layer of thickness x or δ (Hu et al. 2012).
The treatment of stresses at boundaries also differ from that in classical continuum mechanics. Instead of introduced as traction forces, stresses at boundaries are introduced in peridynamics as body forces b within a layer x (Madenci and Oterkus 2014) or δ (Silling and Askari 2005). In practice, to prescribe a stress σ to a material point x, one divides with the grid spacing x, to attain force per unit volume.

Volume correction
The numerical integration of a material point x over its horizon H x includes the entire volume of each material point x within the horizon radius δ (Fig. 13). Therefore, the integration area will be larger than the area of the disc-shaped horizon, and thus, a correction factor is introduced, commonly as a linear variation between 1 and 1 /2 (Parks et al. 2010): where r = x/2, and a uniform cubic lattice is assumed. Eq. (9) will be used in this study and called LAMMPS method. Other correction methods exist and can be found in Seleson (2014) and Seleson and Littlewood (2018).

Evaluation of the peridynamic equation of motion
For numerical implementation, the integral in the peridynamic equation of motion, Eq. (1), is replaced by a summation. For dynamic or quasi-static modeling, the common implementation in practice is by time integration through a finite number of small time steps and looping over the equilibrium equations of the material points. Detailed algorithms are found in Littlewood (2015) and Madenci and Oterkus (2014), and a supportive flowchart can be found in Javili et al. (2019).
For further reading, reviews are offered by Javili et al. (2019) and Diehl et al. (2019). For peridynamic modeling of Mode I loading, which is used in the present study, see Diehl et al. (2016) for investigation on convergence and crack initiation.

The J-area integral as a function of displacement derivatives
The J -integral of a plane homogeneous body is defined as (Rice 1968) where x = x 1 , y = x 2 and u i, The body contains a straight through crack parallel to the x-axis. J is a contour integral evaluated counterclockwise along an arbitrary path enclosing the crack Fig. 2 Region A enclosed by the contour S 2 + S u − + S l with outward unit normal m i . Adapted from Li et al. (1985) tip. T i = σ i j n j is (the components of) the traction vector along , with the outward unit normal vector n j and σ i j stress. u i is a displacement vector and ds is a segment of arc length along , the inner contour in Fig. 2. 3.1 From path to area integral By using the identity dy = δ j1 n j ds and Eshelby's energy momentum tensor P jk = W δ jk − σ i j u i,k (Eshelby 1975), Eq. (10) can be written as (11) Now, we add another path S 2 outside , enclosing the crack tip and traversed counterclockwise so that S 2 , the crack surfaces S u and S l and − together enclose an area A (see Fig. 2). We now introduce a weight function q that takes unit value on and vanishes on S 2 , and is a differentiable function of position in A, including its boundaries. Then J can be written as Note that J vanishes on unloaded, stress free crack surfaces, that is, J vanishes on S u and S l , as both dy and T i vanish on S u and S l . The region enclosed by A does not contain any crack tip singularity and is therefore regular. The theorem of Gauss applied on Eq. (12) then yields as P j1, j = 0, due to absence of body forces. This is the same expression as in Li et al. (1985). Denoting the integrand in Eq. (13) I and expanding yields which is the desired result.

Weight function
Let the weight function q be a truncated square pyramid with base corners at r and s, see Figs. 3 and 4. d i = a gives a full (untruncated) pyramid, while a < d i < d y gives a truncated pyramid, and d i = d y gives a single contour, i.e. the J -integral contour. Before integration of the three regions ( 1 -3 ) of Fig. 3, we need to write the strain energy density and stresses as functions of displacement derivatives.

Displacement derivatives
For a linear elastic material, assumed in this work, the strain energy density Hooke's law provides the following plane stress versus strain relationships in 2D: The strain-displacements relationships are Substitution of Eq. (17) in Eqs. (16a-c) provides the stress-displacements derivative relationships Substitution of Eq. (18a-c) in Eq. (15) yields W as a function of displacement derivatives: where E and ν are Young's modulus and Poisson's ratio, respectively. The second term in Eq. (10) is expanded to σ i j n j u i,1 = (σ 11 n 1 +σ 12 n 2 )u 1,1 +(σ 21 n 1 +σ 22 n 2 )u 2,1 By taking the integration contour as a square around the crack tip, as in Fig. 3, Eq. (20) can be separated into three parts; 1 right hand side, 2 top and 3 left hand side contours: where w 1 , w 2 and w 3 are calculated on the right, top and left sides of the contour, respectively. Substitution of Eqs. (22a-c) in Eq. (14) yields

Integration of regions
The three regions of Fig. 3 ( 1 -3 ) are now ready for integration. The q and I for each region is given in Table 1. The J -area integral, henceforth J A , is then given by: We will now go to finding an exact analytical solution of the center cracked specimen for benchmarking, followed by discretization and the algorithm for J A .

Exact analytical solution of the center cracked specimen
Exact analytical solutions of stresses and displacements for an infinite plate with a straight crack, under uniform, uniaxial remote stress perpendicular to the crack, have been derived by Unger et al. (1983). In-plane (or Mode I) stresses and displacements can be expressed in terms of the real and imaginary parts of the complex potential of Westergaard (1939) where σ ∞ is the remotely applied tensile stress, a is half the crack length, and z = x + iy is a complex coordinate. Unger et al. derived the expressions for stresses and displacements under plane strain conditions by making use of the complex identity of Aifantis and Gerberich (1978): The real and imaginary parts of Eq. (25) can then be found, from which in turn closed form stress and displacement expressions can be derived. We have in an analogue manner derived the stresses and displacements under plane stress conditions where X , Y and A-C are dimensionless variables, α and β are constants for plane stress and plane strain, respectively. The dimensionless variables are given as 28f) α and β are given by: q is calculated using Figs. 3 and 4, and I is calculated using Eq. (23) for plane strain (29c) x and y in Eq. (28a,b) are the spatial coordinates. With Eq. (27) the stresses and displacements can be calculated at any point in the infinite plate. The benefit of this exact solution is that the stresses or displacements can be used as input to peridynamic modeling, or more precisely, as boundary conditions for finite models of the infinite geometry problem.

Discretization of the J-area integral for numerical implementation
The material domain is discretized uniformly to allow for a midpoint integration scheme (one-point Gaussian quadrature) (Silling and Askari 2005). The J -area integral of Eq. (24) can then be recast for numerical implementation by using the trapezoidal rule as follows: where: i: Material point i = 1, . . . , n j , n j : Number of material points in the x/y direction of contour j, j: Contour j = 1, . . . , N , N : Number of contours.
Following the uniformly discretized material domain, x 1 and x 2 refer to the constant grid spacing x. Moreover, d equals d y − d i , as given in Table 1. Last, W i , w i 1 and w i 2 are given by Eqs. (19) and (22a-c) The strains in Eqs. (19) and (22a-c) are approximated by applying the central difference scheme on the displacement field. For a material point i, the strains are given as follows:

Algorithm for evaluating the J-area integral
Recall that the exact solution is of an infinite specimen and the peridynamic counterpart is of a finite one. Thus, exact stresses are calculated at coordinates that corresponds to the boundary of the finite peridynamic specimen and then imposed as boundary conditions in the peridynamic model.
Within the algorithm, evaluation on the exact solution is indicated with a prime, e.g. J A . The algorithm for calculating J A is as follows: -Approximate strains at nodes x and y, using the obtained u and v (Step 3), within the interval Eqs. (30a-c) We use Fortran language for peridynamic modeling and Matlab software for pre-and post-processing. Fortran example codes for peridynamic modeling are provided by Madenci and Oterkus (2014).

Validation of J A with study of PD discretization methods
Validation of J A will be carried out by comparing it to the analytical solution J 0 = K 2 I /E, and to the Jcontour integral as they should give similar results. Thereafter, we will study the J A near and remote to the crack tip, to see the impact of the tip and skin effect. J A should be less sensitive to disturbances. Here, we will also study the influence of the 'corner points' of the rectangular integration domain. Finally, we will make changes to the discretization parameters and evaluate J A 's conformity. As J is evaluated on a contour or region enclosing the crack tip, it can be used as a tool to test the accuracy of peridynamic discretizations, which is demonstrated.

Problem set-up
To facilitate comparison with previous results, e.g. of Stenström and Eriksson (2019) and Hu et al. (2012), we will use a Young's modulus of 72 GPa, a Poisson's ratio of 1 /3, a 10 cm by 10 cm specimen, and we will impose boundary conditions within a layer x, as suggested by Hu et al. (2012). Prescribing the body force within a layer x also maximizes the skin effect.
The stresses are imposed along the boundaries at the top, bottom and right hand side of the test specimen, and symmetry conditions of zero horizontal displacement are imposed on the left hand side of the specimen (Fig. 5).

Fig. 5
The center cracked specimen with crack length a, analytical stresses at the boundary (i.e. equivalent body forces thereof), and symmetry boundary conditions To avoid peridynamic skin effects along the symmetry line, i.e. along the x-axis, the peridynamic model is taken as a square, and not as a symmetry half (Fig. 5).
For a center crack in an infinite plate with a remote loading of σ 0 = 1 MPa perpendicular to the crack, and a half crack length a = 0.05 m, J 0 = K 2 I /E = (σ √ πa) 2 /E = 2.1817 Pa m. The relative difference of the J A a on displacement formulation can then be calculated as (J A − J 0 )/J 0 . By inserting the remote loading σ 0 in the exact solution of the center crack problem, we can calculate the stresses corresponding to the boundary of the specimen shown in Fig. 5; see Fig. 6.

Peridynamic discretization
For validation of J A , we use the discretization method of the code provided by Madenci and Oterkus (2014), here called method PD1. Thereafter we will vary the grid density factor m, the micromodulus c, neighbor volume correction (partial volume, PV) and skin effect / surface correction. See Table 2.
In the table, c ∞ is given by Eq. (8), c(m) is given by Eqs. (A.5) and (A.6) in Appendix A, 'c ∞ with W 0 ' refers to the correction given by Eq. (B.2) in Appendix B, and PV refer to LAMMPS volume correction given by Eq. (9). All four methods can be complemented with skin effect correction but the default value is that it is turned off. As described in Sect. 2, the micromodulus c is commonly calculated assuming a infinite number of neighbor points within the horizon region H x , whose circumference, in turn, divides neighbor elements located at the horizon radius δ into two parts, one part inside of the horizon and one part outside. In method PD0, non of these discretization effects are corrected. In PD1, c = c ∞ is corrected by calibration with the classical strain energy density W 0 and partial volumes are corrected using LAMMPS method. In PD2, c is calculated assuming a finite number of neighbor points and calibrated with use of W 0 . The only difference between c(m) and 'c ∞ with W 0 ' is how W 0 is expressed (see Appendix B). Volumes are corrected in PD2 using LAMMPS method. PD3 differ from PD2 by an increased grid density, δ = 5 x.
PD1 will be used in the next section, and all four methods will be used in the section there after.

Validation of J A
PD1 is used for the validation. There is no skin effect correction and the 10 cm by 10 cm specimen is discretized into 50 2 or 500 2 material points.
An example of integration region (contours no. 8-11) is shown in Fig. 7. The integration domain's 'corner points' are colored orange. The first contour (at the crack tip) comprises four material points and the second twelve. Thus, the contours nearest to the crack tip for which J A can be calculated are contours no. 1 and 2. Further, Fig. 7 shows the deformed configuration of the peridynamic model with displacements magnified 3000 times and color corresponding to the the strain energy density.
For evaluating the J A , we start by plotting the Jcontour integral (J C ) from the second contour from the crack tip to the second last contour; see Figs. 8, 9.  The discretization is of 50 2 and 500 2 material points, respectively. The skin effect that is apparent in the results of the peridynamic model has two origins. The first is the skin effect at the crack faces, which effect all J C results. The second origin is the skin effect at Recalling J 0 = 2.1817, the relative error of the analytical J C at contour 14 in Fig. 8 is 0.032% and that of the corresponding peridynamic model 15%. In Fig. 9, the corresponding errors at contour 126 are 0.005% and 0.38%; the skin effect has less influence on the result. These results are in line with the convergence study presented in the previous article (Stenström and Eriksson 2019). Keeping the number of material points covered by the horizon (m) constant, while decreasing the horizon radius δ, i.e. increases the total number of material points of the model, is called δ-convergence (Bobaru et al. 2009). The skin effect can also be reduced by use of correction methods (Le and Bobaru 2018) or by using the state-based peridynamic method. In 1D, it is possible, in certain situations, to altogether eliminate the skin (i.e. end) effect Stenström 2020, 2021).
The J A gives the same result and error as J C when evaluated at the middle contours; see Fig. 10 for the 50 2 model results. d y − d i is set to 5 x for simplicity, i.e. five contours. It is seen that J A is less sensitive to local disturbances, or in this case, the skin effect.
Considering the contour corners points, it is favorable to include the corner points in the side contours instead of the top contour in the peridynamic model, as slightly more accurate J values are obtained.
As J A has been shown to give similar results as J C , we will now move on to changing the PD discretization method and study J A 's correspondence.

J A and changing PD discretization methods
For studying the discretization methods PD0-3, the methods are first applied to a defect free 10 cm by 10 cm specimen in uniaxial normal stress of 1 MPa, applied on the opposite boundaries in the y-direction. The stress is prescribed as a body force in a layer x. Material parameters are kept the same as for the center cracked specimen. Modeling a defect free specimen allows the study of the stiffness in the bulk of the specimen and at its boundary, and consequently relates it to the discretization method and J A .
Modeling results of the defect free specimen is shown in Fig. 11. The error in displacement between the methods PD0-3 are calculated as e = (u P D −u 0 ) /u 0 . PD1-3 ( Fig. 11b-e) give errors between 0-3%. PD0 (Fig. 11a) gives an error of 7-13%, which can be expected as c and volume corrections are not applied and m is the lowest advisable at 3. Increasing m reduces the error. The difference between PD2 and PD3 is that the m is increased to 5. However, both give similar errors despite the increased m, which is explained by the increased skin effect as the number of material points are kept constant at 50 2 . In PD3, the stiffness in a layer 5 x from the boundary is lowered, which propagates throughout the whole model.
Boundary nodes error in displacement is 4% for PD2 and 18% for PD3 (Fig. 11c-e). When skin effect correction is applied in PD2b and PD3b, theses values become 0.13% and 5.6% ( Fig. 11d and f). However, the error in the bulk increases. Better performing skin effect correction methods are found in Le and Bobaru (2018).
In Fig. 12, the methods PD0-PD3 are applied to the cracked specimen and J A is calculated. The PD1 curve is exact to Fig. 10, i.e. no difference in the problem setup or discretization. PD1 and PD2 give almost the same J A , which is because the correction factors for c differ only in how W 0 is expressed. PD0 gives a J A closest to the analytical curve; however, in a δ-convergence, when the number of material points increases, PD0 will not converge to J 0 = 2.1817 as the other methods. PD3 gives the highest J A which is a result of the increased skin effect and a low number of material points (50 2 ). When the skin effect procedure is applied, both curves of PD2b and PD3b are shifted down. The whole curves are shifted down as the skin effect on the crack faces are present at all contours 1-25.

Conclusion
The J -area integral is an alternative to the corresponding J -contour integral for evaluating the J -integral in peridynamics, or for any other method where the displacements are the principal unknowns, or when a displacement formulation is appealing.
The area integral was found to be less sensitive to local errors, which in this case are due to the skin effect.
Lastly, including the corner points of the integration domain in the side contours somewhat improved the accuracy of the area integral.
A possible future work is to include Mode II loading by comparing peridynamics and FEM, and/or possibly study convergence towards the stress intensity factor K I I of a slanted crack in a biaxial stress field of an infinite plate. Fig. 13 Family members of a quarter horizon region H x of radius 3 and 5 , giving N h = 4·7 and N h = 4·20, respectively

A. Approximation of c(m)
For a finite number of family members within the horizon region N h , Eq. (6) for a 2D body becomes: where v i is a volume correction factor. m = 3 gives N h = 28 family members (see Fig. 13 and thus, c 80 ≈ 1.0697c ∞ .

B. The energy method
The energy method (Madenci and Oterkus 2014) for the special case bond-based peridynamics has been derived by Le and Bobaru (2018). The purpose of the method is to multiply c of material points near boundaries with a factor to attain a stiffness closer the bulk material. The idea is to prescribe a uniaxial tension boundary condition in the x and y directions, consecutively, and for each of the two loading conditions, calculate the strain energy density of that x/y-component. Each point can then be prescribed a strain energy vector W(x) = (w x (x), w y (x)). Points near boundaries with a truncated horizon region will have different W compared to the points in the bulk, whose proportionality leads to a multiplication factor that can be made to vary with bond direction. The strain energy density of a material point is given recalling Eq. (6): Prescribing uniaxial tension boundary condition and calculating the strain energy components of x, leads to the defined multiplication factor For a point in the bulk, w x = w y = w 0 . Since a bond is shared between two points, the multiplication factor is defined as the average of the two: k = [h(x)+h(x )]/2. Between points x and x , k = (k x , k y ) and the unit vector of its bond is n = (n x , n y ). Assuming an ellipsoidal variation with direction, a scalar factor, or scaling constant, for each bond is defined as: k ave = 1 n x k x 2 + n y k y 2 (B.2) k ave will for bonds in the bulk be one and for points near boundaries be more than one. The strain energy vector W can be calculated in the same fashion as in Appendix A, together with a volume correction method.
If the micromodulus c is an approximation, then the component w 0 of points in the bulk will be an approximation. However, w 0 can be set to the classical strain energy density to both calibrate c and adjust bonds in the boundary regions. Thus Eqs. (B.1) and (B.2) will always be more than one. This method of calibrating c is applied in the Fortran codes provided by Madenci and Oterkus (2014), which is used in this study. The plane stress assumption used, of σ 11 = σ 22 , 11 = and 22 = 0, leads to W 0 = w 0 = 9 16 E 2 . For a grid density m = 3, numerically k ave ≈ 1.153, which is close the value given by Eq. (A.5). The only difference between the approximation of c(m) and the energy method in the special case bulk is the formulation of the classical continuum W 0 .