A novel and effective way to impose boundary conditions and to mitigate the surface effect in state‐based Peridynamics

Peridynamics is a nonlocal continuum theory capable of modeling effectively crack initiation and propagation in solid bodies. However, the nonlocal nature of this theory is the cause of two main problems near the boundary of the body: an undesired stiffness fluctuation, the so‐called surface effect, and the difficulty of defining a rational method to properly impose the boundary conditions. The surface effect is analyzed analytically and numerically in the present paper in a state‐based peridynamic model. The authors propose a modified fictitious node method based on an extrapolation with a truncated Taylor series expansion. Furthermore, a rational procedure to impose the boundary conditions is defined with the aid of the fictitious nodes. In particular, Neumann boundary conditions are implemented via the peridynamic concept of force flux. The accuracy of the proposed method is assessed by means of several numerical examples for a state‐based peridynamic model: with respect to the peridynamic model adopting no corrections, the results are significantly improved even if low values of the truncation order for the Taylor expansion are chosen.

• In the most external layer of the body the peridynamic points have an incomplete neighborhood and "feel" the lack of the missing points beyond the boundary. This phenomenon, known as surface effect or skin effect, causes an undesired variation in stiffness near the boundary. In bond-based Peridynamics, the surface effect appears as a softening in the body external layer of thickness , 3,4 whereas in state-based Peridynamics the discrepancy due to the surface effect influences the points closer than 2 to the boundary, with a behavior which is difficult to predict a priori. 5 This behavior in a one-dimensional (1D) state-based peridynamic model is analyzed analytically and numerically in Section 2.2.
• The other issue related to nonlocality near the external surface of the body occurs when imposing the boundary conditions. In bond-based Peridynamics, Dirichlet boundary conditions have been commonly imposed on the peridynamic nodes closest to the boundary, named boundary nodes, whereas Neumann boundary conditions have been applied in several ways. For instance, in Reference 6 a heat flux is applied over the layer of thickness near the boundary assuming a linear variation of temperature there. However, this approach is theoretically inexact if the horizon does not tend to 0. 6 The "worst-case" scenario for the model accuracy near the boundary is given by the imposition of the external load as a force applied solely to the boundary nodes. 5 In state-based Peridynamics, a similar analysis can be applied, but the effects of the approximated boundary conditions are exhibited in the body exterior layer of thickness 2 . 5 Many strategies have been proposed to mitigate the surface effect in bond-based Peridynamics by modifying the stiffness of the bonds near the boundary: the force normalization method, 3 the force density method, 7 the energy method, 8,9 and the volume method. 10 The approach of computing some bond correction factors to recover the proper strain energy density of the points near the boundary under constant deformation, has been employed for state-based peridynamic models in the energy method 8 and in the position-aware linear solid constitutive model. 11 A thorough comparison of these methods has been carried out in Reference 5, highlighting that the boundary inaccuracy is only partially mitigated because these model corrections do not deal with the approximated imposition of the boundary conditions.
The surface effect and the problems of the boundary conditions can be confined in a smaller region if the horizon approaches 0. Therefore, the variable horizon method [12][13][14] can be exploited to decrease the value of the horizon near the boundaries and restrict the external layer of the body affected by the surface effect. 4,[15][16][17] Coupling of Peridynamics with classical continuum mechanics (for instance,  provides an approach to completely eliminate the surface effect and the problems of the imposition of the boundary conditions in a discretized model: the external layer is modeled with a classical continuum mechanics method, for example with FEM, so that there are no peridynamic nodes with an incomplete neighborhood and the boundary conditions are implemented in a local way on FEM nodes. However, this strategy is not suitable for all the applications of interest, for instance, in the case of a crack initiation at the boundary of the body. Furthermore, a crack constitutes a boundary surface of the body which is inevitably modeled with a region of peridynamic nodes.
The method of the fictitious or ghost nodes was introduced in bond-based Peridynamics to overcome the issue of the surface effect by completing the neighborhoods of the nodes near the boundaries with a layer of peridynamic material which actually does not exist. 7,22 The strategy of the fictitious nodes has been later employed successfully also in state-based Peridynamics. 23 This approach has been exploited to impose the boundary conditions in bond-based peridynamic models: the fictitious nodes are forced to move according to some predefined functions in order to obtain the desired boundary displacement or load. In References 24,25 the displacement distribution of the fictitious nodes is chosen ad hoc for each case. In Reference 8 the constraints are imposed by displacing the fictitious nodes with the desired constant value and the external loads are applied in the boundary layer within the body. In References 5,16,26-31 the Dirichlet boundary conditions are implemented by imposing an odd function over the fictitious layer. Moreover, in References 5,16,27,30 the displacements of the fictitious nodes to enforce the desired load are derived from the formulae of classical mechanics. In Reference 17 the odd, polynomial and sinusoidal functions are studied to properly impose the desired constraint by means of an iterative procedure. A more theoretical approach for the fictitious layer method has been developed in  The Taylor-based fictitious node method, 31,35 an approach similar to the polynomial extrapolation function in Reference 17, has been proposed for bond-based Peridynamics. This approach finds, for every real node near the boundary, the corresponding fictitious node along the normal to the local tangent to the boundary and then evaluates its displacement with a linear extrapolation via the first-order truncated Taylor series. The desired boundary conditions, written as equations adapted from the form of a classical continuum mechanics model, are incorporated into the formulae of the Taylor series. 31 There are two main drawbacks using this approach: 31 the procedure of finding the corresponding fictitious node might be rather difficult in complex geometries and the displacement of each fictitious node is not unique because it depends on the real node at which the computation is carried out.
Inspired by the aforementioned works, we propose, for a state-based peridynamic model, a generalized method of the Taylor-based extrapolation over the fictitious nodes considering the nearest-node strategy. This strategy allows to avoid the two described drawbacks of the Taylor-based extrapolation method, since the displacements of the fictitious nodes are determined with a Taylor series expansion centered at the closest boundary node. Therefore, the uniqueness of the displacements of the fictitious nodes is guaranteed because they depend solely on the displacement of the closest boundary node and its spatial derivatives. Moreover, the nearest-node strategy is easy to implement and suitable for any kind of geometry, even very complex. In addition to that, we show how to implement the fictitious node method in state-based Peridynamics with a fictitious layer of thickness (instead of 2 as in Reference 23) via the truncated Taylor series up to a general order n. The surface effect is mitigated by the presence of the fictitious nodes and Dirichlet boundary conditions are applied similarly to what is done in References 31,35. On the other hand, we propose an innovative method to impose Neumann boundary conditions in a "peridynamic way", namely by means of the peridynamic concept of force flux.
The article is organized as follows: Section 2 presents a review of 1D state-based Peridynamics and the analytical analysis of the surface effect in a finite body through the concept of force flux, Section 3 exposes the proposed method and the imposition of the boundary conditions, Section 4 illustrates the discrete model and its implementation, Section 5 shows the comparison between the 1D numerical results obtained by using the proposed method and those produced by the state-based peridynamic model without corrections, Section 6 introduces the extension of the Taylor-based extrapolation method to two-dimensional (2D) problems, and finally Section 7 draws the conclusions of the paper and addresses new possible research developments.

1D STATE-BASED PERIDYNAMICS
The peridynamic theory is based on nonlocal interactions, the so-called bonds, acting between material points. Let Ω ∞ be an infinite 1D body modeled with Peridynamics, as shown in Figure 1. The horizon is defined as the minimum length at which the force interaction vanishes due to increasing distance. Therefore, a point at position x is affected by the points in its neighborhood H( A bond is the relative position vector in the reference configuration between the interacting points: In a deformed configuration, the relative displacement vector is defined as: where u is the displacement along the x-axis. Note that + is the relative position of the points in the deformed configuration.
The dynamics of point x in a homogeneous isotropic infinite body Ω ∞ is governed by the state-based peridynamic equation of motion: 2 where is the material density, T is the force density vector state, dV x ′ is the infinitesimal volume associated to the point at position x ′ and b is the external force density vector. The notation T[x, t]⟨ ⟩ means that the force density scalar state T F I G U R E 1 Reference configuration at time t = 0 and deformed configuration at time t depends on point x and time t and is applied to the bond vector . The integrand on the right-hand side of the equation represents the peridynamic internal force due to the bond acting on the point x (see Figure 2). This article is only concerned with the static solution of the peridynamic problems, [36][37][38] hence, dropping the dependence upon time, the peridynamic equilibrium equation becomes: where, for brevity, In the 1D model, the force density vector state can be defined as: where t⟨ ⟩ is the force density scalar state (magnitude of T) and M⟨ ⟩ is the deformed direction vector state (unit vector in the direction of T).
For later use, the reference position scalar state, representing the bond length in the reference configuration, is defined as: The notation ⟨ ⟩ will be dropped whenever unambiguous. The elongation of the deformed bond is described by the extension scalar state given as: Based on the aforementioned notions, the weighted volume of a peridynamic point x is given in the initial configuration as: where ⟨ ⟩ = is a prescribed spherical influence function.

Correspondence with classical mechanics theory
In order to ensure the correspondence between classical continuum theory and state-based Peridynamics, the elastic deformation energy density computed by both approaches should coincide under the same strain conditions. 2,39 According to classical continuum mechanics, the strain conditions for an isotropic three-dimensional (3D) bar under axial loading in which the cross-sections remain plane and perpendicular to the longitudinal axis are given as: where ij with i, j = 1, 2, 3 are the components of the strain tensor . The x-axis, denoted with 1, is considered the longitudinal axis of the bar. The classical dilatation cl is therefore computed for small strains as: The strain energy density for a linear elastic bar in classical continuum mechanics is computed from the general 3D formula 39 by substituting in it the strains of Equation (9): where k is the bulk modulus, the Poisson's ratio, G the shear modulus, E the Young's modulus, d ij with i, j = 1, 2, 3 are the components of the deviatoric strain tensor d = − ( cl ∕3) ⋅ 1 where 1 is the identity matrix with a size 3 × 3, and ij with i, j = 1, 2, 3 are the components of the stress tensor . Note that in the second line of Equation (11) the strain energy density is not divided into volumetric and deviatoric parts as in the first line. The second line of Equation (11) is useful to be compared with the peridynamic strain energy density to derive the peridynamic constants, although it could be rewritten in a more common form (third line) after some omitted manipulations.
The peridynamic strain energy density of a linear material at point x is computed as Reference 2: wherek andĜ are the peridynamic constants which must be determined by the comparison with the classical strain energy density. and e d are named, respectively, dilatation and deviatoric extension scalar state and their correspondence with the classical values are explained in the next paragraphs. The peridynamic dilatation of the point x is defined to match the classical dilatation under constant deformations: where the dilatation coefficient is c = 1 − 2 . This coefficient is computed through the correspondence of the peridynamic dilatation ∞ of an infinite 1D body Ω ∞ with the classical dilatation cl for an imposed small constant deformation 11 = (the extension scalar state is e = x for every bond): The deviatoric extension scalar state is defined as: The deviatoric extension scalar state e d is the elongation of the bond, which can be viewed as a material element in classical continuum mechanics, due to the deviatoric strain tensor (constant in the whole body Ω ∞ ): 39 where 1 , 2 , and 3 are the components of the bond vector in x, y, and z direction, respectively. In a 1D case the bonds lie on the longitudinal axis of the bar ( 2 = 3 = 0), thus e d = d 11 1 = d 11 x. The peridynamic deviatoric strain energy density W d of a point x in the infinite body Ω ∞ is then given as: The peridynamic elastic constants are derived by equating the classical and peridynamic strain energy densities: The force density scalar state t⟨ ⟩ is found by taking the Fréchet derivative of the peridynamic strain energy density 2 (see Appendix A for details): where k = ( E)∕((1 − 2 )(1 + )) and k e = E∕(1 + ). We showed here how to determine the force density scalar state t by means of the calibration against the classical strain energy density by considering the two simple loading conditions of isotropic expansion and simple shear. The reader can find in Reference 40 an alternative approach.  (5) and (19): where m ′ = m(x ′ ) and ′ = (x ′ ) are, respectively, the weighted volume and the dilatation of the peridynamic point x ′ . Equation (20) provides a conventional form of the equilibrium equation which equates the external force distribution, on the right-hand side, and the internal peridynamic forces, on the left-hand side. The term in the square brackets in Equation (20) represents the magnitude of the total force density in the bond , namely, t⟨ ⟩ + t ′ ⟨− ⟩.

Surface effect
The formulation of the peridynamic theory for a 1D finite body Ω of length l implies the existence of a boundary Ω, as shown in Figure 3. We individuate two sets of points near the boundary which will be useful to describe the surface effect: Ω = {x ∈ Ω ∶ ∀x ′ ∈ Ω, |x ′ − x| < } andΩ = {x ∈ (Ω ⧵ Ω) ∶ ∀x ′ ∈ Ω, |x ′ − x| < 2 }. These domains are shown in Figure 3, where Ω = Ω 0 ∪ Ω l andΩ =Ω 0 ∪Ω l . The derivation of the peridynamic constants k and k e is strictly valid only for points with a complete neighborhood. In a bond-based peridynamic model without any correction, the points in Ω with an incomplete neighborhood suffer a reduction in stiffness due to the lack of bonds with points beyond the boundary. 3 In state-based Peridynamics, the domain of the body which is affected by the surface effect is extended also toΩ because the peridynamic forces acting on point x depend not only on the properties of the bond ( ⟨ ⟩, x⟨ ⟩, e⟨ ⟩) and of the neighborhood of the point itself ( , m), but also on the properties of neighborhoods of the neighboring points ( ′ , m ′ ), as shown in Equation (20). Therefore, the state-based peridynamic solution of a finite body under a small constant deformation would differ from the infinite body solution in (Ω ∪Ω). 5 F I G U R E 3 Domains near the boundary Ω = Ω 0 ∪ Ω l of a finite 1D body Ω of length l: Ω = Ω 0 ∪ Ω l is the nonlocal boundary with a maximum distance from Ω,Ω =Ω 0 ∪Ω l is the set of points with a distance between and 2 from Ω and Γ = Γ 0 ∪ Γ l is the fictitious layer of thickness Let us define H + and H − of a point x as the halves of the neighborhood of that point, respectively, in the same and opposite direction of a unit vector n: The force flux at a point x ′′ in the direction of n (see Figure 4) is defined in 1D state-based Peridynamics as References 1,41: where H ′′ = H(x ′′ ) is the neighborhood of the point x ′′ . The physical meaning of this definition is that the force per unit area acting on the point x ′′ belonging to a plane with normal n is the resultant of the forces per unit area of the bonds which intersect that plane. 41 The force flux is the peridynamic concept closest to the stress in classical mechanics. 1 In order to quantify the discrepancy in stiffness due to the surface effect near the boundaries in a 1D state-based peridynamic model, the force flux given by a small constant deformation 11 = of the finite body Ω is evaluated and then compared with the force flux in the infinite body Ω ∞ under the same deformation . The magnitude of the force density in a bond of a body under a constant deformation (e⟨ ⟩ = x⟨ ⟩, (x) = c ), is given from Equation (20) as: where c k + k e = E. Since the influence function can be chosen arbitrarily, 42 a hyperbolic influence function is adopted: The choice of this influence function helps simplifying the analytical computation of the force flux, which is rearranged in the case of a uniformly deformed body via Equation (22) as: where dV ′ = Adx ′ with A indicating the cross-sectional area. The computation of the weighted volume m ∞ in an infinite body Ω ∞ is carried out according to Equation (8): The weighted volume is constant along the infinite body since all the points have a complete neighborhood. With this result, Equation (24) yields: Note that the force flux computed in each point of the infinite body Ω ∞ under the small uniform strain is equal to the stress 11 = E in classical mechanics. F I G U R E 4 Integration domains for the force flux definition: x varies within H ′′− and x ′ within H ′′+ ∩ H + Dropping the assumption of infinite body, we compute the force flux in a finite bar Ω of length l under a constant deformation . The weighted volume in (Ω ⧵ Ω) is computed exactly as in Equation (25) because in the bulk of the bar the neighborhoods of the peridynamic points are complete. On the other hand, the weighted volume of the points in Ω is smaller due to the lack of some points in their neighborhoods. The computation of the weighted volume in Ω 0 yields: The computation of m in Ω l is carried out analogously. The weighted volume distribution in a finite body normalized with respect to the value of the weighted volume in the infinite body, is reported in Table 1 and plotted in Figure 5 for different values of the horizon . Note that the normalized weighted volume differs from 1 in a more confined domain as the horizon decreases.
The computation of the force flux is carried out by substituting the functions of the weighted volume, given in Table 1, in Equation (24). The force flux in (Ω ⧵ (Ω ∪Ω)) is computed in an identical way in Equation (26) for the case of infinite body since every bond in that domain connects solely points with a full neighborhood. On the other hand, the force flux in the domains near the boundary is computed as follows: The force flux (x ′′ ∈ (Ω l ∪Ω l ), +1) is computed similarly and the result is symmetric with respect to the domain The force flux distribution normalized with respect to the infinite body solution, is reported in Table 2 and plotted for different values of the horizon in Figure 6. As expected, the material points in (Ω ⧵ (Ω ∪Ω)) behave as if the body were infinite, whereas the points in (Ω ∪Ω) "feel" the missing material points beyond the boundaries. Furthermore, it is observed that the variation of the horizon does not affect the values of the force flux, but only the size of the domain over which the surface effect is extended. Intuitively, as ∕l approaches to 0 (see the different curves in Figure 6), the domain which recovers the infinite body stiffness properties is enlarged.
In Figure 7, the normalized force flux near the boundary Ω is plotted with respect to the normalized coordinate x∕ . The characteristic behavior is driven by two counteracting phenomena: • due to the lack of material points and consequent missing bonds, the force flux is reduced because less bonds intersect the section; TA B L E 2 Normalized force flux distribution along the finite one-dimensional body • since the weighted volume decreases close to the boundary, the magnitude of the nearby forces in the bonds is increased (see Equation 19).
InΩ there are no bonds missing and the intensified magnitude of the peridynamic forces increases the stiffness of the material (hardening). On the other hand, both the effects are present in Ω: the intensification of the bond forces overcompensates the lack of bonds far from the boundary, but the effect of the missing bonds prevails near the boundary (softening).
Analyzing Figure 7, there is a maximum ≈ 1.263 ∞ at (x∕ ) ≈ 0.772. In addition, the force flux density matches again the classical stress at (x∕ ) ≈ 0.394. After some numerical investigations (see Figure 8), one can assert that the shape of the force flux distribution changes only slightly with a different, but reasonable choice of the influence function, but the overall behavior near the boundaries is the same: the material exhibits a hardening and then a softening from the bulk towards the boundary.
The imposition of the boundary conditions in Peridynamics is correlated to the problem of the surface effect because of the nonlocal nature of the theory. The nonlocal boundary conditions are often imposed in an approximated way at the boundary points, which causes undesired displacement fluctuations near Ω. In fact, according to Reference 25 the nonlocal constraints in state-based Peridynamics should be imposed over a region of thickness 2 outside the body, whereas the external loads should be applied in the domain (Ω ∪Ω) within the body, as shown in Appendix B. However, the distribution of the variables which should be imposed is not known a priori.
The method proposed in the following section compensates the surface effect and imposes in a proper way the nonlocal boundary conditions without any need to determine the appropriate distribution of the external displacement or force density.

TAYLOR-BASED EXTRAPOLATION OVER THE FICTITIOUS LAYER
As shown in Figure 3, a fictitious layer Γ = Γ 0 ∪ Γ l of thickness is added at the ends of the 1D body Ω. The Taylor-based extrapolation, used to determine the displacements of the fictitious points, is shown hereinafter. According to the nearest-point strategy, the Taylor series of the displacements of the fictitious points should be centered at the closest real point, namely, the boundary point. Therefore, given a maximum order n ≥ 1, the truncated Taylor series of the displacement in x ∈ Γ l about the boundary point x = l yields: With this Taylor-based extrapolation, the displacements of the fictitious points are determined as functions of the unknown displacement of the boundary point and its derivatives.
In state-based Peridynamics, the properties of the points in the fictitious layer should be corrected as well. Therefore, the weighted volume m of the fictitious points in Γ l is fixed to be m ∞ , namely, the weighted volume of a point with a complete neighborhood, and their dilatation is extrapolated by means of the following truncated Taylor series expansion: Note that the dilatation is a measure of the strain, thus the truncation of the relevant Taylor expansion occurs with 1 order less than that of the displacement. Equations (30) and (31) are valid in Γ l , but the same procedure is applied to the fictitious points in the domain Γ 0 . An example of the extrapolation with the Taylor series expansion for n = 1, 2, 3 is shown in Figure 9.

Mitigation of the surface effect
This section shows that, by implementing the method of the Taylor-based extrapolation over the fictitious points with a proper order n, the peridynamic solution of the force flux in a finite body Ω under a predefined smooth deformation recovers that of the infinite body.
The force flux at a point x ′′ in an infinite body Ω ∞ under a constant deformation is given from Equation (26) as ′′ ∞ = E . Consider a finite body Ω under the same constant deformation . Assuming, without loss of generality, u(0) = 0, the displacement imposed to the body is u(x ∈ Ω) = x. Thus, the proposed Taylor-based extrapolation method is applied to Ω with an order n = 1 (dotted line in the example of Figure 9). The displacement in the fictitious layer Γ l is determined from Equation (30) as: Since this Taylor-based extrapolation induces the fictitious layer to deform in the same way as the real body, the extension scalar state of the fictitious bonds is equal to that of the real bonds e = x and the dilatation of the real points near the boundary is computed exactly as Equation (14), yielding (x ∈ Ω l ) = c = ∞ . Therefore, thanks to the extrapolation of the displacement of the fictitious points (Equation 32), the force flux inΩ l is computed as the force flux of an infinite body in Equation (26): The dilatations of the fictitious points are determined by the Taylor-based extrapolation in Equation (31) as: Therefore, in this case the dilatation is constant in the fictitious layer and is equal to the dilatation ∞ of the infinite body under the deformation . Furthermore, the weighted volume of the fictitious points is fixed to be equal to the weighted volume of a point with a complete neighborhood: m(x ∈ Γ l ) = m ∞ = A 3 . Since the extension scalar state e of the fictitious bonds is computed with the same formula as that of the real bonds (from Equation 32) and the properties of the real and fictitious points are equal to those of the points in an infinite body (m(x ∈ (Ω l ∪ Γ l )) = m ∞ and (x ∈ (Ω l ∪ Γ l )) = ∞ ), the force flux near the boundary of the finite body is computed exactly as in Equation (26): The same results as those of Equations (33) and (35) can be obtained at the other end of the finite body through the extrapolation method applied to the fictitious layer Γ 0 .
We have shown that, when the proposed Taylor-based extrapolation method is employed with the order n = 1 (or higher), the force flux in the whole domain (even near the boundaries) of a finite body under a constant deformation is equal to the force flux in an infinite body under the same deformation. Appendix C presents the same procedure carried out for a cubic distribution of displacement imposed to a finite body, to which the proposed method is applied with n = 3. The force flux coincides with the one computed in the infinite body under the same condition in that case as well. Hence, we conjecture that, whenever the order n of the proposed Taylor-based extrapolation method is equal to (or higher than) the order of the imposed distribution of displacement, the points near the boundary of the finite body retrieve the stiffness properties of the points in the infinite body thanks to the fictitious layer, and the surface effect is eliminated.
Since the displacement distribution is rarely known before solving the problems of interest for real application, Section 5 presents some numerical examples in which the displacements are initially unknown and n is chosen differently from the ideal order.

Dirichlet boundary conditions
This section analyzes how to deal with the imposition of the nonlocal displacement boundary conditions by means of the proposed method. The displacement constraints are implemented as in a local model: the desired displacement is imposed at the boundary point and the proposed method autonomously adjusts the displacement of the fictitious points according to the displacement of the boundary point and its derivatives. The equation to impose a Dirichlet boundary condition, for instance, at Ω 0 is given as: where u is the desired value of the constraint. The extrapolation method determines the displacements and the dilatations of the fictitious points similarly to what was shown in Equations (30) and (31): Note that in this case u(0) in Equation (37) is constrained to be u. In this way, the Dirichlet boundary condition is imposed only at the boundary Ω, whereas the fictitious layer Γ mitigates the surface effect, as illustrated in Section 3.1.
For a detailed description of the imposition of the Dirichlet boundary condition in a discretized model, please refer to Section 4.3.

Neumann boundary conditions
This section addresses the problem of applying a load to the nonlocal boundary of the finite 1D body. With the introduction of the fictitious layer, new interactions between real and fictitious points, named fictitious bonds, are generated. These fictitious bonds deform because of the displacement induced by the extrapolation on the fictitious points. Thereby, each fictitious bond exerts a force on both the interacting points (see Figure 2). The forces which the fictitious bonds exert on the fictitious points are "lost", in the sense that they do not contribute to the equilibrium of the real points of the body. Indeed, the fictitious points do not constitute new degrees of freedom of the model since their displacement is extrapolated via Equation (30). On the other hand, the forces of the fictitious bonds acting on the real points contribute to their equilibrium equation and play a role very similar to that of the external force density distribution computed in Appendix B to obtain a constant deformation of the bar: the global effect of all the fictitious bonds is a force density distribution applied to the nonlocal boundary (Ω ∪Ω) capable of compensating the surface effect and imposing the Neumann boundary condition. Therefore, we propose to apply Neumann boundary conditions by means of the concept of force flux (Equation 21) at the boundary points: the resultant of the forces of the fictitious bonds which intersect the boundary is equal to the desired value of the external load (see Figure 13 for the discretized representation of the load boundary condition). In order to impose an outward external force f, for instance, at Ω l , the boundary condition is given as: where H l is the neighborhood of the point x = l. The integral in Equation (39) depends on the elongation (i.e., the extension scalar state) of the real bonds near the boundary and of the fictitious bonds crossing the boundary. Therefore, the Neumann boundary condition written as Equation (39) imposes the elongations of the bonds so that the sum of the forces per unit area in those bonds is equal to the desired value of the force flux at the boundary. The extension scalar state of the bonds involved in Equation (39) is determined by the displacement in (Ω l ∪Ω l ∪ Γ l ), namely, the domain near the boundary of the real body and the fictitious layer. However, the displacement distribution in the fictitious layer Γ l is defined by the Taylor-based extrapolation method. Therefore, the integral in Equation (39) is a function solely of the displacement distribution of the real points near the boundary. The implementation of a Neumann boundary condition in a discretized model, described in Section 4.3, is indeed carried out by adopting a linear combination of the displacements of the real nodes near the boundary. The discretized representation of Equation (39) is shown in Figure 13.
Remark 2. In Reference 30, the zero-traction boundary condition is applied by removing the fictitious layer. By doing so, however, the surface effect reappears near the free boundary. In the case of zero-traction surface, we suggest to use the fictitious layer with an extrapolation method, such as the one presented in this article, and impose the condition (x ∈ Ω, n) = 0. Some numerical examples of the implementation of a zero-traction boundary condition are shown in Section 5.
Remark 3. The Dirichlet boundary condition described in Section 3.2 is imposed at the boundary point Ω 0 of the body. Due to the nonlocal nature of the peridynamic theory, the reaction to the constraint is not the force density applied to that point, but is computed as the force flux at that point multiplied by the cross-sectional area: f u = A ⋅ (0, −1), where f u is the reaction to the constraint imposed with Equation (36).

NUMERICAL IMPLEMENTATION
A mesh-free method is adopted to discretize the domain. 43 For simplicity sake, the peridynamic grid consists of a finite number of equally spaced nodes, each of which is representative of a volume V = AΔ, where A is the cross-sectional area of the 1D body and Δ is the grid spacing (see Figure 10). Since the horizon is often chosen as a multiple of the grid spacing Δ, the m-ratio m = ∕Δ is here considered as an integer number. The parameter m is commonly denoted as m in Peridynamics, but the top bar avoids the confusion that might arise with the weighted volume symbol. The position of the nodes, shown in Figure

Numerical Taylor-based extrapolation
Applying the concepts for the continuum peridynamic formulation of Section 3, the Taylor-based extrapolation is used in the discretized formulation to express the displacements (or dilatations) of the fictitious nodes as functions of the displacements (or dilatations) of the real nodes close to them. Initially, the procedure for the linear extrapolation (n = 1) of the displacements in Γ is carried out and, thereafter, the procedure for a general extrapolation of order n is presented. This procedure is similar to that exposed for the polynomial function in Reference 17. The displacement of a fictitious node f in Γ is evaluated by means of the Taylor series expansion about its closest real node, according to the nearest-node strategy. If the node f , for instance, belongs to the domain Γ 0 , that is, f = −m, … , −1, the linear Taylor series about node 1 yields: The derivative in Equation (40) can be determined by another linear Taylor series about the boundary node 1 for the displacement of the real node closest to node 1: F I G U R E 10 Uniform peridynamic grid with node numeration: The solid dots represent the real nodes, whereas the empty dots represent the fictitious nodes Thus, the displacement u f can be expressed as a function of the displacements of the real nodes close to the boundary by substituting Equation (41) in Equation (40). Now, the procedure to define the displacements of the fictitious nodes in Γ 0 with an extrapolation based on the Taylor series expansion of order n is presented. The real node closest to the domain Γ 0 is node 1, thus the Taylor series expansion should be centered at that node according to the nearest-node strategy. The displacements of the fictitious nodes can be determined as functions of the displacement of node 1 and its n derivatives with the following system of equations: The vector of the displacements of the fictitious nodes in Γ 0 is named u Γ 0 (size: m × 1), the vector of the displacement of node 1 and its n derivatives is named d 1 (size: (n + 1) × 1) and the matrix obtained with the factors of the Taylor series expansion of the displacements of the fictitious nodes about node 1 is named T Γ 0 1 (size: m × (n + 1)). Therefore, Equation (42) can be compactly written as The derivatives in vector d 1 can be determined as functions of the displacements of the real nodes closest to the boundary node 1, according to the nearest-node strategy. Hence, one can write the Taylor series expansions of the displacements of the n + 1 real nodes closest to the boundary Ω 0 about the boundary node 1 as: which is compactly written as u Ω 0 = T Ω 0 1 d 1 , where the displacements of the real nodes closest to the boundary are gathered in vector u Ω 0 (size: (n + 1) × 1) and the factors derived by the Taylor series expansions about node 1 are contained in T Ω 0 1 (size: (n + 1) × (n + 1)). The vector d 1 is written as a function of the part u Ω 0 of the displacement vector u by inverting the Taylor matrix in Equation (43): Thus, the displacements of the fictitious nodes in Γ 0 can be expressed as a function of some of the unknowns of the problem: Note that the matrices in Equation (45) depend only on the coordinates of the fictitious and real nodes. The same procedure can be exploited for the extrapolation of the displacements in Γ l . Moreover, the dilatations in Γ can be analogously extrapolated with the Taylor series expansions of order n − 1.
The numerical procedure for the Taylor-based extrapolation over the fictitious layer Γ can be summarized as follows: • choose the order n of the truncated Taylor series; • find the boundary nodes i = 1, N, which are the real nodes closest to the fictitious layer (nearest-node strategy); • assemble the Taylor matrix T Γ i (size: m × (n + 1)) with the factors of the Taylor series expansion of the displacements u Γ of the fictitious nodes about node i (see Equation 42); • find the n + 1 real nodes closest to the boundary Ω (nearest-node strategy), whose displacements u Ω can be extracted from the displacement vector u; • assemble the Taylor matrix T Ω i (size: (n + 1) × (n + 1)) with the factors of the Taylor series expansion of the displacements u Ω about node i (see Equation 43); ] −1 , which expresses the displacements u Γ of the fictitious nodes as functions of the displacements u Ω of the real nodes close to the boundary: Analogously, the Taylor-based extrapolation of the dilatations Γ (size: m × 1) of the fictitious nodes can be carried out from the dilatations Ω (size: n × 1) of the n real nodes closest to the boundary as: where T Γ i (size: m × n) and T Ω i (size: n × n) are the submatrices obtained by eliminating the rows and the columns of T Γ i and T Ω i related to the node most distant from the boundary and the derivative with the highest order.

Discretized formulation
This section derives the discretized form of the peridynamic equations and shows the approach to combine them with the proposed method of the Taylor-based extrapolation over the fictitious layer. Consider a bond ij, either real or fictitious, connecting node i to node j, as shown in Figure 11. From Equation (1), the relative position vector can be computed for each bond ij of the peridynamic model as ij = x j − x i . The reference vector state x ij = | ij | and the hyperbolic influence function ij = ∕| ij | are computed for each bond, respectively, from Equations (6) and (23). Furthermore, the contribution to the neighborhood of nodes with their volume V only partially within it, is corrected by the volume reduction coefficient ij , computed as the fraction of volume actually involved in the neighborhood. 44 If V is completely inside the neighborhood, then ij = 1.
First, the displacements of the fictitious nodes are approximated with the extrapolation derived by the Taylor-based method previously described. Then, from Equations (2) and (7), the extension scalar state of bond ij is given as: Whenever one of the nodes of the bond is fictitious, its displacement is expressed as a function of some of the unknowns in the displacement vector u with the Taylor-based extrapolation methods exposed in Section 4.1. For instance, following the example of the linear Taylor-based extrapolation for a fictitious node f in Equations (40) and (41), the extension scalar state of a fictitious bond fj, connecting the fictitious node f in Γ 0 to the real node j, is expressed as: where the node indices possibly involved in this case are f = −m, … , −1 and j = 1, … , m (see Figure 10). As just shown with this example, the extension scalar state of all the bonds, both real and fictitious, can be computed as a linear combination of some of the displacements of the real nodes, that is, the unknowns of the problem. This statement is valid also F I G U R E 11 Bond ij between node i and node j in a one-dimensional peridynamic model for a higher order of the truncated Taylor expansion (with n > 1, please refer to the Taylor-based extrapolation procedure exposed in Section 4.1). The weighted volume of a node i is often evaluated in Peridynamics by performing a mid-point Gauss quadrature from Equation (8): where H i is the neighborhood of node i. However, given the choice of the hyperbolic influence function, the weighted volume is here computed from Equation (25) as m = A 3 . Note that, when the proposed method is employed, the weighted volume evaluated at the real nodes near the boundary and also at the fictitious nodes is equal to the one of a node with a complete neighborhood. The dilatation of a real node i is computed by means of the mid-point Gauss quadrature from Equation (13) as: Similarly to the procedure for the displacements, the dilatations of the fictitious nodes are extrapolated via Equation (47) as a linear combination of some of the dilatations of the real nodes. Thereby, the dilatations of all the nodes, both real and fictitious, are determined as functions of previously computed parameters.
The force density scalar state is then computed from Equation (19) as: In the case that the bond ij was fictitious, i and e ij have already been determined by means of the extrapolation method. From the integrand in Equation (20), the peridynamic force in bond ij is given as: where x ij = x ji , ij = ji , ij = ji , and e ij = e ji for symmetry reasons. f ij M ij is the bond force acting on the node i and f ji M ji = −f ij M ij is the same force of opposite direction acting on node j. Finally, the peridynamic equilibrium equation (Equation 4) is rewritten, in a discretized form, for every real node i as: where b i is the external force density applied on node i. Thus, Equation (54) is repeated N times, one for each real node, forming the system of equations to be solved. The right-hand side of the system equation is named body force vector f (size: N × 1). Since the internal bond forces f ij depend linearly on the displacements of the nodes (Equations 48-53), one can decompose the left-hand side of the system as Ku, where K (size: N × N) is the stiffness matrix and u (size: N × 1) is the displacement vector. Therefore, the system of equations can be written in the standard form: Appendix D shows the step-by-step assembly of K following the equations of this section. The stiffness matrix K embeds the stiffness correction obtained with the additional internal forces due to the fictitious bonds, thus the surface effect is compensated.

Application of boundary conditions in discretized form
The concepts about the peridynamic boundary conditions presented in Sections 3.2 and 3.3 are here translated in discretized terms: the boundary conditions matrix B is assembled to embed the system of equations derived from the boundary conditions. The Neumann boundary condition can be transformed in a relation among nodal displacements, achieving an equation similar to that of a constraint. Therefore, the boundary conditions are written in the form B u = c, where B (size: 2 × N) is the boundary condition matrix and c (size: 2 × 1) is the vector containing the known values of the boundary conditions. B and c have 2 rows in the 1D case because there is 1 condition for each boundary.
In order to show how to assemble B and c, suppose that the following boundary conditions are imposed to a finite body of length l: a constraint u at x = 0 and a force f at x = l. As shown in Figure 12, the desired constraint is not applied to any peridynamic node because each peridynamic node is centered with respect to its representative volume. However, the notions exposed in Section 4.1 can be exploited to extrapolate u(0) = u as a function of the nodal displacements. For instance, if a linear extrapolation was chosen (order n = 1), the displacement boundary condition would be written from Equations (40) to (41) as: Equation (56) can be rearranged as a row of matrix B and vector c as follows: In a similar way, the constraint for a Taylor series expansion of order n can be expressed as: where T Ω 0 1 (size: 1 × (n + 1)) is the matrix containing the factors of the Taylor series expansion of the displacement of the constraint at x = 0 about node 1, whereas T Ω 0 1 and u Ω 0 were derived in Equation (43). A similar equation can be written for the other boundary (x = l) if a displacement constraint is present there.
The extrapolation carried out in Equation (58) involves only the real nodes near the boundary and the coordinate of the boundary Ω 0 . However, the fictitious nodes play a role in the equilibrium equations of the real nodes near the boundary. The contribution of the fictitious layer Γ 0 , which mitigates the surface effect, is already taken into account in the stiffness matrix K (see Appendix D).
Consider now the imposition of an external load f at the boundary Ω l , as shown in Figure 13(A). The Neumann boundary condition is imposed through the concept of force flux: Therefore, the resultant of the forces of the fictitious bonds crossing the boundary is enforced to match the external force f (see Figure 13 Note that the matrix S l is defined in order to sum all the force of the fictitious bonds crossing the boundary Ω l . Equation (60) is a linear combination among the displacements of the real nodes near the boundary.
Hence, the boundary conditions for the current example are implemented as: Since Equation (61) expresses the boundary conditions as a relation among the displacements of the real points of the body, the technique of the Lagrange multipliers is very advantageous to include those boundary conditions in the system of equations given by the peridynamic equilibrium of the points (Equation 55). The Lagrangian function is introduced in the system as follows: 45 {

1D NUMERICAL EXAMPLES
Some benchmark problems are presented to verify the reliability and accuracy of the proposed method of the Taylor-based extrapolation over the fictitious nodes. The numerical peridynamic results are compared with the classical continuum mechanics solution, considered as the reference solution. The discrepancy between the numerical results and the classical solution is evaluated at node i by the relative percentage "error" : where u cl stands for displacement solution derived with classical continuum mechanics. As demonstrated in Reference 46, the classical solution coincides with the peridynamic one when the third-order and higher derivatives of the displacement are equal to 0 or when the horizon approaches 0. Since the numerical simulations must necessarily consider a finite value of , the classical solution can be considered the exact solution for the numerical peridynamic case solely if the highest displacement derivative is, at most, a second-order derivative. When this is not respected, a part of the difference between classical solution and numerical peridynamic solution is due to the different formulations of the two theories.
The data used in the benchmark problems is reported in Table 3. The value of the m-ratio m is chosen to be 3, as a convenient compromise between accuracy and computational cost. First, each peridynamic model is solved without the use of any correction method and, then, by means of the Taylor-based extrapolation method with different orders. In the peridynamic model adopting no corrections, a Dirichlet boundary condition is imposed by constraining the displacement of the closest node to its analytical value derived by classical continuum mechanics, whereas a Neumann boundary condition is imposed by applying the desired force to the node closest to the boundary. The imposition of the boundary conditions in the models which employ the Taylor-based extrapolation method, is described in Section 4.3.

Clamped bar under traction
The method is tested on a clamped bar with a traction load f = 10 4 N imposed at the end x = l of the bar, as shown in Figure 14. Solving the classical differential equation for an axially loaded bar with those boundary conditions yields:  Figure 15 shows the numerical peridynamic results of the clamped bar under traction loading, obtained without adopting any correction method and by implementing the proposed Taylor-based extrapolation method with the order n = 1. The hardening/softening behavior due to the surface effect near the boundaries is evident in the peridynamic model without corrections. On the other hand, when the method of the extrapolation over the fictitious layer is adopted, the results match perfectly the reference solution and the relative "error" , shown in Figure 16, is very close to the machine precision.

Clamped bar under uniformly distributed load
Consider now a clamped bar, with a free end at x = l, subjected to a uniformly distributed load b 0 = 10 6 Nm −3 , as shown in Figure 17. The solution derived by the classical continuum mechanics in this case is given as: Since d k u cl dx k = 0 with k ≥ 3, the classical solution coincides with the peridynamic one. 46

F I G U R E 15
Numerical results obtained with the peridynamic models without any correction (in blue) and with the Taylor-based extrapolation method with order n = 1 (in red) for the clamped bar under traction. The dashed black line represents the reference solution derived analytically by the classical continuum mechanics F I G U R E 16 Relative percentage "errors" for the peridynamic models without any correction (in blue) and with the Taylor-based extrapolation method with order n = 1 (in red) for the clamped bar under traction Figure 18 shows the numerical results for the peridynamic models with no correction and with Taylor-based extrapolation of orders n = 1, 2. The results of the former model exhibit a significant discrepancy with respect to the classical solution, whereas the results obtained by employing the proposed method are considerably closer to the reference solution.
Since the reference solution is a quadratic function, the order n = 1 for the Taylor-based extrapolation is not high enough to follow exactly that solution. However, the results are still greatly improved with respect to the model without fictitious layer (see blue and red lines in the relative "error" plot of Figure 19). This result is of paramount importance: in most cases the order of the truncated Taylor series to properly describe the variation of the displacement near the boundaries will be unknown, but a sufficiently high order will nonetheless provide a good approximation. On the other hand, choosing the order n = 2 (or higher) entails the numerical results to coincide with the analytical solution.

F I G U R E 17
Boundary conditions for the clamped bar under uniformly distributed load, which is acting in the direction of the bar F I G U R E 18 Numerical results obtained with the peridynamic models without any correction (in blue) and with the Taylor-based extrapolation method with orders n = 1, 2 (respectively, in red and orange) for the clamped bar under uniformly distributed load. The dashed black line represents the reference solution derived analytically by the classical continuum mechanics F I G U R E 19 Relative percentage "errors" for the peridynamic models without any correction (in blue) and with the Taylor-based extrapolation method with orders n = 1, 2 (respectively, in red and orange) for the clamped bar under uniformly distributed load

Clamped bar under linearly distributed load
The Taylor-based extrapolation method is tested in the case of a clamped-free bar under a linearly distributed load b 1 x, where b 1 = 10 6 Nm −4 , as shown in Figure 20. The solution of the nonhomogeneous differential equation of the axially loaded bar in classical continuum mechanics is given as: Since d 3 u cl dx 3 ≠ 0, the peridynamic solution is different from the classical one. 46 Therefore, a part of the difference, which is computed with Equation (64) as relative "error" between the numerical peridynamic results and the analytical classical solution, is due to the different formulations of the theories. Figure 21 shows the numerical results for the peridynamic models with and without employing the proposed method. As observed also in the plot of the relative "errors" of Figure 22, there is a great improvement in the accuracy of the results when adopting the Taylor-based extrapolation method of any order n. However, even if the same order of the expected solution is chosen (n = 3), the numerical results cannot match perfectly the classical solution because of the different formulation of the theories. Nonetheless, the choice of n = 3 is arguably the best one since the results do not exhibit any undesired fluctuation near the boundary due to the surface effect or the approximated imposition of the boundary conditions, whereas the results of the peridynamic models with n ≤ 2 do (see the right-hand side of Figure 22).

Clamped bar under sinusoidally distributed load
The proposed method has been proven to work really well with polynomially distributed loads. In order to further investigate the performance of the method, consider a clamped-free bar under a sinusoidally distributed load b s sin  Figure 23. The classical analytical solution in this case is given as: Since d k u cl dx k ≠ 0 with k ≥ 3, the peridynamic solution does not coincide with the classical one. 46 Therefore, a part of the difference is surely due to the different formulations of the theories. Figure 24 shows the numerical results obtained for the peridynamic model without adopting the proposed method and with the Taylor-based extrapolation of orders n = 1, 3, 5. Again, the proposed method greatly improves the accuracy of the solution (see the left-hand side of Figure 25). Furthermore, increasing the order n of the Taylor-based extrapolation leads to a gradual reduction of the result fluctuations due to the boundary issues, as one can observe in the right-hand side of Figure 25. Therefore, 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.

EXTENSION TO 2D PROBLEMS
In this section, we show how to straightforwardly apply the 1D concepts of the Taylor-based extrapolation method to 2D peridynamic problems. In 2D the fictitious nodes surround the real body as shown in Figure 26. In 2D models we can have corner nodes that do not exist in 1D. The numerical results of a 2D example are also shown hereinafter. This section is not meant to provide a complete treatment of the Taylor-based extrapolation method in 2D peridynamic models (which will be discussed in future developments), but just to show that the proposed method yields very accurate results also in 2D problems.

Taylor-based extrapolation method in 2D models
The proposed method is based on the nearest-node strategy, which entails the search for the nearest real node for each fictitious node. Clearly, all the real nodes found by this procedure are boundary nodes (i.e., real nodes closest to the physical boundary). Moreover, most of the fictitious nodes lies on the same line as their own closest boundary nodes, as shown F I G U R E 22 Relative percentage "errors" for the peridynamic models without any correction (in blue) and with the Taylor-based extrapolation method with orders n = 1, 2, 3 (respectively, in red, orange, and green) for the clamped bar under linearly distributed load F I G U R E 23 Boundary conditions for the clamped bar under sinusoidally distributed load, which is acting in the direction of the bar

F I G U R E 24
Numerical results obtained with the peridynamic models without any correction (in blue) and with the Taylor-based extrapolation method with orders n = 1, 3, 5 (respectively, in red, green, and yellow) for the clamped bar under sinusoidally distributed load. The dashed black line represents the reference solution derived analytically by the classical continuum mechanics F I G U R E 25 Relative percentage "errors" for the peridynamic models without any correction (in blue) and with the Taylor-based extrapolation method with orders n = 1, 3, 9 (respectively, in red, green, and yellow) for the clamped bar under sinusoidally distributed load in Figure 26(A). Therefore, the extrapolation for those fictitious nodes can be reduced to a 1D Taylor series expansion, which employes the same formulae already described in Section 4.1 (the same matrices can be used for displacements both in x and y direction). However, there are still the fictitious nodes near the corner of the body that require a 2D Taylor-based extrapolation procedure. We show here, briefly, the Taylor series expansion of a fictitious node i near the corner of the body (see Figure 26(B)) with a maximum order n = 2. This procedure can be easily generalized also for a higher order of the Taylor expansion. The displacements in x and y direction, named, respectively, u and v, of node i can be determined by the Taylor series expansion about the corner node j as: where x = x i − x j and y = y i − y j . The derivatives of the displacements of the corner node j can be found with a classical finite difference method in two dimensions (correspondent to the computation of the derivatives carried out in Section 4.1 for 1D models).
Remark 4. The 1D Taylor-based extrapolation procedure can be also exploited to impose Dirichlet boundary conditions by following the lines perpendicular to the boundary as shown in Figure 26(A) (please refer to Section 4.3 to see how to impose constraints by means of the 1D Taylor-based extrapolation). In particular, the matrices containing the factors of the Taylor series expansion (see Equation 58) can be used for the displacements both in x and y directions. Therefore, the procedure for imposing the constraints at the edges of the 2D body is the same used in 1D models, but applied twice: once to u and once to v. However, there is an ambiguity at the corner nodes since two lines perpendicular to the boundary can be followed, as shown by node j in Figure 26(A). In order to prevent the corner node from being overconstrained, we suggest to impose the Dirichlet boundary conditions only for the displacement perpendicular to the boundary (i.e., constrain the displacement in x direction for a boundary parallel to y-axis and the displacement in y direction for a boundary parallel to x-axis).

Clamped plate under uniformly distributed load
The advantages of using the proposed method in 2D models to reduce the surface effect near the boundaries can be shown with the following numerical example of a clamped plate under uniformly distributed load in x direction. The boundary conditions are shown in Figure 27 and the used data is reported in Table 4. As in the previous section, we solve the peridynamic problem first without adopting any corrections at the boundary and second by employing the proposed method with an increasing order n of the Taylor-based extrapolation. These numerical results are compared with the reference solution (see Figure  Δx FEM = Δy FEM = Δx∕2. The reference solution does not coincide with the peridynamic solution (see References 46,47) because numerical simulations must adopt a finite value for , but this is the closest available solution. Therefore, the percentage "errors" of the displacements are defined at each node i as where u PD i and v PD i are the displacements of node i in x and y direction computed with Peridynamics, u FEM i and v FEM i are the displacements of the same node obtained with the reference solution and u FEM is the displacement vector obtained with the reference solution at all the nodes.
The numerical results for the percentage "errors" computed with Equation (70) are shown in Figures 29 and 30. Note that the nodes near the corners are those with higher fluctuations. This is reasonable because the neighborhoods of the nodes near the corners are smaller than those of the nodes near the edges. Nonetheless, the numerical results are evidently improved by employing the Taylor-based extrapolation with respect to the case without corrections at the boundaries (except for v of only one node near the corner, see Figure 30(B)). Furthermore, the increasing of the order n for the Taylor-based extrapolation allows to mitigate more effectively the fluctuations near the edges and even near the corners of the body.

CONCLUSIONS
Two problems arising from dropping the infinite body assumption in Peridynamics have been addressed: the surface effect and the imposition of the boundary conditions. The surface effect in the 1D state-based peridynamic theory has been analyzed both analytically and numerically: the peridynamic points distant at least 2 from the boundary retrieve the same mechanical properties of the points in an infinite body, whereas those near the boundary experiences a stiffness fluctuation which generates a discrepancy with respect to the infinite body solution. The characteristic hardening/softening behavior from the bulk towards the boundary due to the surface effect is driven by two counteracting phenomena: the lack of some bonds near the boundary reduces the material stiffness but the decrease of the weighted volume of the points in the most external layer of the body partially compensate the previous phenomenon.
The other issue is that there is no standard procedure to impose boundary conditions in peridynamic models and the currently adopted strategies, derived often from classical mechanics concepts, imply some kind of approximation. The authors proposed a new version of the Taylor-based extrapolation method with the introduction of a fictitious layer in order to mitigate the discrepancy due to the surface effect. The displacements of the fictitious points are expressed as a function of the displacement of the boundary point and its derivatives via the truncated Taylor series expansion. Therefore, the fictitious layer keeps deforming as the real body near the boundary and completes in a rational way the partial neighborhoods near the boundary. Furthermore, the boundary conditions are implemented in a "peridynamic way": a Dirichlet boundary condition is imposed by constraining the boundary point and the fictitious layer corrects accordingly the stiffness of the real points near the boundary, and a Neumann boundary condition is implemented via the peridynamic concept of force flux, which involves the forces of all the fictitious bonds in the process.
The accuracy of the Taylor-based extrapolation method has been assessed by the several numerical examples, dealing with both 1D and 2D problems. The proposed method significantly improves the numerical results with respect to the peridynamic model without adopting any correction, even if low values of the Taylor expansion order n are chosen. Furthermore, it has been shown that the undesired fluctuation due to the aforementioned boundary issues can be made to decrease to a negligible level by increasing the order n.
Since bond-based Peridynamics is a special case of the state-based one, 2 the proposed method can be successfully used also in bond-based models. The authors will conduct further studies in order to complete the extension to the 2D case of the new Taylor-based extrapolation method by including the algorithm to impose Neumann boundary conditions.

ACKNOWLEDGMENTS
The term in square brackets is the desired force density scalar state, which can be rearranged through Equation (15) in a form which allows an easier computational implementation: where k = ( E)∕((1 − 2 )(1 + )) and k e = E∕(1 + ) are the newly defined constants, which were derived from the other peridynamic constants c ,k andĜ. Note that the dilatational term vanishes for = 0.

APPENDIX B. FORCE DENSITY DISTRIBUTION FOR THE NONLOCAL BOUNDARY CONDITION
The force density distribution b(x) needed to impose a constant deformation 11 = on a finite body of length l (see Figure 3) has been evaluated with Equations (4) and (22) solving the integrals in a way similar to that applied to the case of the force flux. The analytical result is reported in Table B1 and plotted in Figure B1 for different values of the horizon . This force density distribution also compensates the surface effect. As reasonably expected, the resultant of force density distribution is equal to the force which should be imposed at the boundary in classical mechanics. The obtained function changes slightly if another influence function is chosen, as one can verify with further numerical analyses.

APPENDIX C. ELIMINATION OF THE SURFACE EFFECT IN A FINITE BODY UNDER CUBIC DISPLACEMENT
Here, we aim to compute the force flux , defined in Equation (21), in an infinite body Ω ∞ under a smooth deformation, by following a procedure similar to that illustrated for the bond-based Peridynamics in appendix B.1 of Reference 46. Then, we compare the obtained solution with the solution derived with the proposed Taylor-based extrapolation method in a finite body Ω under a predefined displacement distribution.
In an infinite body Ω ∞ , the weighted volume m ∞ is constant in the whole domain and is equal to A 3 when the hyperbolic influence function = ∕| | is chosen (Equation 25). For later use, the dilatation ∞ of a point x in an infinite body Ω ∞ is computed by means of a Taylor series expansion of u(x ′ ) = u ′ about x: The dilatation converges to the classical mechanics dilatation cl if the third-order and higher derivatives of the displacement are negligible or if the horizon approaches 0.
Consider now the notation of Figure C1. The force flux (x ′′ ∈ Ω ∞ , +1) = ′′ ∞ in the infinite body is therefore computed as: where the values of the influence function and the weighted volume have been substituted. We proceed by solving the integrals separately.
where d k u ′ dx k and d k u ′′ dx k are the derivatives of the displacement computed at points x ′ and x ′′ , respectively. The second integral in Equation (C2) is computed by means of a Taylor expansion of the displacements u and u ′ about x ′′ :  . (C5) Note that the force flux converges to the classical stress 11 if the third-order and higher derivatives of the displacement are negligible or if the horizon approaches 0, as in References 46,47. Now we show that, by implementing the method of the Taylor-based extrapolation over the fictitious points, the peridynamic solution of the force flux in a finite body Ω under a predefined smooth deformation matches that of the infinite body. We choose to show the case in which the fourth-order or higher derivatives of the imposed displacement are equal to 0. The dilatation ∞ and the force flux ′′ ∞ in the infinite body Ω ∞ are therefore given as: Substituting Equations (C9) and (C11) in Equation (C8) and solving the integrals, the following result is obtained: