A new method based on Taylor expansion and nearest-node strategy to impose Dirichlet and Neumann boundary conditions in ordinary state-based Peridynamics

Peridynamics is a non-local continuum theory which is able to model discontinuities in the displacement field, such as crack initiation and propagation in solid bodies. However, the non-local nature of the theory generates an undesired stiffness fluctuation near the boundary of the bodies, phenomenon known as “surface effect”. Moreover, a standard method to impose the boundary conditions in a non-local model is not currently available. We analyze the entity of the surface effect in ordinary state-based peridynamics by employing an innovative numerical algorithm to compute the peridynamic stress tensor. In order to mitigate the surface effect and impose Dirichlet and Neumann boundary conditions in a peridynamic way, we introduce a layer of fictitious nodes around the body, the displacements of which are determined by multiple Taylor series expansions based on the nearest-node strategy. Several numerical examples are presented to demonstrate the effectiveness and accuracy of the proposed method.


Introduction
The propagation of cracks in solids and structures is one of the most common problems in structural engineering. In recent years, a new non-local continuum theory able to simulate crack propagation, named peridynamic theory, attracted the attention of many researchers. Each point in a body modelled with peridynamics interacts with all the neighboring points within a distance 훿. The non-locality of the peridynamic theory is essential to describe fracture phenomena in solid bodies without ad hoc criteria. Firstly, the so-called "bond-based peridynamics" was developed [40], which however has a limited capability of prescribing the Poisson's ratio. This shortcoming is avoided by the second formulation of the theory, named "state-based peridynamics" [42].
The non-local nature of the theory leads to two interrelated problems near the boundary of the body: the "surface effect" and the difficulty to impose the boundary conditions [9]. The surface effect, sometimes also called "skin effect", is due to the fact that peridynamic points near the boundary lack some neighboring points, leading to an undesired variation of the stiffness properties in the most external layer of the body [16,19]. Bond-based and state-based peridynamic models exhibit respectively a softening and a hardening-softening behavior near the boundary [2,20,35].
Imposing boundary conditions in a peridynamic model is not a trivial task to accomplish. The application of the boundary conditions to the points on the boundary, as one would do in a local model, leads to large fluctuations of the solution near the boundary [16]. In [21] it is suggested that external loads and constraints should be imposed on a layer of finite thickness respectively inside and outside the body. This strategy is surely closer to a non-local concept, but it is not really clear the proper procedure to "distribute" the boundary conditions over the finite layers. In the following, we present some of the most commonly used methods to mitigate the surface effect and impose the boundary conditions in a peridynamic model. A possible approach is to couple peridynamics with classical continuum mechanics: peridynamics is employed only in the interior of the body and the layer of material near the boundary is modelled, for instance, with the Finite Element Method [11,25,38,45,47,48], with the Carrera Unified Formulation [30], with the Extended Finite Element Method [13,14] or with the Meshless Local Exponential Basis Functions method [39]. In this way, the surface effect and the problem of the imposition of the boundary conditions in peridynamics is avoided. However, if cracks initiate or propagate near the boundary, those regions must inevitably be modelled with peridynamics and the coupling approach is not suitable to avoid the boundary issues. Furthermore, there are some spurious effects at the interface of the coupling region due to the different formulations of peridynamics and classical continuum mechanics (see the computation of outof-balance forces in [26]).
The maximum distance of interaction, namely 훿, is a measure of the non-locality of the theory. Therefore, the external layer of the body which is affected by the surface effect is thinner as 훿 approaches 0. Similarly, the imposition of the boundary conditions in a local way (constraints and loads applied only to the points closest to the boundary) becomes a better approximation if 훿 tends to 0. Since the number of nodes is bound to increase as 훿 decreases, the computational effort may become excessive. In this case, the variable horizon method can be employed to decrease the value of 훿 near the boundaries [2,3,6,[31][32][33]44]. However, this approach of reducing the non-local nature of the peridynamic theory is solely capable of confining the solution fluctuation in a smaller region, but never of completely correcting it.
The approach of modifying the stiffness properties of the bonds near the boundary has been proposed in many methods: the force normalization method [20], the force density method [15], the energy method [21,27], the volume method [1] and the position-aware linear solid constitutive model [23]. The comparison of these methods, carried out in [16], highlights that there are still some residual fluctuations of the solution near the boundary because they do not cope with the problem of the imposition of the boundary conditions in a non-local way. Another recently devised approach consists in modifying the peridynamic formulation in points which are affected by the surface effect in order to recover the classical mechanics solution for 훿 → 0 [4,46]. Nevertheless, the treatment of the boundaries becomes much more complex.
The method of the "fictitious nodes" consists in adding around the body some nodes which provide the previously lacking interactions near the boundary, mitigating in this way the surface effect [12,15,34]. The fictitious nodes have been employed also to impose the boundary conditions: the displacements of the fictitious nodes are extrapolated by means of various types of functions, such as constant, linear, polynomial, sinusoidal or odd functions, in order to obtain the desired value of the constraint or load [6][7][8]16,21,22,28,29,31,49]. Moreover, the displacements of the fictitious nodes can be determined also by means of the formulae of classical continuum mechanics to enforce the desired load at the boundary [16,22,28,31,49]. However, these procedures to impose the boundary conditions are casedependent and are applicable only for simple geometries and boundary conditions.
We proposed a new version of the "Taylor-based extrapolation method" adopting the nearest-node strategy [35]: the displacements of the fictitious nodes are determined as functions of the displacements of their closest real nodes by means of multiple Taylor series expansions truncated at a general order 푛 푚푎 푥 . The surface effect is sensibly reduced by this effective method. Moreover, the boundary of the body is discretized by a new type of nodes, named "boundary nodes". As the fictitious nodes, the boundary nodes do not constitute new degrees of freedom in the model because their displacements are determined by means of the Taylor-based extrapolation method. Dirichlet boundary conditions are included in the Taylor series expansion of the displacements of the boundary nodes about their closest real nodes, whereas Neumann boundary conditions are imposed through the peridynamic concept of force flux.
The paper is organized as follows: Sect. 2 presents a brief review of the ordinary state-based peridynamic theory, particularly focusing on the peridynamic stress tensor, the force flux, the surface effect and the imposition of boundary conditions; Sect. 3 illustrates the Taylor-based extrapolation method and the imposition of boundary conditions in a peridynamic model; Sect. 4 shows the discretization of the peridynamic model, the numerical evaluation of the peridynamic stress tensor and of the force flux, and the numerical implementation of the proposed method; Sect. 5 compares the numerical results of several meaningful 2-dimensional examples obtained without any corrections at the boundary and by using the proposed method; Sect. 6 shows the differences that may arise in crack propagation near the boundaries between corrected and uncorrected models; Sect. 7 draws the conclusions.

Review of peridynamic theory
Peridynamic points interact with each other, even within finite distance, through entities named "bonds". A bond is identified by the relative position vector in the reference configuration as where x and x ′ are the position vectors of two points in a body B modelled with peridynamics. The bond vanishes if 1 Body modelled with ordinary state-based peridynamics in the reference configuration B 푟 and deformed configuration B 푑 : a pairwise force density f arises in the bond due to the deformation of the body the distance between the interacting points exceeds the value 훿, called "horizon". A point x therefore interacts with all the points x ′ inside its neighborhood, which is defined as where B 푟 is the body in the reference configuration. Point x is named "source point" and the points within H x are named "family points". In the deformed body configuration B 푑 at time 푡, the relative displacement vector is defined as where u is the displacement field. Note that + is the relative position of points x and x ′ in the deformed configuration. The peridynamic equation of motion of a point x within the body B is given by [40,42]: where 휌 is the material density, u is the acceleration field, f is the pairwise force density, d푉 x ′ is the differential volume of a point x ′ within the neighborhood H x and b is the external force density field. The pairwise force density represents the force (per unit volume squared) in a bond. The peridynamic equilibrium equation is derived from Eq. 4 by dropping the dependence on time: where b x = b(x). f (x, x ′ ) is the force density applied to point x due to the interaction with a point x ′ inside its neighborhood. Conversely, point x belongs to the neighborhood H x ′ , thus a force density f (x ′ , x) = −f (x, x ′ ) is applied to point x ′ (see Fig. 1). The formulae to compute the pairwise force density depending on the deformation of the body are shown in the following section.

Ordinary state-based peridynamics
In state-based peridynamics, the pairwise force density is defined as [42] f (x, where T is the force density vector state. T[x] and T[x ′ ] − depend respectively on points x and x ′ , and they respectively operate on bonds and − .
In an ordinary peridynamic material, the force density vector state is aligned with the corresponding bond for any deformation, as depicted in Fig. 1, and it can be written as where 푡 is the force density scalar state (magnitude of T) and M is the deformed direction vector state (unit vector in the direction of T), defined as Note that M = −M − . Furthermore, under the assumption of small deformation ( ≪ ), the deformed direction vector state can be approximated with the bond direction unit vector in the reference configuration: Therefore, the pairwise force density can be rewritten as The reference position scalar state 푥, representing the bond length in the reference configuration, and the extension scalar state 푒, describing the elongation (or contraction) of the bond in the deformed body configuration, are respectively defined as The influence of the neighborhood H x on a source point x is expressed by two non-local properties of that point, the weighted volume 푚 and the dilatation 휃, which are defined as where 휔 is a prescribed spherical influence function and 푐 휃 is a peridynamic constant. We adopt the Gaussian influence function since it assures a smooth convergence of the numerical integration [36]. The weighted volume describes the "fullness" of the neighborhood: a neighborhood completely full of peridynamic points results in the maximum value of 푚, whereas the weighted volume of an incomplete neighborhood has a lower value. This lack of neighboring points is the origin of stiffness fluctuations, the so-called "surface effect" [16], near the boundary of the body. The surface effect is further analyzed in Sect. 2.4.
On the other hand, the dilatation represents the volumetric deformation of the neighborhood. Consider a point x subjected to a homogeneous, isotropic and small deformation 휀, so that 푒 = 휀 푥 for each bond. The peridynamic dilatation 휃 of point x corresponds to the dilatation 휃 푐푙 in classical continuum mechanics under the same deformation if the constant 푐 휃 is chosen as [17,35,42] where 휈 is the Poisson's ratio.
The force density scalar state can be computed as [35] 푡 where 푘 휃 and 푘 푒 are the peridynamic stiffness constants. These constants are derived by equalizing the peridynamic strain energy density in a point x with a complete neighborhood under homogeneous deformation, with the classical continuum mechanics strain energy density in a point sub-jected to the same deformation [17,35,42]: where 퐸 is the Young's modulus. By substituting Eq. 17 in Eq. 10, the pairwise force density is computed as Note that the magnitude of the pairwise force density in ordinary state-based peridynamics depends on the neighborhood properties (푚 and 휃) of the two points x and x ′ connected by the bond. Hence, the resultant of all the bond forces in a point x, obtained with the integral of the peridynamic equilibrium equation (Eq. 5), depends on the deformation of the points within a 2훿-distance from x.

Peridynamic stress tensor
The peridynamic stress tensor, introduced in [18], is defined in a point x with a complete neighborhood as where 훺 is a unit sphere centered in x and d훺 m is the differential solid angle on 훺 in any bond direction m. The points x − 푠m and x + (푟 − 푠)m are connected by a bond passing through point x, and we respectively name them x ′ and x ′′ . Therefore, 푠 = x ′ − x and 푟 = x ′′ − x ′ , as shown in Fig. 2. 푠 is the distance between points x and x ′ , whereas 푟 is the length of the bond between x ′ and x ′′ . The definition of the integration domain allows to take into account all the bonds passing through point x. Note that each bond passing through x (between x ′ and x ′′ ) has a corresponding bond in the opposite direction (between x ′′ and x ′ ), so that the same pairwise force density is integrated twice in Eq. 21. This is the reason why the factor 1/2 appears at the beginning of the formula. The integral over the unit sphere 훺 is not affected by the variables 푠 and 푟, but it depends only on the bond direction m. On the other hand, the integrals related to d푠 and d푟 are interdependent, as shown by the integration domain depicted in Fig. 3. For later use, the peridynamic stress tensor in a point x with a complete neighborhood is rewritten by changing the order of the integrals: where d푉 x ′′ = 푟 2 d푟 d훺 m . Under the assumption of homogeneous deformation, the bonds with the same length and direction have the same pairwise force density in any position of the body. This means that, for each bond of length 푟 and direction m, its pairwise force density does not depend on 푠 anymore. Therefore, x can be simplified from Eq. 22 as follows: Note that, since the value of 푠 does not affect f (x ′ , x ′′ ) in a body under homogeneous deformation, we conveniently choose the pairwise force density for 푠 = 0, i.e., f (x, x ′′ ). We want to compare the peridynamic stress tensor with the stress tensor in classical continuum mechanics for the same deformation conditions. For simplicity sake, we choose the where 휀 푖푠표 and 휀 푠ℎ are the values of the imposed deformations and 휎 푖푠표 and 휎 푠ℎ are the corresponding stresses.
In the following analysis of the peridynamic stress tensor, only points with complete neighborhoods are considered. The inclination of a bond with respect to the 푥-axis is called 휙. Therefore, the bond direction in a 2-dimensional model can be written as m = {cos 휙, sin 휙} ⊤ . Furthermore, the weighted volume of a point x with a complete neighborhood is given from Eq. 13 by where ℎ is the thickness of the 2-dimensional body.
In the case of a body subjected to a small isotropic deformation, any extension scalar state is 푒 푖푠표 = 휀 푖푠표 푥 and the corresponding dilatation in a point x with a complete neighborhood is 휃 푖푠표 x = 푐 휃 휀 푖푠표 . The peridynamic stress tensor under this condition is given from Eq. 23 as The obtained peridynamic stress tensor yields the same result of the stress tensor computed with classical continuum mechanics in a point under isotropic deformation. Note that only a tensile stress 휏 11 = 휏 22 = 휎 푖푠표 arises from the imposed deformation 휀 푖푠표 and there is no shear stress (휏 12 = 0). In the case of a body subjected to a small shear deformation 휀 푠ℎ , the extension scalar state can be computed by substituting = { cos 휙, sin 휙} ⊤ and = {휀 푠ℎ sin 휙, 휀 푠ℎ cos 휙} ⊤ in Eq. 12: where the formula is simplified under the assumption of sufficiently small deformation by dropping the second order terms and employing the Taylor series expansion for the square root. The corresponding dilatation in a point x with a complete neighborhood is 휃 푠ℎ x = 0 given the anti-symmetry of the integrand and the symmetry of the integration domain. The peridynamic stress tensor under this condition is given from Eq. 23 as The obtained peridynamic stress tensor yields the same result of the stress tensor computed with classical continuum mechanics in a point under simple shear deformation. Note that only a shear stress 휏 12 = 휎 푠ℎ arises from the imposed deformation 휀 푠ℎ and there is no tensile stress (휏 11 = 휏 22 = 0). We showed that the peridynamic solution for the stress tensor corresponds to that of the classical continuum mechanics for homogeneous and small deformations, as shown also in [26] for bond-based peridynamic models. However, this statement is not valid near the boundaries of the body due to the surface effect.

Force flux
The force flux (x, n) at point x in the direction of the unit vector n (see Fig. 4) is derived from Eq. 21 as [18]: where x ′ = x − 푠m and x ′′ = x + (푟 − 푠)m. As in the definition of the peridynamic stress tensor, a factor 1/2 is required since the integration domain takes into account the magnitude of the pairwise force density of each bond twice (for both direction m and −m). We briefly recall the mechanical interpretation of the force flux [18]. Consider a plane P with normal n passing through point x, as shown in Fig. 5. Points x ′ and x ′′ respectively lie in the different half-spaces generated by plane P. The differential volumes of those points are d푉 x ′ = 푟 2 d푠 d훺 m and d푉 x ′′ = 푟 2 d푟 d훺 m . The differential area of point x ′ , which is perpendicular to the bond direction m, is the portion of a sphere centered in x ′′ with a radius 푟 which subtends the differential solid angle d훺 m , namely d퐴 x ′ = 푟 2 d훺 m . By the same token, the differential area d퐴 x ′′ on a sphere centered in x ′ with a radius 푟 is equal to d퐴 x ′ . As shown in Fig. 5, the The differential pairwise force acting through the bond between points x ′ and x ′′ is f (x ′ , x ′′ ) d푉 x ′ d푉 x ′′ . Therefore, the differential pairwise force per unit area on plane P is given by Note that the integrand in Eq. 30 is the pairwise force per unit area on plane P. This provides the mechanical interpretation of the force flux as the resultant of the pairwise forces per unit area of all the bonds intersecting P in x.

Surface effect
The non-local formulation of the peridynamic theory exhibits some issues near the boundaries due to the incomplete neighborhoods of points close to free surfaces. The peridynamic constants 푘 휃 and 푘 푒 in Eqs. 18 and 19 are derived for points with a complete neighborhood. Therefore, the points near the boundaries, whose neighborhood is lacking some bonds, have different stiffness properties with respect to the points in the bulk. This phenomenon is called "surface effect" [16].
x m x 5 Differential variables involved in the computation of the force flux (x, n): the differential volumes of points x ′ and x ′′ are respectively is the projection of d퐴 x ′ and d퐴 x ′′ on the plane P In ordinary state-based peridynamics, there are two nonlocal properties of a point which may contribute to the surface effect: the weighted volume 푚 and the dilatation 휃. The latter (Eq. 14) is independent from the neighborhood "fullness" because it is normalized by the value of the weighted volume. Therefore, we focus on the value of 푚. We define the value 푑 푏 as the minimum distance of a peridynamic point from any boundary of the body. The weighted volume has its maximum value when 푑 푏 ≥ 훿, and it decreases gradually from points with 푑 푏 = 훿 towards points with 푑 푏 = 0 on the boundary. Moreover, points approaching corners, with respect to those approaching edges or surfaces, exhibit a steeper reduction in the weighted volume and a lower minimum value at the boundary.
The equilibrium of a peridynamic point x (Eq. 5) is determined by the sum of the pairwise forces of all its bonds. Therefore, x primarily interacts with the points inside its neighborhood H x . However, the magnitude of the pairwise force density (Eq. 20) depends on the weighted volumes and dilatations of both point x and point x ′ within H x . This means that x secondarily interacts with points up to a distance of 2훿 from itself. Thus, as shown in Fig. 6, we can discriminate 3 types of points depending on 푑 푏 : 6 Types of state-based peridynamic points depending on the distance 푑 푏 from the closest boundary: points are of type-I, also named points in the "bulk", if The source point x interacts primarily with the family points within the neighborhood H x (dashed line) and secondarily with all the points in the neighborhoods of the family points (dotted line) Type-I points are said to be in the "bulk" of the body and they are the only ones which are not affected by the surface effect.
As the weighted volume of one or both the points of a bond decreases, the pairwise force density of that bond increases according to Eq. 20. As shown in Fig. 6, type-II points interact with at least one point with a partial neighborhood, so that the peridynamic forces applied to those points increase. Therefore, in the layer of the body where 훿 ≤ 푑 푏 < 2훿 the peridynamic material is stiffer and exhibits a hardening behavior. The pairwise forces applied to type-III points increase even more. However, a type-III point is affected by less bonds than type-I or type-II points due to the lack of at least one family point. Therefore, the most external layer of the body (푑 푏 < 훿) exhibits a hardening-softening behavior towards the boundary. This hardening-softening behavior can also be observed in the analytical solution of a 1-dimensional state-based body subjected to a homogeneous small deformation [35]. However, we expect that the stiffness fluctuation would be amplified near the corners of the body because the points in those regions have the smallest partial neighborhood. Figures 7 and 8 show the components of the peridynamic stress tensor x in a 2-dimensional body subjected respectively to a isotropic deformation 휀 푖푠표 and to a simple shear deformation 휀 푠ℎ . x is computed numerically with a relatively high density of nodes within each neighborhood (푚 = 10) and normalized with the analytical solutions derived in Sect. 2.2. Please refer to Sect. 4.3 for the numerical procedure to compute the peridynamic stress tensor. The numerical result for the points in the bulk of the body is really close to the analytical solution, whereas there are large differences for the points near the boundary, especially near Components of the peridynamic stress tensor x for every point in a 2-dimensional body subjected to a isotropic deformation 휀 푖푠표 . The plots are normalized with the analytical solution of the tensile stress 휎 푖푠표 for a peridynamic point with a complete neighborhood the corners. Moreover, the points near the corners, due to the asymmetry of their neighborhood with respect to both 푥and 푦-axis, have a non-zero value of the peridynamic stress even without the corresponding deformation: 휏 12 ≠ 0 in the case of isotropic deformation 휀 푖푠표 and 휏 11 = 휏 22 ≠ 0 in the case of simple shear deformation 휀 푠ℎ .

Imposition of the boundary conditions
Another issue in peridynamics, which is related to the surface effect, is the proper definition of the boundary conditions. The easiest method to impose the peridynamic boundary conditions would be to assign the desired value to the boundary points, as in classical continuum mechanics. However, this method does not consider the non-local nature of the theory and results in additional fluctuations of the solution near the application of the boundary conditions.
A widely used method suggests that external loads and constraints should be imposed on a layer of finite thickness The plots are normalized with the analytical solution of the shear stress 휎 푠ℎ for a peridynamic point with a complete neighborhood respectively inside and outside the body [21]. The finite thickness is defined to be 2훿 in state-based peridynamics [37]. Since this method involves type-II and type-III points in the boundary conditions, it is undoubtedly more accurate than the previous one. However, it is not really clear the exact procedure of "distributing" the boundary conditions over the finite layer.
We propose in the next section a novel method capable of reducing considerably the surface effect and of imposing the boundary conditions in a "peridynamic way".

Taylor-based extrapolation method
A fictitious layer F of thickness 훿 is added around the body B [12], as shown in Fig. 9. The neighborhoods of the family points of type-II points are completed thanks to the additional fictitious points, so that type-II points can be considered as points in the bulk (type-I points). Similarly, the neighbor- Fig. 9 Types of state-based peridynamic points depending on the distance from the closest boundary 푑 푏 in a body with a fictitious layer of thickness 훿: there is no difference between type-I and type-II points anymore, whereas type-III points lack some secondary interactions in the neighborhoods of the family points hoods of type-III points are completed by the fictitious layer, but some of the neighborhoods of their family points are not. However, we assign to the fictitious points the value of the full weighted volume. In this way, all the points inside the body B behave as points in the bulk. The next section shows a procedure to evaluate the displacement and dilatation fields over the fictitious layer.

Extrapolation procedure to mitigate the surface effect
The displacements and the dilatations of the fictitious points are determined by means of the Taylor-based extrapolation method [35]. Consider the displacement We name u 푏 = u(x 푏 ) the displacement of the boundary point with the minimum distance from that fictitious point (nearest-point strategy). The Taylor series expansion of u 푓 about x 푏 = {푥 푏 , 푦 푏 , 푧 푏 } ⊤ truncated at the maximum order 푛 푚푎 푥 ≥ 1 is given by where 푛 is the global order (푛 = 1 is related to the gradient, 푛 = 2 to the Hessian matrix, etc.) and 푛 1 , 푛 2 and 푛 3 are the orders respectively in 푥, 푦 and 푧 directions. Similarly, consider the dilatation 휃 푓 = 휃 (x 푓 ) of a fictitious point and the dilatation 휃 푏 = 휃 (x 푏 ) of the boundary point closest to x 푓 . The Taylor series expansion of 휃 푓 about x 푏 truncated at the maximum order 푛 푚푎 푥 − 1 is given by Note that the dilatation is a measure of the strain, thus its truncation of the Taylor expansion occurs with 1 order less than that of the displacement.
This method allows to determine the displacement field and the dilatation field in the fictitious layer F as a function of the respective fields in the body B. Since the displacement and dilatation fields in F are approximated by means of a Taylor series expansion, more accurate results are obtained by increasing the truncation order 푛 max or by reducing the thickness of the fictitious layer.
The new bonds between real and fictitious points, called "fictitious bonds", are the interactions that are lacking in the peridynamic models without fictitious layer. The Taylorbased extrapolation method provides the displacement and dilatation values of the fictitious points, which are required to compute the pairwise forces of the fictitious bonds. In this way, the proposed method is able to mitigate the surface effect.

Peridynamic boundary conditions
We propose hereinafter a novel method to impose the boundary conditions in a peridynamic way when using the previously described fictitious layer method [35]. The desired boundary conditions are applied solely on the boundary points, exactly as in classical continuum mechanics. However, the influence of the boundary conditions on the body is non-local thanks to the Taylor-based extrapolation method. This concept is explained for Dirichlet and Neumann boundary conditions in the following.
A constraint u imposed in a boundary point x 푏 is simply given as This boundary condition determines the displacement field in the fictitious layer through the Taylor-based extrapolation method (by substituting Eq. 35 in Eq. 33). Therefore, the influence of the constraint can be seen as distributed in the whole thickness of the fictitious layer, as suggested in [21, pp. 29-30]. An external load per unit area p applied to a boundary point x 푏 is expressed by means of the peridynamic concept of force flux (see Eq. 30): where n is the unit vector perpendicular to the boundary in x 푏 . By definition of force flux, (x 푏 , n) is the sum of the pairwise forces (per unit area) of all the bonds passing through x 푏 . Since point x 푏 lies on the boundary, all the bonds involved in Eq. 36 are fictitious bonds. On the one hand, the pairwise forces of those bonds applied to the fictitious points, which do not constitute new degrees of freedom, are ignored. On the other hand, the corresponding pairwise forces applied to the real points are the only ones "perceived" by the body. For how the magnitude of those forces is computed (see Eq. 20), the boundary condition in point x 푏 affects the displacement in a sphere of radius 2훿 centered in x 푏 . Therefore, the external load, expressed by means of the definition of the force flux, is distributed on the points in a layer of thickness 2훿 within the body, as suggested in [21, pp. 30-32].
The proposed method for imposing the boundary conditions, which makes use of the Taylor-based extrapolation on the fictitious layer, defines a peridynamic way to distribute the constrains or the loads in the non-local region near the boundary.

Remark 1
The reaction force acting on a boundary point x 푏 , due to a constraint imposed as in Eq. 35, can be computed as the force flux in x 푏 in the direction of the unit vector n perpendicular to the boundary, i.e., (x 푏 , n).

Remark 2
The zero-traction boundary condition is sometimes applied by removing the fictitious layer in literature [22]. However, in order to maintain the correction of the surface effect, we suggest to keep the fictitious layer and impose the condition (x 푏 , n) = 0 to all the points of that boundary.

Numerical implementation
In order to discretize the domain, a mesh-free method is adopted [36,41]. For simplicity sake, the peridynamic grid consists of a finite number of equally-spaced nodes, as shown in Fig. 10. Each peridynamic node is representative of a finite volume 훥푉 = ℎ훥푥 2 , where 훥푥 = 훥푦 is the grid spacing and ℎ is the thickness of the body. Please note that the most external real nodes do not lie exactly on the boundary of the body since the nodes are positioned at the center of the volume 훥푉. The ratio between the horizon and the grid spacing is defined as 푚-ratio: 푚 = 훿/훥푥. The value of this parameter determines the density of peridynamic nodes within a neighborhood. Furthermore, the fictitious layer (empty dots in Fig. 10) is added to the real body to complete the neighborhoods of the nodes near the boundary.

Numerical Taylor-based extrapolation method
The Taylor-based extrapolation procedure in the discretized model aims to determine the values of the variables of the fictitious nodes (displacements 푢 and 푣 respectively in 푥 and 푦 directions and dilatations 휃 in the case of a peridynamic 2dimensional body). The numerical procedure to determine, for instance, the displacement 푢 푖 in a fictitious node 푖 with coordinates (푥 푖 , 푦 푖 ) is carried out as follows: -find the real node of index 푗 closest to node 푖; -perform a Taylor series expansion of 푢 푖 about node 푗 with coordinates (푥 푗 , 푦 푗 ): where 푢 푗 and 휕 푛 1 +푛 2 푢 푗 휕푥 푛 1 휕푦 푛 2 are the displacement of node 푗 and its derivatives, 푛 푚푎 푥 is the maximum order of the truncated Taylor series, 푛 1 and 푛 2 are the orders respectively in 푥 and 푦 directions and 푛 is the global order so that 푛 2 = 푛 − 푛 1 .
Since the coordinates of nodes 푖 and 푗 are known, the displacement 푢 푖 in Eq. 37 is written as a function of the displacement 푢 푗 and its derivatives. However, we aim to express 푢 푖 as a function solely of the displacement of the real nodes. The total number of derivatives of 푢 푗 , considered before truncating the Taylor series, is 푛 푑 = (푛 푚푎 푥 (푛 푚푎 푥 + 1)/2) − 1.
They can be determined as functions of the displacements of the 푛 푑 real nodes near node 푗 by following another Taylorbased extrapolation procedure: -find the 푛 푑 real nodes with indices 푗 푘 closest to node 푗, where 푘 = 1, . . . , 푛 푑 (see Remark below for the conditions on the node search); -for each of those nodes with coordinates (푥 푗 푘 , 푦 푗 푘 ), perform a Taylor series expansion of their displacements 푢 푗 푘 about node 푗: -solve the system of equations in Eq. 38 to obtain the derivatives of 푢 푗 as a function of the displacements 푢 푗 and 푢 푗 푘 : Therefore, by combining Eqs. 37 and 39 , the displacement of a fictitious node is a function only of the displacements of some real nodes. Note that the adopted nearest-node startegy is really simple to implement also for complex geometries. This procedure can be applied to determine the displacements 푢 and 푣 and the dilatations 휃 of all the fictitious nodes.
i j 5 j 4 j 2 j 3 j 1 j Fig. 11 Taylor-based extrapolation method for a fictitious node 푖 near a corner of the body: node 푗 is the real node closest to node 푖 and nodes 푗 푘 with 푘 = 1, . . . , 5 are the real nodes closest to node 푗 In the case of the dilatations, the truncation order 푛 푚푎 푥 must be replaced by 푛 푚푎 푥 − 1.

Remark 3
There might be some cases in which the system of equations (Eq. 38) is not solvable for the nodes 푗 푘 that are the closest to node 푗. For instance, if we want to determine the second derivative in 푥 direction, the nodes 푗 푘 must include at least two 푥 푗 푘 coordinates different from each other and from 푥 푗 (see example in Sect. 4.5). However, given the adoption a uniform grid in which nodes on the same lines share the same coordinates, this condition is not always met when searching for nodes 푗 푘 via the closest-node strategy without any additional condition. Therefore, in order for the system of equations to be solvable, the nodes 푗 푘 should comprise at least 푛 1 different 푥 푗 푘 coordinates and 푛 2 different 푦 푗 푘 coordinates (excluding 푥 푗 and 푦 푗 ) for each derivative of order 푛 1 in 푥 direction and 푛 2 in 푦 direction.
In the following we present an example for determining the displacement 푢 푖 of a fictitious node 푖 near a corner of the body by means of the Taylor-based extrapolation method with 푛 푚푎 푥 = 2. As shown in Fig. 11, the node 푗 near the corner is the real node closest to node 푖. Thus, the displacement 푢 푖 can be given via a Taylor series expansion about node 푗 (see Eq. 37) as where 푙 푥 = 푥 푖 − 푥 푗 and 푙 푦 = 푦 푖 − 푦 푗 . Note that the number of derivatives of the displacement 푢 푗 is 푛 푑 = 5. As shown in Fig. 11, 푗 푘 with 푘 = 1, . . . , 5 are the 5 indices of the real nodes closest to node 푗. Note that, in order to be compliant with the condition given in Remark 3, the search for the closest nodes should be carried out in terms of Manhattan distance. A system of 5 equations is written by performing a Taylor series expansion of 푢 푗 푘 about node 푗 (see Eq. 39): The factors of the Taylor series expansions, which are multiplied by the derivatives of 푢 푗 , are easily derived from Fig. 11. After some manipulations, the system in Eq. 41 yields: Therefore, by substituting Eq. 42 in Eq. 40, the displacement 푢 푖 of the fictitious node is expressed as a function solely of the displacements of the real nodes. This procedure can be repeated for the required variables of all the fictitious nodes.

Numerical formulation of peridynamics
Consider a real node 푖, as shown in Fig. 12. The neighborhood H 푖 of node 푖 embeds the complete volume of the nearest nodes and the partial volume of the nodes near the horizon limit. Therefore, all the nodes with at least a portion of their own volume within the horizon limit are considered part of the neighborhood H 푖 . For each family node 푗, the volume correction coefficient 훽 푖 푗 ≤ 1 is computed as the fraction of volume actually contained in the neighborhood [36]. If 훥푉 of node 푗 is completely inside the neighborhood, then 훽 푖 푗 = 1. The bond 푖 푗, which connects node 푖 to node 푗, could be either a real bond or a fictitious bond. In both cases, 12 The neighborhood H 푖 of a node 푖 is constituted by the nodes (black dots) with at least a part of their volume inside the horizon (gray line). The volume correction coefficient 훽 is the fraction of the volume of the family nodes 푗 within the horizon limit its reference scalar state 푥 푖 푗 , Gaussian influence function 휔 푖 푗 and inclination 휙 푖 푗 with respect to the 푥-axis, can be computed from the coordinates of the two nodes. Under the assumption of small displacements, the extension scalar state of bond 푖 푗 is given as where 푢 푖 and 푢 푗 are the displacements in 푥 direction respectively of nodes 푖 and 푗 and 푣 푖 and 푣 푗 are the displacements in 푦 direction respectively of nodes 푖 and 푗. If the family node 푗 is fictitious, Eq. 43 holds and 푢 푗 and 푣 푗 are determined as a function of the displacements of the real nodes by the Taylor-based extrapolation method exposed in Sect. 4.1.
The weighted volume of node 푖 is evaluated by performing a mid-point Gauss quadrature from Eq. 13: Since the neighborhoods of all the real nodes are complete thanks to the presence of the fictitious nodes, the weighted volume is constant in the whole body. Furthermore, the value of the weighted volume of the real nodes is assigned also to the fictitious nodes, as dictated by the Taylor-based extrapolation method.
Similarly, the dilatation of a real node 푖 is computed from Eq. 14 as On the other hand, the dilatations of the fictitious nodes are determined as a function of the dilatations of the real nodes by means of another Taylor-based extrapolation, as illustrated in Sect. 4.1.
The magnitude of the pairwise force density in bond 푖 푗 is given from Eq. 20 as Note that the constants 푘 휃 and 푘 푒 are determined by the constitutive modelling of the peridynamic theory, the parameters 푚 푖 , 푚 푗 , 휔 푖 푗 , 푥 푖 푗 and 훽 푖 푗 depend only on the geometric coordinates of the nodes in the reference configuration and the variables 푒 푖 푗 , 휃 푖 and 휃 푗 can be written as functions of the displacements of the real nodes (by using the proposed Taylor-based extrapolation method for the variables of the fictitious nodes). Therefore, by combining Eqs. 43-46 together, one can write an equation for each bond 푖 푗, either real or fictitious, to relate the magnitude f 푖 푗 of its pairwise force density to the displacements of the real nodes.
Finally, under the assumption of small deformation, the peridynamic equilibrium equation (multiplied by the node volume 훥푉) is written for every real node 푖 as where m 푖 푗 = {cos 휙 푖 푗 , sin 휙 푖 푗 } ⊤ is the bond direction in the reference configuration and b 푖 is the external force density vector applied to node 푖. The system of equations in Eq. 47 can be rewritten in the standard form where K is the peridynamic stiffness matrix (size: 2푁 × 2푁), u is the displacement vector (size: 2푁 × 1) and f is the force vector (size: 2푁 × 1). 푁 is the number of real nodes. The stiffness matrix K includes the contributions of the fictitious bonds, thus it embeds the correction of the surface effect by means of the Taylor-based extrapolation method.

Numerical evaluation of the peridynamic stress tensor
This section deals with the numerical procedure to compute the peridynamic stress tensor. The theoretical background can be found in Sect. 2.2.
In the discretized peridynamic model, we name the nodes corresponding to the points x, x ′ and x ′′ (see Fig. 2) respectively as 푖, 푗 and 푘. Under the assumption of point 푖 being in the bulk of a body subjected to a homogeneous deformation (see Eq. 23), the peridynamic stress tensor can be computed as where 푟 푖푘 is the length of the bond 푖푘. Equation 49 provides a good approximation also in the case of non-homogeneous deformations if the horizon 훿 is sufficiently small, as shown in [10] for bond-based peridynamics. However, Eq. 49 is not valid for nodes near the boundary which are affected by the surface effect (if no fictitious layer is employed).
In order to compute numerically the peridynamic stress tensor 푖 in a general node 푖, also not in the bulk of the body, the integrand of Eq. 22 should be evaluated for each bond 푗 푘 between node 푗 and node 푘. The differential volume d푉 x ′′ corresponds simply to the finite volume of node 푘, i.e., 훥푉. On the other hand, we must distinguish two types of bonds, which are shown in Fig. 13, to determine 훥푠 as the corresponding of the differential length d푠 of point x ′ in the direction of the bond m: where the trigonometric functions are within the absolute value because 훥푠 > 0. In the former case bond 푗 푘 is a type-A bond and, in the latter, a type-B bond, which are respectively shown in Fig. 13a and b. A type-A bond contributes to 푖 if it intersects the area 훥퐴 A 푖 = ℎ훥푦 passing through node 푖 perpendicular to 푥 direction, where ℎ is the thickness of the 2-dimensional body. Similarly, a type-B bond contributes to 푖 if it intersects the area 훥퐴 B 푖 = ℎ훥푥 passing through node 푖 perpendicular to 푦 direction. Therefore, the peridynamic stress tensor in a node 푖 is given as where f 푗 푘 is the magnitude of the pairwise force density of bond 푗 푘 obtained with Eq. 46, m 푗 푘 = {cos 휙 푗 푘 , sin 휙 푗 푘 } ⊤ is the bond direction and 훼 푗 푘 is a correction coefficient, given as: where 휕 퐴 푖 is the boundary of the area 훥퐴 푖 , which can be referred either to type-A bonds (훥퐴 A 푖 ) or to type-B bonds (훥퐴 B 푖 ). The different cases in Eq. 52 are illustrated in Fig. 14: -in the case with 푗 = 푖 (see Fig. 14a), since only half of the length 훥푠 related to node 푗 is on the opposite side of 훥퐴 푖 with respect to node 푘, then 훼 푗 푘 = 1/2; -similarly, in the case with 푘 = 푖 (see Fig. 14b), since only half of the volume 훥푉 of node 푘 is on the opposite side of 훥퐴 푖 with respect to node 푗, then 훼 푗 푘 = 1/2; -in the case that the bond intersects the boundary of 훥퐴 푖 , i.e., 푗 푘 ∩휕 퐴 푖 ≠ ∅, since 휕 퐴 푖 overlaps the boundary of the area of another node (node 푞 in Fig. 14c), the magnitude of the pairwise force density of bond 푗 푘 is equally shared between those nodes and, therefore, 훼 푗 푘 = 1/2; -in the case that the bond intersects 훥퐴 푖 not on its boundary, i.e., 푗 푘 ∩ ( 훥퐴 푖 \ 휕 퐴 푖 ) ≠ ∅ (see Fig. 14d), the magnitude of the pairwise force density of bond 푗 푘 contributes entirely to 푖 and, therefore, 훼 푗 푘 = 1; Moreover, to improve the computational efficiency, one might remove the factors 1/2 in Eq. 51 and consider each bond just once (for example consider only bond 푗 푘 and not bond 푘 푗). The proposed numerical procedure to compute the peridynamic stress tensor is used to highlight the surface effect in a 2-dimensional body in Sect. 2.4. Figures 7 and 8 show that the numerical computation of the peridynamic stress tensor is very close to the analytical solutions, obtained in Sect. 2.2, for nodes with a complete neighborhood.

Numerical evaluation of the force flux
Consider a finite area 훥퐴 which constitutes one of the sides of the volume cell of a node, as shown in Fig. 15. The numerical procedure to compute the peridynamic force flux through the finite area 훥퐴 is exposed in this section.
The force flux (x, n) of a point x in a direction n is interpreted in Sect. 2.3 as the sum of the pairwise forces per unit area of all the bonds intersecting the differential area d퐴 x on the plane P in x, where P is the plane passing through x perpendicular to n (see Figs. 4 and 5 ). Therefore, the force (a) (b) Fig. 13 Examples of type-A and type-B bonds for the computation of the peridynamic stress tensor 푖 in node 푖. The bonds contribute to 푖 only if they intersect the corresponding area 훥퐴 푖 flux through 훥퐴 can be discretized from Eqs. 30 and 32 as where x and n are respectively the centroid and the normal of 훥퐴, f 푗 푘 m 푗 푘 훥푉 2 is the pairwise force of any bond 푗 푘 intersecting 훥퐴 and 훼 푗 푘 is the correction coefficient given by Eq. 52 (see also in Fig. 15a, b the possible cases in the com- putation of the force flux). Note that, in order to improve the computational efficiency, the factor 1/2 is removed because only bonds satisfying the condition m 푗 푘 · n > 0 are considered. Since we are mostly interested in the force flux computed at the boundary of the body to impose Neumann boundary conditions, this concept is further analyzed in the next section.

Numerical implementation of the peridynamic boundary conditions
The real nodes closest to the boundary are not exactly on the boundary of the body (see Fig. 10). Therefore, by following the concepts exposed in Sect. 3.2, the boundary conditions in the discretized model should be imposed on the sides of the volume cells which overlap the boundary. We introduce a new cathegory of nodes, called "boundary nodes", at which the boundary conditions are imposed. Each boundary node is positioned at the centroid of the side of the volume cell 16 Boundary nodes at the boundary of the body: each node 푏 is representative of a finite area 훥퐴 푏 and is associated to the normal n 푏 external to the body of the nodes closest to the boundary and is representative of the finite area 훥퐴 푏 of that side, as shown in Fig. 16. As the fictitious nodes, the boundary nodes do not constitute new degrees of freedom because their displacements are determined as a function of the displacements of the real nodes by means of the Taylor-based extrapolation method. Suppose that the problem requires a constraint 푢 for the displacement 푢 푏 in 푥 direction of a boundary node 푏, condition given as 푢 푏 = 푢. The Taylor-based extrapolation method is applied to the boundary node exactly as done for the fictitious nodes in Sect. 4.1. The following procedure is valid also for 푢 = 0. The displacement 푢 푏 of the boundary node 푏 with coordinates (푥 푏 , 푦 푏 ) is determined by a Taylor series expansion about node 푗 with coordinates (푥 푗 , 푦 푗 ) as where node 푗 is the real node closest to node 푏. The 푛 푑 derivatives of 푢 푗 can be expressed as functions of the displacements of the 푛 푑 real nodes close to node 푗. Therefore, the Dirichlet boundary condition can be written as a function of the displacements of some real nodes: For example, we consider the case shown in Fig. 17 with a truncation order 푛 푚푎 푥 = 2. The Taylor series expansion of the displacement 푢 푏 about node 푗 is given as In order to determine the two derivatives in 푥 direction of Eq. 56, we find nodes 푗 1 and 푗 2 , shown in Fig. 17, as the nodes closest to node 푗 having 푥 coordinates different from each other and from 푥 푗 (see Remark 3). We write the system b u b = u j j 1 j 2 Fig. 17 Example of the Taylor-based extrapolation method used on a boundary node 푏 of equations given by the Taylor series expansions of 푢 푗 1 and 푢 푗 2 about node 푗, as in Eq. 38, and solve it to obtain the needed derivatives: By substituting Eq. 57 in Eq. 56, the constraint condition 푢 푏 = 푢 is given as a function of the displacements of some real nodes: More in general, the proposed method allows to write the Dirichlet boundary conditions as functions of the displacement vector u: Suppose now that an external force per unit area p is applied to a boundary node 푏, as shown in Fig. 18. The Neumann boundary condition is written in terms of force flux through the area 훥퐴 푏 associated to node 푏: where x 푏 is the position of node 푏 and n 푏 is the unit vector perpendicular to 훥퐴 푏 external to the body. Since the pairwise force density of any bond can be expressed as a function of the displacements of the real nodes (see Sect. 4.2), also (x 푏 , n) is a function of those displacements. Thus, Neumann boundary conditions can be written as functions of the displacement vector u: where B is the matrix of the boundary conditions (size: 2푁 푏 × 2푁), u is the displacement vector (size: 2푁 × 1) and c is the vector of the known terms (size: 2푁 푏 × 1). 푁 푏 is the number of boundary nodes. In order to include the boundary conditions (B u = c) in the system of equations derived from the equilibrium of the real nodes (K u = f), we conveniently use the technique of the Lagrange multipliers [5]. The vector of the Lagrange multipliers (size: 2푁 푏 × 1) is introduced in the system as The displacement vector u, extracted from the vector { u, } ⊤ , is the solution to the system of equilibrium equations which satisfies the imposed boundary conditions.

Numerical examples
Several examples are presented to verify the reliability and accuracy of the proposed method. Whenever possible, the numerical peridynamic results are compared with the reference solutions derived from classical continuum mechanics. The reference solution coincides with the peridynamic solution only in the limit of the horizon 훿 approaching 0 [26,43]. Therefore, the "difference" between these solutions includes two components: a discrepancy due to the different (local and non-local) formulations of the theories and the actual error given by the discretization and the implementation of the peridynamic model (either with or without the proposed method). The difference (in percentage) of the displacements between the peridynamic numerical results and the reference Poisson's ratio solution is computed at each node 푖 as where 푢 푖 and 푣 푖 are the displacements of node 푖, 푢 푟 푒 푓 (x 푖 ) and 푣 푟 푒 푓 (x 푖 ) are the displacements obtained with the reference solution at the position x 푖 of node 푖 and u 푟 푒 푓 is the displacement vector obtained with the reference solution at all the nodes. The reference solution is defined for each example by providing either the analytical solution (if possible), or the results obtained with the Finite Element Method. For simplicity sake, we consider a plate under plane stress conditions with different boundary conditions. The parameters adopted for the simulations of the plate are reported in Table 1. Firstly, we solve each example without adding the fictitious nodes. In this case, the boundary conditions are implemented by assigning the desired value of the constraints or loads to the most external nodes of the plate. In particular, Dirichlet boundary conditions are imposed by assigning to the nodes closest to the boundary the value of the displacement computed with the reference solution. Then, the same examples are solved by adopting the proposed Taylor-based extrapolation method and by implementing the boundary conditions as described in Sect. 4.5.

Plate under traction
The boundary conditions of the first example are shown in Fig. 19. The analytical solution is given by classical continuum mechanics as where 푝 is the traction load.
x y p Fig. 19 Boundary conditions for the plate under the traction 푝 = 1MPa The plots in Fig. 20 show the difference of the displacement field, computed without adopting any correction to the peridynamic model, with respect to the reference solution. The surface effect and the approximated way of imposing the boundary conditions lead to large errors near the boundary of the plate, especially near the corners. On the other hand, there are no fluctuations in the displacement field when the proposed Taylor-based extrapolation method with 푛 푚푎 푥 = 1 is employed. The differences of the displacement field of the corrected model (see Fig. 21) decrease sensibly with respect to those obtained without corrections at the boundary. Similar results are obtained by choosing higher orders for the Taylor-based extrapolation.
In the case adopting the proposed method, the error can be further reduced by increasing the accuracy of the integration over the neighborhoods, i.e., by increasing 푚 [36]. In order to show this, we compute the relative difference (in percentage) at a node 푖 as In the case employing the Taylor-based extrapolation on the fictitious nodes, the relative difference of every node inside the body is the same. Therefore, we gather in Table 2 the relative errors for different values of 푚. We can observe that the relative differences decrease significantly as the value of 푚 increase.

Plate under shear load
We present another example considering a plate under a shear load 푝. Figure 22 shows the boundary conditions of this case. Classical continuum mechanics yield the following analytical solution: (67) The differences of the displacement field, computed without corrections at the boundary, with respect to the reference solution are shown in Fig. 23. We observe that the differences are lower with respect to the case of the plate under traction in Sect. 5.1. However, as highlighted in Eq. 67, one would expect the displacements 푣 in 푦 direction to be 0 in the whole body. This fact is not verified in the numerical simulation because of the surface effect (see plot of the component 휏 22 of the peridynamic stress tensor in Fig. 8). This problem is completely solved by implementing the proposed Taylor-based extrapolation method with 푛 푚푎 푥 = 1, as shown in Fig. 24. Also, the differences of the displacements 푢 in 푥 direction decreases with respect to those obtained without corrections at the boundary. Similar results are obtained by choosing higher orders for the Taylor-based extrapolation.
As in the case of the plate under traction, we compute the relative differences 휖 푟 푒푙 푢 with Eq. 66 for different values of 푚 when implementing the proposed method in the case of the plate under shear load. For each 푚, the relative differences of the displacements 푢 are constant in the whole body also in this case and they are reported in Table 3. The numerical results show a significant reduction of the differences when the numerical integration is improved.

Plate under sinusoidal load
In the previous examples, the linear variation of the displacement field was properly captured by the Taylor-based extrapolation with the order 푛 푚푎 푥 = 1 (or higher). We investigate now an example in which the order of the Taylor-based extrapolation does not match the order of variation of the displacement field. The boundary conditions of this example are shown in Fig. 25. The force density applied thoughout the plate is given as where 푏 = 10 6 N/m 3   , where the origin of the reference system is at the center of the plate. Given the symmetry of the boundary conditions, the displacements 푣 in 푦 direction on the 푥-axis are fixed to be 0 coordinates with some of the FEM nodes, so that the peridynamic results can be compared with the reference solution. We expect that, by increasing 푛 푚푎 푥 , the Taylor-based extrapolation would approximate better the displacements of the fictitious layer. Figures 27 and 28 show the differences between the numerical results and the reference solution for the numerical models either without employing any correction for the

Crack propagation near the boundaries
A qualitative study is hereinafter conducted to investigate the behavior of crack growth near the boundaries of the body. We compare the results provided by the proposed method with the solution obtained with peridynamics when no corrections for the surface effect are adopted and the boundary conditions are imposed in a local way, i.e., constraints and load are applied only at the nodes closest to the boundary. The surface effect near the new boundaries generated by the crack growth is not corrected in the present paper, but it will be dealt with in future works. In order to model fracture phenomena, we introduce the scalar 휇 which yields the status of the bond (unbroken or broken) [41]: where 푠 is the stretch of the bond and 푠 푐 is the critical stretch for plane stress conditions. These quantities are respectively computed as [25] where 퐺 0 is the energy release rate. The scalar 휇 is historydependent since a broken bond cannot be restored. The equilibrium equation is therefore modified as and the damage at each node can be evaluated as [41] For the quasi-static crack propagation, the sequentially linear analysis used in [24] is employed: not only the stiffness matrix should be modified accordingly by removing the contributions of the broken bonds, but also the same external tension should be applied solely through the residual unbroken fictitious bonds.

Crack propagation due to Dirichlet boundary conditions
The geometry and the boundary conditions of a plate with a pre-existing crack are shown in Fig. 29. The properties of the plate under plane stress conditions are: thickness ℎ = 0.005푚, Young's modulus 퐸 = 1GPa, Poisson's ratio 휈 = 0.2 and energy release rate 퐺 0 = 196J/m 2 . The constraint is given as 푢 = 0.001m. A grid spacing 훥푥 = 0.0025m is used and 푚 = 3 is chosen. The results obtained by means of the peridynamic model with no corrections at the boundary and with the proposed method are compared in Fig. 30, in which the difference of the displacements in a node 푖 is computed as where u 푐표푟푟 and v 푐표푟푟 are the displacement fields respectively in 푥 and 푦 directions. The superscript 푢푛푐표푟푟 stands for "uncorrected" (peridynamic model with no corrections at the boundaries) and 푐표푟푟 for "corrected" (peridynamic model with the Taylor-based extrapolation method with 푛 푚푎 푥 = 1). It can be noticed that there are non-negligible differences near the crack tip that may lead to different crack paths. As shown by the damage of the nodes in the final configuration of the model with no corrections (see Fig. 31), the crack branches and reaches the upper edge of the plate in two separated paths. On the other hand, when the Taylorbased extrapolation method is adopted, the crack propagates to the upper edge in a unique path and there is no branching phenomenon. Hence, the crack path may change near the boundaries if the surface effect is mitigated and the boundary conditions are imposed in a "peridynamic way".

Crack propagation due to Neumann boundary conditions
The geometry and the boundary conditions of a plate with a pre-existing crack are shown in Fig. 33.  36 Crack propagation in the pre-cracked plate for the case with Neumann boundary conditions when the Taylor-based extrapolation method is used As in Sect. 6.1, the differences between the displacements obtained with and without the proposed method are computed as in Eq. 74. The differences near the tip of the preexisting crack (see Fig. 34) are non-negligible and they may lead to different crack behaviors. The crack paths, shown in Figs. 35 and 36 respectively for the case without and with boundary corrections, are indeed different between each other. Therefore, the crack behavior may be modified by the mitigation of the surface effect and the proper imposition of the Neumann boundary conditions.

Conclusions
Two issues arising near the boundary of a body modelled with ordinary state-based peridynamics are addressed: -the surface effect, i.e., the stiffness fluctuation near the boundary; -the current lack of standard strategies to impose the boundary conditions.
The surface effect has been studied numerically by evaluating the peridynamic stress tensor with a novel discretization method (see Sect. 4.3) and the characteristic hardening/softening behavior towards the boundary has been highlighted (see Figs. 7 and 8 ). This issue has been addressed by introducing a fictitious layer that completes the partial neighborhoods of the nodes near the boundary. We proposed a new version of the Taylor-based extrapolation method adopting the nearest-node strategy: the displacements of the fictitious nodes are expressed as functions of the displacements of the closest real nodes by means of multiple Taylor series expansions. In this way, the surface effect is mitigated.
The fictitious layer is also exploited to impose the boundary conditions in a peridynamic way. The boundary of the body is discretized by the so-called "boundary nodes". As the fictitious nodes, the boundary nodes do not constitute new degrees of freedom because their displacements are obtained with the Taylor-based extrapolation method. On the one hand, Dirichlet boundary conditions are implemented by constraining the boundary node and, accordingly, the fictitious layer mitigates the surface effect. On the other hand, Neumann boundary conditions are applied, via the numerical computation of the force flux at the boundary, to the bonds involving the fictitious nodes. Therefore, the boundary conditions are imposed in a "peridynamic way".
Several numerical examples were presented to verify the accuracy of the proposed method. The numerical results obtained with the Taylor-based extrapolation method show a great improvement with respect to the peridynamic models without corrections at the boundary. Furthermore, the order of the Taylor-based extrapolation can be increased until the undesired fluctuations of the numerical results become negligible for the application of interest. It is also shown that the numerical integration of the peridynamic equilibrium equation plays a fundamental role, so that the numerical results are improved even further by increasing the value of the 푚-ratio.
Moreover, we carried out a qualitative study on crack propagation near the boundaries by comparing the results obtained by means of the proposed method with those of the peridynamic model without boundary corrections. We presented two numerical examples in which the crack paths are different because of the difference in the displacement fields. This highlights the importance of mitigating the surface effect and of imposing properly the peridynamic boundary conditions.

Declarations
Conflicts of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.