The J-contour integral in peridynamics via displacements

Peridynamics is a nonlocal formulation of solid mechanics capable of unguided modelling of crack initiation, propagation and fracture. Peridynamics is based upon integral equations, thereby avoiding spatial derivatives, which are not defined at discontinuities, such as crack surfaces. Rice’s J-contour integral is a firmly established expression in classic continuum solid mechanics, used as a fracture characterizing parameter for both linear and nonlinear elastic materials. A corresponding nonlocal J-integral has previously been derived for peridynamic modelling, which is based on the calculation of a set of displacement derivatives and force interactions associated with the contour of the integral. In this paper, we present an alternative calculation of the classical linear elastic J-integral for use in peridynamics, by writing Rice’s J-integral as a function entirely of displacement derivatives. The accuracy of the proposed J-integral on displacement formulation is investigated by applying it to the exact analytical displacement solution of an infinite specimen with a central crack and comparing the exact analytical expression of its J-integral KI2/E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_I^2/E$$\end{document}. Further comparison with a well-known peridynamic crack problem shows very good agreement. The suggested method is computationally efficient and further allows testing of the accuracy of a peridynamic model as such.


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). The J-integral has ever since been extensively used to study crack propagation conditions in both linear and nonlinear elastic materials, perhaps foremost concerning ductile materials. It is since long a firmly established and an essential supplemental tool to classic continuum solid mechanics.
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. The expression derived by Hu et al. is a special case of the more general expression by Silling and Lehoucq. The peridynamic J-integral formulation requires computation of a set of 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. Breitenfeld et al. (2014) used the classical J-integral in a study of the accuracy of 2D peridynamic crack tip stress and displacement fields. The J-integral computation scheme is however not described in detail.
In this paper, we present a calculation of the classical linear elastic J-integral for use in peridynamics, by writing Rice's J-integral as a function entirely of displacement derivatives. In peridynamic problem formulation, displacements are generally the principal unknowns to be determined, from which other quantities subsequently can be obtained. The J-integral on displacement formulation can thus be directly obtained from the displacement field of a peridynamic model. This formulation requires less computation time. It depends only on the adjacent neighboring material points of the integration path and is therefore more efficient. At the same time, it also provides a method to test the accuracy of a peridynamic model as such.
The peridynamic theory is briefly introduced in the next section. In the section thereafter, we derive the J-integral on displacement formulation. The accuracy of the method is investigated by comparison with the exact analytical solution for an infinite specimen with a central crack. 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 where is the domain of the body, u is the displacement vector field, ρ is the mass density and b is a prescribed body force field. f is the pairwise force function (a vector) per unit volume squared, denoting the force the material point at x exerts on the material point at x. This interaction between pairs of material points is called bond, or spring in case of a linear elastic material. The integral is defined over a region H x , of radius δ, called the horizon, Fig. 1. The horizon can be seen as a sphere, disk or range, for 3D-, 2D-and 1Dmodels, respectively. A suitable horizon size is chosen and the material body discretized in accordance with problem geometry, loading and desired accuracy of the results. Convergence studies may be performed to justify the selection of horizon and discretization. The relative grid density factor m = δ/ x, where x is the uniform grid spacing, should in a plane square lattice arrangement have a ratio of at least 3 (Silling and Askari 2005;Madenci and Oterkus 2014) and in many cases 4 or higher (Ha and Bobaru 2010;Henke and Shanbhag 2014;Dipasquale et al. 2016), to provide grid independent crack growth patterns.
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 mechanical problems. Other kernels are possible and the selection influences the nonlocality, convergence, and thus the discretization applied (Chen et al. 2016). 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 less computationally expensive.

The J-integral as a function of displacement derivatives
The J-integral of a plane homogeneous body is defined as (Rice 1968) where x = x 1 and y = x 2 are 2D Cartesian coordinates with origin at the crack tip. 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 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 an element of arc length along . The strain energy density is where i j is strain. For a linear elastic material, assumed in this work, Eq. (10) reduces to Hooke's law provides the following plane stress versus strain relationships in 2D: The strain-displacements relationships are Substitution of Eq. (13) in Eq. (12a-c) provides the stress-displacements derivative relationships Substitution of Eq. (14) in Eq. (11) 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. (9) is expanded to By taking the integration contour as a square around the crack tip, Fig Substitution of Eq. (14) in Eq. (17) yields where w I , w I I and w I I I are calculated on the right, top and left sides of the contour, respectively. The factor E/(1 − ν 2 ) is left out for the time being, to be invoked later. Due to symmetry, the J-integral (Eq. 9) can then be calculated as twice the sum of the three contour parts: The factor 2 on the right hand side in Eq. (19) accounts for the J-integral contour starting from the lower fracture surface, going around the crack tip counterclockwise and ending on the upper fracture surface. Equation (19) will be compared, in what follows, with the exact analytical solution of 'the central crack' problem and applied to a corresponding peridynamics model. We also note that Eq. (19) can be applied to any other method solving for displacements.

Exact analytical solution of the central 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) Z I : 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. (20) 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 α and β are given by: for plane strain (24c) x and y in Eq. (23) are the spatial coordinates. The geometry of the material test specimen is shown in Fig. 2. The dimensions are chosen to 10 cm by 5 cm in size, with a half crack length a of 5 cm, in order to facilitate comparison with previous results, e.g. Hu et al. (2012). With Eq. (22) the stresses and displacements can be calculated at any point in the infinite plate. The displacements can be used as input for calculating the Jintegral on displacement formulation, and the displacements or stresses can be used as input to peridynamic modelling, or more precisely, as boundary conditions for finite models of the infinite geometry problem. The J-integral of a central crack in an infinite plate in Mode I loading is given by The relative difference of the J-integral on displacement formulation, Eq. (19), can then be calculated as (J − J 0 )/J 0 .

Numerical implementation of the classical J-integral based on displacements
The strains in Eqs (15) and (18) can be approximated by applying the central difference scheme on the displacement field. For a material point i, the strains are given as follows: The material domain is discretized uniformly to allow for a mid-point integration scheme (one-point Gaussian quadrature) for the domain integral (Silling and Askari 2005). Furthermore, we can approximate the integrals in Eq. (19) with the trapezoidal rule as follows: W i and w i are given by Eqs. (15) and (18). n j is the number of material points along the contours I -I I I . n 1 appears due to the identity dy = n 1 ds (n 1 is the x-axis component of the outward unit normal vector of the contour) which makes the contribution of W I I vanish.
The algorithm for calculating the J-integral is as follows: The 10 cm by 5 cm specimen is discretized into xand y-coordinates. The analytical stresses and displacements are calculated for each coordinate as per Eq. (22). The analytical stresses, or displacements, along the boundary of the specimen are inputs to the peridynamic model (Fig. 3). To avoid peridynamic boundary effects along the symmetry line, i.e. along the x-axis, the peridynamic model is taken as a 10 cm by 10 cm domain. Thus, the 10 cm by 5 cm specimen half is reflected in the x-axis. The output from the peridynamic model is the displacement of material points. The J-integral on displacement formulation, Eq. (27), can now be computed on a contour enclosing the crack tip. We use Fortran language for peridynamic modelling and Matlab software for pre-and post-processing. Fortran example codes for peridynamic modelling are provided by Madenci and Oterkus (2014).
The peridynamic micromodulus of Eq. (8) was calibrated on the classical strain energy density at a point x interacts with its neighboring points, the family of x, within a radius δ. Therefore, material points that are located less than a distance δ from a free surface will have a truncated horizon, resulting in a smaller domain of integration. The effect is a softer material and larger strains near boundaries. This skin effect will be further discussed in the results section.
The peridynamic equation of motion is a nonlinear integro-differential equation in time and space, but does not contain any spatial derivatives. Instead of local constraints, nonlocal volume constrains are applied (Silling and Askari 2005;Du et al. 2012), commonly introduced within a layer x or δ (Hu et al. 2012).
Stresses at boundaries also differ from that of the 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 this work, we impose boundary conditions within a layer x, as suggested by Hu et al. (2012). 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. 3).
The numerical integration of material point x over its horizon H x includes the entire volume of each material point x within the radius δ. For a horizon of δ = 3 x, the integrated area would be larger than the area of the disc-shaped horizon. Therefore, a correction factor is used in the peridynamic code as per Madenci and Oterkus (2014).

Results and discussion
To facilitate comparison with previous results of others, e.g. Hu et al. (2012), we will use a Young's modulus of 72 GPa and a Poisson's ratio of 1 /3 for modelling. We will further keep the number of material points covered by the horizon m constant at three, i.e. δ = 3 x.
For a central crack in a plate with a remote loading of σ 0 = 1 MPa and a half crack length a = 0.05 m, the J-integral J 0 = 2.1817 Pa m, Eq. (25), is our reference value.
By setting the remote loading σ 0 to 1 MPa in the exact analytical solution of the central crack problem, we can calculate the boundary stresses of the specimen shown in Fig. 3. The boundary stresses, Fig. 4, are further used as inputs to the peridynamic model.
Resulting strain energy densities, Eq. (15), calculated with the exact analytical displacements and of the corresponding peridynamic formulation are shown in Fig. 5a-c, respectively. The results in Fig. 5a, b are for a 250 × 250 material points model. The color map shows close agreement between the results of W from exact analytical displacements and from the peridynamic model. Skin effects appear as brighter blue colored semi-continuous lines along the boundaries of Fig. 5b, c, but are too thin to be properly visualized in Fig. 5b. Figure 5c shows a 50 × 50 points model, rendering individual material points visible. The results illustrate the deformed state with a displacements magnification of 3000 times (for visualization purpose).
The strain energy densities on the contour of the Jintegral is further shown in Fig. 6. Moreover, the results of the second term of the J-integral, Eq. (18), are shown in Fig. 7. The resolution of the underlying models is 500 × 500 material points.
Common convergence studies used in peridynamics are δ-, m-and δm-convergence studies (Bobaru et al. 2009). The δ-convergence study is performed by keeping the number of material points covered by the horizon constant (m), while decreasing the horizon radius δ, i.e. increases the total number of material points of the model. In the m-convergence study, the horizon radius is kept constant while increasing m. For δmconvergence, δ is decreased and m increased, with m increasing faster than δ decreases. In this study, we use δ-convergence for comparing the classical J-integral (on displacement formulation) results, of the exact analytical displacements and the corresponding peridynamic formulation. This is carried out because the skin effect diminishes as the material points resolution increases (δ decreases). However, it should be noted that δ-convergence depends on the peridynamic kernel, micromodulus and the discretization chosen. In certain cases, peridynamic results do not δ-converge to the corresponding classical solution; the δm-convergence may be required in such case (Chen et al. 2016). Nevertheless, δ-convergence is used in this study as it is less computationally costly compared to δm-convergence. Furthermore, the triangular micromodulus has shown to improve the convergence rate compared to the constant micromodulus applied in this study (Bobaru et al. 2009). See Table 1 for figures of the δ-convergence study. Recall that the dimension of the test specimen is 10 cm by 10 cm in size when reflected in the x-axis and that δ = m x.
The results of the calculations of the classical Jintegral (a) on displacement formulation, Eq. (27), (b) for the exact analytical displacements and (c) for the corresponding peridynamic formulation, are shown in Fig. 8. The results are presented as relative difference with respect to the reference value J 0 = K 2 I /E = 2.1817 Pa m. We use the term difference instead of error as the convergence of nonlocal peridynamic models to the corresponding classical solution depends on the kernel, micromodulus and discretization, which is beyond the scope of this study.
A calculation of the J-integral on displacement formulation using the exact analytical displacements (J A ) gives a relative difference smaller than 0.1% in comparison to the reference value J 0 , for the 50 × 50 discretization scheme. The difference may be attributed to numerical errors and the approximation of strains, Eq. (26). However, approximating strains with the central difference scheme applied to the displacement field generally results in small errors.
The J-integral on displacement formulation applied to the bond-based peridynamic model approaches J 0 as the horizon is taken towards zero under δ-convergence. A horizon radius δ of 0.6 mm or smaller with the relative grid spacing of m = 3, or about 500 material points in each spatial dimension, is required for a relative difference less than a percent. This relatively fine discretization required is a consequence of the skin effect, which is present at the specimen boundaries including the crack surfaces and the crack tip. The softer material points within a distance δ from free surfaces influence their nearby material points through weaker pair-wise forces, which influence propagates throughout the model. There are several surface correction methods available for reducing the skin effect. Le and Bobaru (2018) studied some known methods and benchmarked them using the peridynamic J-integral of Hu et al. (2012) in conjunction with FEM calculations. A relative grid density factor m = 3 has been reported necessary for general crack growth modelling. A smaller m typically results in crack growth along grid rows or columns (Silling and Askari 2005). However, the choice of the m depends on the aim of the modelling. An m of just 1 may in certain cases be sufficient (Gerstle 2015), e.g. simulation of self-similar crack growth, but such a horizon radius is a very special case of the peridynamic theory, as nonlocality is removed. In many cases m = 4 or higher is required to provide grid independent crack growth patterns independent of the grid, as mentioned in Sect. 2.
Peridynamics could be viewed as an upscaling of molecular dynamics, as both methods are nonlocal, accounting for the effects of long-range forces (Seleson et al. 2009). From this particular point of view, 500 × 500 material points for a test specimen is small in relation to the billions of atoms commonly simulated in molecular dynamics. The connection between peridynamics and molecular dynamics in fracture mechan-ics has been further discussed by Madenci and Oterkus (2014) and Tong and Li (2016).
The peridynamic nonlocal J-integral studied by Hu et al. (2012), has also been studied by Oterkus (2015) and Le and Bobaru (2018). The specimen dimensions and model setup in this work are similar to the work by Hu et al. (2012) except for that Hu et al. used nonlocal J-integral and a triangular micromodulus. By comparing the results in Fig. 8 with the results of Hu et al., we find that the relative difference at δ = 3 mm is in a similar range, of 5-10%. Le and Bobaru (2018) have reduced this difference by compensating for the skin effect through certain corrections of the skin layer stiffness.
The J-integral on displacement formulation has several advantages. Besides material parameters, only displacement derivatives are required as input, let alone pairwise forces. The J-integral on displacement formulation can therefore be calculated directly from the displacement fields obtained with bond-based or statebased peridynamic models, as well, of course, with classical continuum mechanics models. In addition, the J-integral calculation is based only on the nearest neighboring material points of the contour integral, which makes the procedure computationally efficient and simple to implement.

Conclusion
In this study, we have presented an alternative calculation of Rice's contour J-integral for a linear elastic material, formulated as a function of displacement derivatives entirely, a formulation which is suitable for direct implementation in peridynamics modeling. The exact analytical solution of stresses and displacements in a central cracked specimen, derived by Unger et al. (1983) under plane strain conditions, has been extended to include plane stress conditions, and the latter case is used as a reference solution in this work. The J-integral has been calculated using the exact analytical displacements and the displacements of the corresponding peridynamic solution of the central cracked specimen problem (Fig. 8).
The J-integral calculation with the exact analytical displacements resulted in a relative difference smaller than 0.1% in comparison to the reference value J 0 = K 2 I /E, for the problem considered. For the bond-based peridynamic model, the J-integral approached J 0 as the spatial discretization was made finer. The relative difference in relation to J 0 shrank to less than one percent as the horizon was shrunk towards zero under δ-convergence (δ = 0.6 mm and m = 3). Thus, the J-integral and the strain energy part thereof, W , on displacement formulations can be used as energy density and stress intensity parameters in peridynamics modelling. Further, as shown, the J-integral on displacement formulation can as well be used for testing the accuracy of a peridynamic model as such.
Assessing the accuracy of an approximate calculation by comparison with an exact analytical solution is more accurate, and therefore preferable, than comparison with some other approximate method, e.g. FEM, whenever possible.
In testing accuracy, the exact analytical stress and displacement fields of the central cracked specimen are ideal as boundary conditions for a finite peridynamic model. The accuracy of a model is then obtained by comparison of the computational results with some reference value.
An obvious future study is to compare results of computations of the J-integral on the present displacement formulation with those of the peridynamic nonlocal J-integral derived by Silling and Lehoucq (2010) or the one by Hu et al. (2012).