Bond-based peridynamic modelling of singular and nonsingular crack-tip fields

A static meshfree implementation of the bond-based peridynamics formulation for linearly elastic solids is applied to the study of the transition from local to nonlocal behavior of the stress and displacement fields in the vicinity of a crack front and other sources of stress concentration. The long-range nature of the interactions between material points that is intrinsic to and can be modulated within peridynamics enables the smooth transition from the squareroot singular stress fields predicted by the classical (local) linear theory of elasticity, to the nonsingular fields associated with nonlocal theories. The accuracy of the peridynamics scheme and the transition from local to nonlocal behavior, which are dictated by the lattice spacing and micromodulus function, are assessed by performing an analysis of the boundary layer that surrounds the front of a two dimensional crack subjected to mode-I loading and of a cracked plate subjected to far-field tension.


Introduction
Classical continuum mechanics (CCM) has provided the foundation for models of the response of solids subjected to mechanical, thermal, electrical and other types of loading. The theoretical underpinnings provided by CCM have enabled the development of rigorous and robust analytical and computational techniques that capture complex behavior in solids and structures including large deformations and timedependent nonlinear relationships between stress and strain. However, because of the intrinsic features and requirements of its governing equations, CCM is severely limited in its ability to describe problems involving geometric discontinuities such as cracks. For example the classical (local) theory of elasticity, developed within the CCM framework, cannot account for the initiation of a crack in an originally uncracked solid, without resorting to special treatments such as the introduction of potential crack initiation lines modeled as cohesive zones that obey certain traction-separation laws [5,18,40,42]. In fact, the partial differential equations associated with boundary value problems of elasticity are not defined along the line that defines the crack front, and therefore they predict singular strain and stress fields. Classical plasticity, which is also within the framework of CCM, is yet another example of a theory that is limited in its scope. Because it does not consider internal length scales representing the underlying microstructure and the evolution of defects, it cannot account for localization and other important features of plastic deformation. The limitations of CCM have been remedied to some extent, either by introducing into the theories supplementary conditions, or by generalizing the theories to include material length scales that are capable of capturing additional phenomena. For example, in order to deal with problems involving cracks in nominally linear elastic solids and structures, it became necessary to develop an entirely new branch of mechanics called linear elastic fracture mechanics (LEFM) [4,27,29]. LEFM, which is applicable when inelastic deformation is limited to a sufficiently small region surrounding the crack front, has proved to be successful in predicting the propagation of existing cracks at known locations within a structure. However, LEFM cannot predict crack initiation, and semi-empirical conditions must be introduced to determine the magnitude of the applied load that causes crack extension and the path taken by the crack as it grows. In order to take into account the effects of microstructure on the macroscopic mechanical response of materials in regions experiencing high strain gradients, numerous mathematical models of the so-called generalized (micropolar/micromorphic) continuum were introduced in the 1960's [37], along the lines set out by the Cosserat brothers [17]. This approach, which involves internal length scales introduced into the model by involving strain gradient terms in the constitutive laws [2,38,45], can be formulated in gradient or integral forms [22]. It has been extensively applied to eliminate the physically inadmissible stress and strain singularities associated with cracks as predicted by the classical theory of elasticity. In this context, Aifantis proposed a gradient elasticity model with only one additional constant in the constitutive equation that eliminates strain singularities along dislocation lines and crack fronts [1]. Other analytic crack-tip strain and stress fields have been derived using an implicit enhanced gradient elastic theory [3,30]. While gradient-elasticity theories provide a satisfactory description of the non-local behavior, they suffer from the problem of how to interpret and satisfy the strain gradient boundary conditions. In the integral model, the stress-strain relations present additional integral terms involving an appropriately weighted strain field [22,23,45]. Eringen et al. [23] applied this idea to obtain a nonlocal elasticity solution of the Griffith crack problem that exhibited a finite stress along the crack front. However, the integral type nonlocal theory retains the problem of requiring displacement field continuity [20]. Recently the first mathematical formulation of a continuum without spatial derivatives was proposed by Silling [48]. Inspired by the original molecular theory of elasticity by Navier [41], Cauchy [15] and Poisson at the beginning of XIX century [11,44], the Peridynamic (PD) theory replaces the differential equations of CCM with integro-differential equations. Solid materials are considered as a collection of infinitesimal particles held together by infinitesimally small forces that are functions of relative displacement between two interacting particles within a material horizon of radius d. The size of the horizon changes depending on the nature of the problem. For problems that do not exhibit significant nonlocal behavior the radius is assumed to approach zero, rendering PD equivalent CCM [51]. On the other hand, for problems that are expected to have nonlocal physical behavior [35], an appropriate nonzero value of the horizon is required. In summary, the most important features of PD are that it accounts for long-range forces and nonlocal physical behavior, and that the governing equations are valid even when geometric discontinuities evolve within the solid. The PD model as originally proposed (bond-based PD, i.e. BB-PD) is a central force model, and therefore (as for the Navier and Poisson elasticity theory), the Poisson ratio is a fixed value m ¼ 1=4. However, Silling introduced a more general formulation, coined state-based PD (SB-PD), which eliminates this limitation of the BB-PD [53] and the non-ordinary state-based model (NOSB), a tool for adapting classical material models for use with peridynamics. PD has been used to model a wide variety of problems involving fracture [8,12,25,35,49], plasticity and viscoelasticy [36,39,55]. In the context of crack problems, a bond-based formulation has been used for computing the PD version of the J-Integral [28] and a NOSB peridynamic model has been proposed to approximate the stress and displacements solutions associated with LEFM [9]. However, the authors are not aware of any systematic study on the transition from local to nonlocal behavior of a crack in a nominally elastic material.
In this paper a peridynamic formulation that includes two definition of stress within a solid is used to study the transition from local to nonlocal behavior of the stress and displacement fields in the vicinity of a crack front and other sources of stress concentration. Some theoretical aspects of PD are treated and the numerical implementation of the model is discussed. In particular, an implicit BB-PD scheme is used and an analytical formulation of the stiffness matrix is derived. Convergence of the PD formulation to the results predicted by LEFM, in the limit of the horizon going to zero, is presented first. This is followed by an investigation of the non-local behavior predicted by the PD model (m-convergence [6]). Finally, an application of the model to a centrally cracked twodimensional plate subjected to uniform far-field tension is presented and discussed within the context of network materials.
2 The discrete linearized peridynamic model PD can be considered a macroscopic scale version of molecular dynamics (MD), in that the material is modeled as points whose time-dependent position are determined by integration of their equations of motion. Therefore discontinuities such as cracks evolve naturally as a consequence of the movements of these points. PD can be summarized as follows. The acceleration of any particle at X in the reference configuration at time t is found from q € uðX; tÞ À Z Hx fðu 0 À u; X 0 À XÞdV X 0 À bðX; tÞ ¼ 0 where b is the domain occupied by the body, while X 0 À X ¼ n and u 0 À u ¼ g are the relative position (i.e. reference bond) and the relative displacement between the material points X and X 0 . The body force vector is b, and f is the pairwise force function in the peridynamic bond. 1 The integral is defined over the horizon region H x of radius d (see Fig. 1). For static problems the discretized form of the (linear momentum) equilibrium equations take the form [57]: where subscripts j denote the particles within the horizon of particle i. Thus the sum in Eq. 2 is taken over all nodes j such that jX j À X i j d (i.e. neighboring particles of particle i). In the case of a microelastic PD material, a linear relationship is adopted between pairwise forces and displacements, and the previoius equation reduces to where cðjnjÞ is the micromodulus function. In this work we adopted two micromodulus functions; one is constant (cylindrical) and the other is conical. They are written as The latter exhibits faster convergence to the classical elasticity solutions (as the horizon approaches zero) [6,7]. The expression of cðjnjÞ is obtained by equating, for a homogeneous deformation, the microelastic energy density of peridynamic continuum with the strain energy density of classical linear elastic homogeneous body [26] /ðXÞ ¼ WðXÞ where s is the bond stretch which in a linearized theory is written as where g n is the component of g directed along the undeformed bond of unit vector n=jnj. For plane-stress conditions we obtain where t is the thickness and E is the Young's modulus. The computer code performs the discretization of the PD body using mesh generation tools developed for finite-element analyses. Thus the spatial domain is discretized into a set of subvolumes, each of which has a single PD material point located at its centroid. Subsequently an algorithm determines the neighboring particles (i.e. the connections) of each particle of the discretization. The linear relationship between the magnitude of the pairwise force and the bond relative elongation is taken as (see Fig. 2 We implemented a quadrature scheme in which partial neighbor intersections are also considered. This aspect adds complexity to the creation of neighbor lists (see Fig. 3) but it increases the precision of the quadrature. The results of the partial neighbor intersection computation is calculated using the value of the volume correction coefficient a [33], defined as where Dx is the grid spacing. The linear problem is solved using an implicit scheme where the stiffness matrix is defined at bond level i.e. ½K bond as where n x and n y are the components of the bond vector between the particles i and j. This expression of ½K bond can be obtained considering the definition of the macroelastic potential energy for a PD bodŷ In the case of a single bond of length jnj between two particles i and j (see Fig. 4), we can write in a discrete form the balance of the total macroelastic energy and the work U done by the external nodal forces fpg (in the hypothesis of small deformations) aŝ where c is the micromodulus of the bond, fug is the vector representing the particle displacements, DV i and aDV j are the volumes of the bond particles (DV j is reduced according to Eq. 10), while b T defined as and relates the bond particles displacements with the bond stretch in such a way that The factor 1/2 in Eq. 12 is included because the energy stored in each bond is associated equally with the two particles i and j, thus in the case of a uniform grid it is removed from Eq. 13 in order to avoid doublecounting of the bonds i À j and j À i. In this way, following Eq. 13, ½K bond is given by Following the assembly of the global structure stiffness matrix, the displacements and forces can be determined in the usual manner. It is worth noting that in BB-PD the Poisson ratio of the material assumes the value 1/3 for (2D) plane stress and 1/4 for (2D) plane strain and 3D analyses. This is a result of the pairwise central force nature of the interaction between particles. The restriction on Poisson ratio is removed in the improved state-based PD [52,53]. The state-based PD formulation eliminates the Poisson ratio restriction,  but it comes with a computational penalty, normally increasing the cost by a factor of at least two compared with the BB-PD [16]. Since we are interested in the PD modelling capabilities of nonsingular near-front fields, and since the specific value of the Poisson ratio does not influence the qualitative nature of the solutions, in this paper we use the BB-PD formulation.

Convergence schemes, local and nonlocal behavior
One can dictate whether a PD converges to either local or nonlocal material response. As shown schematically (see Fig. 5), in d-convergence the horizon size is reduced while the number of nodes is maintained constant inside the horizon region (thus increasing the overall grid density).
With m-convergence, on the other hand, the horizon size is fixed while the grid density is increased (m ¼ d=Dx is increased). In d-convergence the nonlocal solution is expected to converge to the classical, local solution. On the other hand, under m-convergence the exact nonlocal solution for the particular horizon size selected is approached (i.e. a parameter related with the internal length of the material). The value of d=Dx ¼ horizon/grid spacing ¼ 3 (i.e. the value of m) represents a good compromise between computational costs and accuracy of the solution in PD simulations involving materials with classical local behavior. Thus in this work the minimum value chosen for the horizon is d ¼ 3Dx, (m ¼ 3). If m ¼ 1, d ¼ 1Dx and the PD model reduces to a truss-type system, in which each node is connected to only its nearest neighbors [33]. The modelling of materials which exhibit evident long range effects is more delicate. It is worth noting that the stiffness matrix for a nonlocal model contains nonzero entries for each pair of nodes that interact directly, and is inherently more dense than that of a local model. Moreover, large increases in computational costs occur as the number of elements within each quadrature point's horizon increases leading to higher integration costs as well as higher linear or non-linear solver costs due to an increase in the overall size of the system. Figure 6 shows the structure of the stiffness matrix for d ¼ 3Dx and d ¼ 12Dx (where Dx is the grid spacing), for a discretized circle made of 14,400 equally spaced particles.

Peridynamic asymptotic local near crack-front solution
In this section, the ability of the PD solver to capture the near-crack front displacement and stress fields of classical elasticity (in the limit of the horizon going to zero, i.e. d-convergence) is investigated. The universal nature of the asymptotic fields near the crack front eliminate the need to model finite geometry specimens and specific loadings when studying the near-crack front region. The dominance of the asymptotic solution allows a boundary layer analysis involving only the near-front region, to which the effects of the loading and specimen geometry are transmitted by prescribing to its boundary the displacement field from this elastic solution and the associated stress intensity factors [21,43]. The performed analysis setup is shown schematically in Fig. 7, where a circular shape of diameter 2a, unit thickness, and crack of length a (note that this model represents a crack whose length is infinite compared to the size of the modeled region), is subjected to displacement boundary conditions along its entire circumference corresponding to the plane stress ''K-field'' solution (Williams [56]) where the polar coordinates (r; h) are defined at the crack-tip, l ¼ E=2ð1 þ mÞ is the shear modulus, and k ¼ ð3 À mÞ=ð1 þ mÞ for plane stress. The boundary layer analysis thus greatly reduces the simulation volume size, and allows the results to be applied to finite geometry configurations whose near-crack front regions are dominated by their own (known) stress intensity factors [21]. The strategy used in this work to model the pre-existing crack is to include in the proximity search procedure a geometric entity through which bonds cannot pass. Any bond that would pass through the crack plane is excluded from the model. In this way the crack has zero thickness which does not depends on the specific discretization of the problem, and some particles near the crack line would have a truncated sphere of interaction H x . If instead the crack is modeled simply eliminating a row of particles, the consequence is that the crack thickness depends on the discretization. In what follows, E ¼ 70 GPa, m ¼ 1=3, and the loading amplitude is specified by prescribing K I ¼ 1, with the displacement boundary conditions given by Eq. 18 applied over a boundary region along the circumference. Four values of grid spacing (Dx ¼ 1=15a; 1=30a; 1=60a; 1=120a) and two different micromodulus functions (constant and conical) have been used with a constant value of m equal to 3, as specified in the previous section. In this way, keeping m ¼ 3 fixed, the value of d was progressively decreased and set equal to d ¼ 0:2a; 0:1a; 0:05a; 0:025a, respectively. Eight different linear analyses have been performed using the PD code. The error norms (Table 1) show that even for d ¼ 0:1a the computed values of the crack opening displacement Du y (Fig. 8) and the u x and u y displacement values are in perfect agreement with the analytical values obtained with Eq. 18. In fact, for the models with d ¼ 0:1a and d ¼ 0:05a (see Fig. 9), the average difference between the computed values of the vertical displacement field and those provided by the asymptotic formulas are less than 0.35% and 0.18%, respectively. Moreover the conical micromodulus function reduces the average difference in vertical displacements with respect to the values obtained adopting a constant micromodulus function (see Table 1).
Even though the PD theory is more general than the conventional (elasticity) theory, it is possible to relate the PD bond forces and their distributions, to the notion of stress in the continuum theory [24]. According to a physically-based interpretation of stress similar to that given by Cauchy [10,13] in the early discussions on elasticity, 2 and adapting this concept to peridynamics, the PD traction vector with respect to plane p and with outward pointing unit normal n Ã at point X can be obtained by computing all the interactions between particles arranged on a line normal to the area of interest and particles beneath this surface. This concept can be expressed in integral form as where H þ X is the positive part of the horizon region H X (i.e. the family of X) and L s is the set of collinear being  Equation 19 represents the original definition of the peridynamic traction vector given by Silling [48] in the context of homogeneous deformations. Other more sophisticated notions of the traction vector are possible for collections of particles that interact through finite distances [32,34]. However, as Silling asserts, all such notions are somewhat artificial, since the idea of a stress is fundamentally tied to the idea of particles that are in direct contact with each other [48]. The formal discretized equation of the peridynamic traction vector in the case of regular grid of particles can be written as where i is the index of the particle for which the traction vector is computed, S À is the total number of collinear particles X k determined by the limit of the influence zone of the i-particle in the negative side, and Hþ is the number of particles in the positive side of the i-particle's horizon. A i denotes the area corresponding to the i-particle, as Fig. 10 shows. According to Love [34], when defining a stress measure in a system of particles, the linear dimensions of the surface of interest on p are supposed to be such that the forces between particles, whose distance apart are greater than these linear dimensions, are negligible. Hence, in our case the intersection between H X and p is a circle of radius d. The normal component of the traction vector defined above is the normal stress It is worth noting that in principle, the fact of choosing a positive side instead of a negative side of the discretized domain H X is arbitrary. Moreover due to the discrete nature of Eq. 22 and following the considerations made in [24,32] we calculate r n Ã relatively to plane p and with outward pointing unit normal n Ã at X i heuristically as the average of the two values obtained integrating only in the positive part and then in the negative part of the horizon region. For a more precise measure of the stress, the partial intersections and the partial particles volumes should also be considered in the computations. A specific subroutine which reads the displacement values as input and computes the peridynamic traction at X summing up the forces which satisfy the definition just given, has been implemented. Clearly according to this definition of stress, there is a less number of interactions between particles j and k at zones close to the problem boundaries. Particles in these zones would have a truncated sphere H X with radius d, where particles within this volume interact.
The full-field LEFM classical r yy near-tip solution is compared with the PD solution obtained with the different models adopted. Given the specific conditions of the problem here analyzed, r yy obtained with this procedure is in good agreement with the analytical formula, as the horizon size decreases. In fact Figs. 11, 12, 13 and 14 clearly show the accuracy of the proposed algorithm in computing r yy near the tip, even in this case of non-homogeneous deformation field. It is noted that since we calculate the PD stress r yy at each point X i , we refer to r yy ; h ¼ 0 as the stress An important consideration is that due to what just said and due to non-local actions, there exists a transition zone in which r yy is reduced from a maximum value ahead the crack-tip, to zero behind the tip, where no bonds can cross, of course, the crack line [19]. Given a specific mesh, the size of this transition zone depends on the non-locality of the model, thus in the case of a fixed grid spacing Dx, it increases as the m-parameter increases (refer to Fig. 15 for a graphical explanation of this concept). The differences between the solutions obtained considering a cylindrical and conical micromodulus function appear to be more evident in the stress map than in the displacement map (see for example Table 1; Fig. 13). For each discretization considered, the analytical asymptotic local near-front r yy is quite better approximated by the conical function solution than the cylindrical one, specially along h ¼ 0. In fact a deeper analysis of the results show that the differences between the analytical solution and the PD solution results, in the case of conical micromodulus function, are more influenced by the orientation h than those obtained with the constant function.
At this point it is interesting to compare the definition of stress used in this work, to that introduced in [32], which is strictly related to an original stress definition due to Saint Venant [46] and accepted later by Cauchy [32,54]. According to this definition [14], the total stress on an infinitesimal element of a plane taken within a deformed elastic body is defined as the resultant of all the actions of the molecules situated on one side of the plane upon the molecules on the other, the directions of which intersect the element under consideration. Replacing the notion of molecule with the PD particle results in a mechanistic definition that is consistent with the well known stress concept of classical continuum theory [32]. Considering an arbitrary plane p passing through X which has a normal vector n Ã and divides the family region H X into two pieces, the force which one part exerts on the other can be expressed as The line segment X 00 À X 0 given by the couple of interacting points intersects the dividing plane p at a unique point X. The line segment X 00 À X has the length f and points in the outer direction v, i.e. v Á n Ã [ 0. The line segment X À X 0 has the length t and points in the opposite direction such that X 00 ¼ X þ fv; The integration over all interacting couples ½X 00 ; X 0 can be rewritten as a surface integral over the contact plane p of a corresponding surface density, i.e. the traction vector tðX; n Ã Þ. Hence, the PD traction vector with respect to plane p, with outward pointing unit normal n Ã at point X is now defined by [32] tðX; n Ã Þ ¼ 1 2 where L denotes the unit sphere, and dx v denotes a differential solid angle on L in the direction of any unit vector v. The factor of 1 / 2 appears in Eq. 27 because the integral sums the forces on X 0 due to X 00 and those on X 00 due to X 0 [32,51]. Equation 27 can be simplified and written in discrete form as where H À is the number of particles in the negative side of the i-particle's horizon.
The summation involves only the set of bonds passing through or ending at the cross section A i from the positive side, as Fig. 16 shows. Obviously, following Eq. 27 the same operation involving this time the bonds from the negative side is also required. The sought traction vector is then given by the sum divided by two of these two values.
The stress algorithms described by Eqs. 22 and 28 lead both, to the classical LEFM solution as the horizon goes to zero (see Fig. 17). However each definition of stress is characterized by distinctive features so that different non-local solutions can be obtained, as the horizon size increases. In fact the stress definition derived from Eq. 27 is differently influenced by long-range interactions with respect to Eq. 22. Both stress definitions can be traced back to original concepts developed in the early days of elasticity, however the stress definition corresponding to Eqs. 27-28 seems to be strictly related to the classical and local continuum concept of stress. In particular, the normal stresses near the crack surfaces are almost vanishing with Eq. 28 while one observes non-vanishing normal stress near the crack surfaces using the approach derived from Eq. 22 close behind the tip, as previously discussed . Indeed, r yy obtained with Eq. 28 also experiences a transition zone in Fig. 15 Definition of the H þ X region valid for the calculus of r yy and corresponding to the first three particles siting behind of the crack-tip: Example with m ¼ 5 Fig. 16 Definition of the PD traction vector at point X i according to Eq. 28 with respect to plane p of outward pointing unit normal n Ã which the r yy stress goes from a maximum value ahead the crack-tip, to zero behind the tip. However such transition zone is not as smooth as that corresponding to Eq. 22, as Fig. 17 shows. Moreover higher maximum non-local stress along the crack plane is obtained with Eq. 28. 3 Other definitions of stress could lead to other less or more different non-local solutions. This should not be surprising because the bond-based peridynamics microelastic model is built considering only the pairwise force acting between particles as work-conjugate of the corresponding bond stretch of the ligament. Thus the equilibrium Eq. 1 together with a constitutive f-s relationship essentially contain the totality of the peridynamic model for linear and nonlinear materials [50] and no concept of stress and strain are strictly needed or directly used within the theory.

Peridynamic nonlocal near-front solution
For an assigned grid spacing Dx, an increase in both the horizon size and the value of m results in the stress field changing from the values predicted by the singular solution to finite (non-singular) values consistent with the previously mentioned non-local solutions. This aspect is emphasized by considering different values of m (note that d is not held fixed so at the same time also the value of the horizon is increased) are adopted for a given discretization Dx ¼ 1=60a. The influence of the micromodulus function adopted (cylindrical and conical) is also investigated. Since we do not have the exact solution for the nonlocal problem, we compare the results with those of the classical solution for illustration purposes. It can be noted that the non-local stress solutions shown in this paragraph are consequence of the specific heuristic stress computation algorithm used in this work and which is applied to the study of the transition from local to nonlocal behavior of the r yy stress field. As seen in Fig. 18, as the horizon increases (and thus the m value), the difference between the PD model and the classical singular model of the stress component r yy increases. This effect is less evident but still present for the case of conical influence function. In fact, the computed magnitude of stress ahead of the crack tip changes from a certain value to zero within a very short distance. In the case of d ¼ 0:2a and m ¼ 12 a strong non-local behavior appears. The estimated stress has a maximum at some distance in-front of the tip and then drops towards a lower magnitude at the tip, a feature that is in accordance with other non-local (gradient and integral) theories and experiments in network materials (see Fig. 19 for details). This particular behavior has been widely observed in network materials and implies that the physical stress field in heterogeneous fiber-based materials is finite and experiences nonlocality: the fibers introduce longranging actions that limit the size of gradients in strain and stress around defects and cracks [30]. Differences are also evident in the near-tip horizontal displacement field between the two extreme cases of d ¼ 3 with conical micromodulus function and m ¼ 12 with cylindrical function (see Fig. 20). Moreover, the results in the sequences shown in Figs. 21 and 22 suggest that as the level of nonlocality increases, the shape of the crack-front becomes more and more   blunt. This type of blunting is a phenomenon observed experimentally in microstructured materials such as non-woven felts and paper. It is caused by bending deformation of fiber segments in the vicinity of the crack-front [30]. In a bond-based PD model instead, this behavior is due to a macroscopic bending deformation of the horizon region corresponding to particles siting near the crack-tip. The blunting effect may be observed visually: when loaded, blunted cracks have rectangular shapes rather than sharp elliptical as expected in linearly elastic fracture mechanics. PD can easily model this important feature of such type of materials because of its implicit definition of an internal length related to the horizon d. However at this point it is important to know what is the value of m for an assigned d and crack length a such that the non-local PD solution is well approximated numerically for that situation (i.e.   intensity), depending on the specific horizon-crack length ratio considered. The effect of increasing d=l keeping m ¼ 12 fixed is twofold: the maximum stress near the tip is reduced and the crack surface profile becomes more blunt (see Figs. 25,26,27). These procedures could be applied in future works for modelling the size effect in non-woven felts, paper, fiber composites, textiles and their alike, i.e. random fiber networks. In fact, such materials are heterogeneous on multiple scales [47] and especially in regions of high-strain gradient show a strong non-local behavior. When a macroscopic crack is present, fibers introduce long-range microstructural effects and distribute stresses near the crack-front, thus reducing the probability of failure immediately ahead of the crack since forces are transferred to remote regions [31]. Moreover a characteristic length (i.e. the average fiber length as demonstrated by [30]) could be related to the PD horizon d and to a specific micromodulus function calibrated on experimental results (see Fig. 28).

Conclusions
In this study a meshfree static bond-based PD formulation was applied to conduct an in-depth evaluation of the numerical method's capability to capture the local and non-local stress and displacement fields near the tip. The influence of grid spacing, micromodulus function and horizon size on the numerical results was investigated, and the overall accuracy of the method in approximating the local LEFM solution was demonstrated for both the crack opening displacement and the near-tip r yy stress and displacement fields. Given the specific conditions of the problem here analyzed, the r yy map obtained using the heuristic approach proposed seem to be in good agreement with the analytical formula, as the horizon size decreases. This is supported by numerous plots and figures which clearly show the accuracy of the proposed algorithm in computing the stress intensity near the tip even in this case of quite complex deformation field. Then, through an m-convergence study (i.e choosing an appropriate value of d=Dx), the non-local solution of the problem was investigated. The maximum non-singular hoop stress in the near-tip  region is substantially well approximated by all the discretization adopted, whereas the corresponding values obtained adopting a conical weight function result to be higher than those computed using a cylindrical function. The distance from the crack-tip to the point of maximum is roughly 0.25 times the parameter d, which is associated with a characteristic length and controls the range of nonlocal actions. Results demonstrated that m ¼ 12 is a good compromise between accuracy of the non-local peridynamic solution for a given horizon d and computational efficiency. Moreover a comparison between the stress computation algorithm used in this work and a PD stress definition strictly related to the familiar concept of stress of elasticity was made. The two stress definitions both lead to the asymptotic LEFM solution as the horizon goes to zero, of course. However, being each definition of stress characterized by distinctive features, quite different non-local solutions, as the horizon size increases, are obtained. Finally, the capability of the PD solver in modelling the transition from local to nonlocal behavior of the stress and displacement fields near the tip accounting for size effect has been demonstrated by benchmarking it with a classical 2D Griffith problem. The effect of increasing d/crack-length, keeping m ¼ 12 fixed is twofold: The maximum stress near the tip is reduced and the crack surface profile became more blunt, a well known experimental evidence in heterogeneous materials like random fiber networks.