On the treatment of boundary conditions for bond-based peridynamic models
OOn the treatment of boundary conditions for bond-basedperidynamic models
Serge PrudhommeDepartment of Mathematics and Industrial Engineering, Polytechnique Montr´eal,C.P. 6079, succ. Centre-ville, Montr´eal, Qu´ebec H3C 3A7, [email protected] DiehlCenter for Computation and Technology, Louisiana State University,Digital Media Center, 340 E. Parker Blvd, Baton Rouge, LA 70803, USAhttp://orcid.org/0000-0003-3922-8419
Abstract
In this paper, we propose two approaches to ap-ply boundary conditions for bond-based peridynamicmodels. There has been in recent years a renewedinterest in the class of so-called non-local models,which include peridynamic models, for the simulationof structural mechanics problems as an alternativeapproach to classical local continuum models. How-ever, a major issue, which is often disregarded whendealing with this class of models, is concerned withthe manner by which boundary conditions should beprescribed. Our point of view here is that classi-cal boundary conditions, since applied on surfaces ofsolid bodies, are naturally associated with local mod-els. The paper describes two methods to incorporateclassical Dirichlet and Neumann boundary conditionsinto bond-based peridynamics. The first method con-sists in artificially extending the domain with a thinboundary layer over which the displacement field isrequired to behave as an odd function with respect tothe boundary points. The second method resorts tothe idea that peridynamic models and local modelsshould be compatible in the limit that the so-calledhorizon vanishes. The approach consists then in de-creasing the horizon from a constant value in the inte- rior of the domain to zero at the boundary so that onecan directly apply the classical boundary conditions.We present the continuous and discrete formulationsof the two methods and assess their performance onseveral numerical experiments dealing with the sim-ulation of a one-dimensional bar.
In recent years, novel non-local modeling theorieshave been developed for the simulation of problemsin structural and fracture mechanics. Peridynamictheory [1], which can be viewed as a generalization ofclassical continuum mechanics, provides for instancesuch a class of non-local models. Peridynamics doesnot involve differentials of displacement fields, whichmakes it an attractive framework for the modelingand simulation of fracture mechanics applications.Validation of the peridynamic theory has been under-taken using experimental data from a variety of wavepropagation applications as well as crack initiationand propagation experiments, see e.g. [2]. However,major issues identified in that review include the so-called surface or skin effect [3, 4, 5] and the manner bywhich boundary conditions should be modeled within1 a r X i v : . [ c s . C E ] A ug non-local model. This was recognized early on dur-ing the development of peridynamics as Silling andAskari in [6], for instance, explained that displace-ment and traction boundary conditions should beprescribed within a layer of finite thickness under thesurface of a solid body. Unfortunately, such bound-ary data cannot be assumed to be known in practicalapplications of interest. Several methods have sincebeen proposed, see e.g. [7, 8, 9, 10, 11, 12, 13], inorder to circumvent the lack of data and correct thespurious surface effects that are created as a result inthe vicinity of the boundaries. Starting from the factthat data is only available at the surface of a solidbody, we believe that boundary conditions for peri-dynamic models should be prescribed in a mannersimilar to that for local classical continuum models.We study in this paper two approaches to incorporateclassical Dirichlet and Neumann boundary conditionsinto a one-dimensional linearized bond-based peridy-namic model [1]. Our main objective in this paper isto mathematically analyze the resulting models, andmore specifically, to verify that their solutions con-verge to that of the compatible linear elasticity modelwith respect to the so-called peridynamic horizon anddiscretization parameters.The first approach, referred here to as the ex-tended domain method (EDM), relies on the ideato artificially extend the domain by a thin layer ofthe size of peridynamic horizon. The idea is notnew and is similar to that proposed by the authorsin [14, 15, 16, 17, 18, 19], among others. In [14, 15, 16]though, the thin layer was introduced within the solidbody and subjected to Dirichlet-like non-local bound-ary conditions, that is, the displacement field in thelayer was either arbitrarily set to zero or was con-strained to match the values of the exact solutionin the case of a manufactured problem. Zhou andDu [20] proposed to constrain the solutions to be ei-ther odd or even from either side of the boundary, inorder to mimic homogeneous Dirichlet or Neumannboundary conditions, respectively. We propose hereto constrain the displacement field in the artificiallayer to satisfy odd extensions with respect to theboundary points so that one can handle both Dirich-let and Neumann non-homogeneous boundary condi-tions. The second approach, referred here to as the vari-able horizon method (VHM), proposes to vary theperidynamic horizon from a constant value in the in-terior of the domain to zero as one approaches theboundary. The approach was considered in [3, 21]for the treatment of boundary conditions along withadaptive mesh refinement. The method is motivatedby the fact that the local model and non-local modelare consistent and compatible in the limit that thehorizon vanishes, so that one can justify the use ofclassical Dirichlet and Neumann boundary conditionsto be prescribed at the boundary of the domain.Varying the horizon size necessitates however a scal-ing of the material parameter for the peridynamicmodel. We propose here an approach similar to [3, 21]but with a scaling in 1D that appears to be differentfrom that suggested in [3]. We also mention that thevariable horizon approach was utilized in [15, 22, 23]for the treatment of interface problems in which dif-ferent peridynamic models or a continuum model anda peridynamic model are coupled together.Performance and comparison of the two methodsare investigated on several numerical examples in-volving a one-dimensional linear elastic bar. In par-ticular, we carry out a δ -convergence analysis [3, 16]by considering four test cases in which the manu-factured solutions are linear, quadratic, cubic, andquartic, in order to put in evidence the differencesand similarities between the methods. The EDM andVHM solutions are also compared with the solution ofthe local linear elasticity model (LLEM). All modelsare discretized here using the finite difference methodin order to allow for consistent comparison and areassessed using the maximum relative error with re-spect to the exact solution of the linear elastic localproblem. In all cases, we verify that the numericalrates of convergence at which the peridynamic so-lutions converge to the linear elasticity solution areconsistent with the theoretical rates. We also presentadditional examples in order to study the effect of acorrecting factor for the extended domain method, tocarry out a m -convergence analysis, and to considerthe case of a solution with a very steep gradient inthe vicinity of one boundary point.The paper is organized as follows. Following theintroduction, we present in Section 2 the model2roblem based on the classical linear elasticity the-ory and corresponding linearized bond-based peridy-namic model for the simulation of a one-dimensionalbar held fixed at one end and subjected to a trac-tion at the other end. We proceed in Section 3 withthe description of well-known spurious effects nearboundaries when using non-local models and assessa correction method to alleviate the issue. We intro-duce in Section 4 and Section 5 the continuous anddiscrete formulations of the extended domain methodand the variable horizon method, respectively, andassess their performance in Section 6 on several one-dimensional numerical examples. We complete thepaper in Section 7 with concluding remarks and per-spectives. The model problem consists in finding the static equi-librium of a bar, of unit length and constant cross-sectional area, held fixed at one end and subjectedto a longitudinal traction at the other end. We shallsuppose here that the deformations in the bar are in-finitesimally small and can be adequately describedby the theory of linear elasticity. We present belowthe local model based on classical elasticity theoryand the non-local model counterpart built upon thelinearized bond-based peridynamic theory [1].
Let Ω = (0 , u in the bar such that: − EAu (cid:48)(cid:48) = f b , ∀ x ∈ Ω , (1) u = 0 , at x = 0 , (2) EAu (cid:48) = g, at x = 1 , (3)where E and A are the constant modulus of elastic-ity and cross-sectional area of the bar, respectively, g ∈ R is the traction force applied at end point x = 1,and the scalar function f b = f b ( x ) is the externalbody force density (per unit length). We suppose here that f b is chosen sufficiently smooth so that reg-ularity in the solution is not an issue when comparingthe solutions from the two models. Note that we usethe notation u for the solution to the classical elas-ticity problem in order to emphasize that it may bedifferent from the peridynamic solution u introducedbelow. Let δ > H δ ( x ) = ( x − δ, x + δ ) bethe subdomain of the neighboring particles within thehorizon. We observe that for any given point x in theinterval Ω δ = ( δ, − δ ), we have that H δ ( x ) ⊂ Ω. Wealso introduce the boundary domains B δ = (0 , δ ) and B δ = (1 − δ,
1) such that Ω = B δ ∪ Ω δ ∪ B δ , as showin Figure 1. The general formulation, in one or higherdimension, of the linearized microelastic model [1] isgiven by: (cid:90) H δ ( x ) ∩ Ω κ ξ ⊗ ξ (cid:107) ξ (cid:107) ( u ( y ) − u ( x )) dy + f b ( x ) = 0 , (4)where κ is the parameter that characterizes the stiff-ness of the “bonds” between point x and the neigh-boring points y ∈ H δ ( x ), ξ is the vector betweentwo material points in the reference configuration, i.e. ξ = y − x , (cid:107) ξ (cid:107) is the Euclidean norm of vector ξ , and u ( x ) is the displacement of x in the deformed config-uration. In the case of the one-dimensional bar, theabove integral at a point x ∈ Ω δ can be rewritten as: (cid:90) x + δx − δ κ u ( y ) − u ( x ) | y − x | dy + f b ( x ) = 0 . (5)Ideally, the material parameter κ could be identi-fied from experimental data. However, one may beinterested in the value of κ such that the solution ofthe linearized microelastic brittle peridynamic modelis the same as that of the continuum local model.One may do so by matching the strain energy of thetwo models [1]. Alternatively, one can try to recoverequation (1) from (4) when taking the limit δ → u ( x ) is sufficiently3 δ x − δ x x + δ − δ δ H δ ( x ) B δ B δ Figure 1: Definition of the computational domainsfor the peridynamic model.smooth, one can write, for all y ∈ Ω δ : u ( y ) − u ( x ) = u (cid:48) ( x )( y − x ) + 12 u (cid:48)(cid:48) ( x )( y − x ) + 16 u (cid:48)(cid:48)(cid:48) ( x )( y − x ) + O ( | y − x | ) , (6)so that the integral in (4) reduces to: (cid:90) x + δx − δ κ u ( y ) − u ( x ) | y − x | dy = κδ u (cid:48)(cid:48) ( x ) + O ( δ )) . (7)Substituting the integral in (4) with the above resultyields: − κδ u (cid:48)(cid:48) + O ( δ )) = f b , ∀ x ∈ Ω δ . (8)By taking the limit when δ →
0, one then recoversthe differential equation (1) pointwise whenever κ ischosen as: κδ EA, that is κ = 2 EAδ , (9)in agreement with [24]. In other words, the peridy-namic model can be used to describe the behavior inthe bar as one is able to identify the value of parame-ter κ that makes it compatible to the linear elasticitymodel. The main and well known issue with peridy-namics is whether the model can also be employed inthe region Ω \ Ω δ , i.e. in the region of size δ along theboundary. A solution to this issue could be to trans-pose the Dirichlet condition (2) and Neumann condi-tion (3), prescribed at the boundaries of the bar, intovolumetric boundary conditions to constrain the ma-terial in B δ and B δ . However, it is still unclear how this can be achieved, and if it were, how to identifythe corresponding boundary data. The point of viewthat we will follow here is that the notion of bound-ary conditions, which are prescribed on the surface ofsolid bodies, is naturally associated with local mod-els. Our goal in the following section is to presentsome approaches to adapt the Dirichlet and Neumannconditions to the peridynamic problem. In this section, we study the integral in (4) at a pointin B δ , i.e. near the left boundary of the bar. Theanalysis is similar for a point in B δ near the otherboundary. Let x ∈ B δ . Since H δ ( x ) extends beyondΩ, the integral becomes: (cid:90) x + δ κ u ( y ) − u ( x ) | y − x | dy = (cid:90) δ − x κ (cid:16) u (cid:48) ( x ) ξ | ξ | + u (cid:48)(cid:48) ( x )2 ξ | ξ | + u (cid:48)(cid:48)(cid:48) ( x )6 ξ | ξ | + . . . (cid:17) dξ = κu (cid:48) ( x ) (cid:90) δ − x ξ | ξ | dξ + κ u (cid:48)(cid:48) ( x ) (cid:90) δ − x | ξ | dξ + κ u (cid:48)(cid:48)(cid:48) ( x ) (cid:90) δ − x | ξ | ξdξ + . . . = κ ( δ − x ) u (cid:48) ( x ) + κ δ + x u (cid:48)(cid:48) ( x )+ κ δ − x u (cid:48)(cid:48)(cid:48) ( x ) + . . . (10)where we have used the change of variable ξ = y − x . When comparing the above result to (7), it isclear that truncation of the set H δ ( x ) close to theboundary induces forces involving the first derivativeof u . This term vanishes in the case of free surfacessince zero traction, i.e. EAu (cid:48) = 0, is prescribed on theboundary. Moreover, the coefficient of u (cid:48)(cid:48) is differentfrom that in (7). Since 0 ≤ x ≤ δ , the coefficient isactually smaller than for points in the interior of thedomain, meaning that the material becomes softer asone approaches the free boundary, which is referredto as the so-called skin effect [3].4everal methods to correct the skin effect, suchas the volume correction method [25, 26], the forcedensity method [27, 18], the energy method [27, 18],the force normalization method [28], or the position-aware method [29] have been compared in [5]. Wejust note that some of the methods seem somewhatheuristic and that the force normalization approachof Macek and Silling [28] appears to be valid only fortwo- and three-dimensional peridynamic models, asit involves an integral that is singular in the case ofone-dimensional problems.Looking at (10), a simple approach would be toreplace the parameter κ by a position-dependent pa-rameter ¯ κ ( x ) such that the coefficient of the secondderivative matches that obtained for a point insideΩ δ , see (7), or equivalently, allows one to recover thelinear elasticity model, that is:¯ κ ( x )2 δ + x κδ EA, which leads to: ¯ κ ( x ) = κ δ δ + x . (11)It is notable that this result is equivalent to the cor-rection factor that one obtains by the volume correc-tion method. In fact, one would recover the sameresult by the energy method, at least in 1D. Indeed,the method consists in identifying the parameter ¯ κ ( x )such that the strain energy computed at a point x near the boundary matches that computed at a pointinside Ω δ for a given deformation. The strain energyper unit length of a “bond” is given by: ω ( u ) = ¯ κ ( x )2 ( u ( y ) − u ( x )) | y − x | . It follows that the strain energy at a point x ∈ B δ isprovided by: (cid:90) x + δ ¯ κ ( x )2 ( u ( y ) − u ( x )) | y − x | dy ≈ ¯ κ ( x )2 [ u (cid:48) ( x )] (cid:90) δ − x | ξ | dξ ≈ ¯ κ ( x )4 [ u (cid:48) ( x )] ( δ + x ) . For a point inside the domain, i.e. x ∈ Ω δ , the strainenergy reads: (cid:90) x + δx − δ κ u ( y ) − u ( x )) | y − x | dy ≈ κ u (cid:48) ( x )] (cid:90) δ − δ | ξ | dξ ≈ κ u (cid:48) ( x )] δ . Matching those two energies for a given deformation u (cid:48) naturally leads to (11).Replacing κ by ¯ κ ( x ) in (10) and using (11), oneobtains the new expression of the integral: (cid:90) x + δ κ u ( y ) − u ( x ) | y − x | dy = κ δ ( δ − x ) δ + x u (cid:48) ( x ) + κδ u (cid:48)(cid:48) ( x )+ κ δ ( δ − x ) δ + x u (cid:48)(cid:48)(cid:48) ( x ) + . . . = κδ (cid:18) − ( x/δ )1 + ( x/δ ) δ u (cid:48) ( x ) + u (cid:48)(cid:48) ( x )+ 1 − ( x/δ ) x/δ ) δ u (cid:48)(cid:48)(cid:48) ( x ) + . . . (cid:19) (12)Therefore, if u (cid:48) ( x ) = 0, such as near free surfaces, theperidynamic model leads to: − κδ u (cid:48)(cid:48) + O ( δ )) = f b , ∀ x ∈ B δ . (13)In other words, due to the loss of symmetry of theintegration domain, the peridynamic model providesan approximation of the classical linear elasticity the-ory of order one in δ in the vicinity of free surfaces, tobe compared to O ( δ ) in the interior of the domain,see (8).Nevertheless, we emphasize here that the abovecorrection approach is still unsuitable for prob-lems involving Dirichlet boundary conditions or non-homogeneous Neumann boundary conditions, sincethe first derivative of the solution does not necessarilyvanish near the boundaries in these cases. We presentin the next two sections alternative approaches toalleviate boundary effects, namely the extended do-main method, if the horizon must be kept constanteverywhere, and the variable horizon method, if thesize of the horizon is allowed to vary in the domain.5 Extended domain method
Let (cid:98) Ω δ = [ − δ, δ ] ⊃ Ω be the domain Ω augmentedby the boundary domains D δ = [ − δ,
0) on the left-hand side and D δ = (1 , δ ] on the right-hand side,see Figure 2. The goal here is to construct an exten-sion of the solution u in Ω to the whole domain (cid:98) Ω δ .In the previous section, the integral in (10) can beviewed as the integral over [ x − δ, x + δ ] obtained byextending u to u ( x ) = 0, ∀ x ∈ D δ . The issue withsuch an extension is that it introduces a term involv-ing the first derivative arising from the discontinuityin u (cid:48) at x = 0. We propose here, given a continuousfunction u in Ω, to consider odd extensions of u withrespect to x = 0 and x = 1 as follows: u ( x ) = 2 u (0) − u ( − x ) , ∀ x ∈ D δ , (14) u ( x ) = 2 u (1) − u (2 − x ) , ∀ x ∈ D δ . (15)We first note that the above conditions imply thatthe function u is continuous at x = 0 and x = 1 andsatisfies the following properties: u ( n ) ( x − ) = u ( n ) ( x + ) , for odd n > ,u ( n ) ( x − ) = − u ( n ) ( x + ) , for even n > , (16)where x = 0 or 1, u ( n ) ( x − ) and u ( n ) ( x + ) denote thelimits of u n when one approaches x from the leftand from the right, respectively. In other words, theodd derivatives of u are continuous at x = 0 and x = 1 while the even derivatives are opposite at thosepoints.Using the above properties and Taylor expansionsof u at x = 0, one can show, following lengthy calcu-lations, that for x ∈ (0 , δ ): (cid:90) x + δx − δ κ u ( y ) − u ( x ) | y − x | dy = κδ (cid:20)(cid:18) xδ − (cid:16) − xδ (cid:17) x δ (cid:19) u (cid:48)(cid:48) (0 + ) + xu (cid:48)(cid:48)(cid:48) (0) + . . . (cid:21) = κδ (cid:20)(cid:18) xδ − (cid:16) − xδ (cid:17) x δ (cid:19) u (cid:48)(cid:48) ( x ) + O ( δ ) (cid:21) (17) − δ x − δ x x + δ δ (cid:98) Ω δ Ω H δ ( x ) D δ D δ Figure 2: Definition of domains for the extended do-main approach.where we have used the facts that u (cid:48)(cid:48) (0 + ) = u (cid:48)(cid:48) ( x ) − xu (cid:48)(cid:48)(cid:48) ( x ) + . . . and 0 < x < δ . This is a remarkableresult which shows that the leading term is now aterm that involves the second derivative. The termcontaining the first derivative vanishes here thanksto the fact that the first derivative of the extendedfunction is continuous at x = 0. If one had chosen aneven extension, this term would not have cancelled.Moreover, the peridynamic model converges to theclassical linear elasticity model as δ tends to zero ifthe material property κ is replaced, for x ∈ (0 , δ ), by:¯ κ ( x ) = κ (cid:18) xδ − (cid:16) − xδ (cid:17) x δ (cid:19) − , (18)and, similarly for x ∈ (1 − δ, κ ( x ) = κ (cid:18) − xδ − (cid:16) − − xδ (cid:17) (1 − x ) δ (cid:19) − . (19)However, the convergence is only linear in δ . As anexample, if x/δ = 1 /
2, one has ¯ κ ≈ k/ .
9, that is¯ κ ≈ . κ .We now formulate the peridynamic problem for theextended domain method. For the sake of simplicityin the presentation, we shall keep the value of κ con-stant over the whole domain, just keeping in mindthat its value can be corrected in the intervals (0 , δ )and (1 − δ,
1) according to (18) and (19), respectively.The problem consists then in finding u ∈ (cid:98) Ω δ such6hat: − (cid:90) x + δx − δ κ u ( y ) − u ( x ) | y − x | dy = f b ( x ) , ∀ x ∈ Ω , (20) u ( x ) − u (0) + u ( − x ) = 0 , ∀ x ∈ D δ , (21) u ( x ) − u (1) + u (2 − x ) = 0 , ∀ x ∈ D δ , (22) u = 0 , at x = 0 , (23) κδ u (cid:48) = g, at x = 1 , (24)where (9) has been used to modify the Neumannboundary condition (3) to (24). The fact that thedisplacement is extended in the boundary domainsby odd extensions allows one to prescribe Dirich-let or Neumann boundary conditions. Indeed, ifthe Neumann condition at x = 1 is replaced by anon-homogeneous Dirichlet boundary condition, onewould simply substitute (24) by: u = u , at x = 1 , where u ∈ R is a prescribed displacement.Although the solution to the problem (20)-(24)converges to the solution of the local linear elasticitymodel as the horizon δ approaches zero, the formula-tion nevertheless presents several drawbacks:A the problem is solved on a larger domain thanthe domain occupied by the bar, which increasesthe number of degrees of freedom for a given dis-cretization parameter h . This is in particulartrue when δ/h is large and in higher dimensions;B odd extension of the displacement field outsideof Ω may become difficult to carry on in twoand three dimensions for domains with complexgeometries, especially when boundaries involvecorners and edges;C conceptually, there may still be an ambiguity inthe application of the Neumann boundary condi-tion, which is derived from the local model, whilethe behavior in the vicinity of the boundary isgoverned by the non-local model.The approach proposed in Section 5 is an attempt atcircumventing those issues. We describe here the discretization of Problem (20)-(24) using a collocation approach. For a given δ , weintroduce a uniform grid spacing h chosen such that,as it is customarily done in the literature [6, 25], δ is a multiple of h , i.e. δ/h = m , with m a positiveinteger. For the sake of simplicity, but without lossof generality, we shall choose here m = 2. In thiscase, choosing h = 1 /n , with n a positive integer, thediscretization of domain (cid:98) Ω δ is determined by the gridpoints: x i = ih, i = − , − , , . . . , n, n + 1 , n + 2 , (25)where the points x − , x − , x correspond to the gridin D δ and the points x n , x n +1 , x n +2 to the one in D δ . Moreover, let u i , i = − , − , . . . , n + 2, denotethe discrete displacement at the grid points x i .Approximation of the integral in (20) is obtainedby classical quadrature formula using the grid points x i . Note that other more advanced quadrature rules,especially developed for non-local models, could al-ternatively be used, see e.g. [30, 31]. We then have,using the second-order trapezoidal integration rule, κδ h (cid:20) − u i − − u i − + 54 u i − u i +1 − u i +2 (cid:21) = f b ( x i ) , i = 1 , . . . , n − , (26)or, in a more compact form, with α = κδ / (16 h ), − αu i − − αu i − + 10 αu i − αu i +1 − αu i +2 = f b ( x i ) , i = 1 , . . . , n − . (27)The constraints (21) and (22) to enforce odd exten-sions on the displacements read: u − − u + u = 0 , (28) u − − u + u = 0 , (29) u n − − u n + u n +1 = 0 , (30) u n − − u n + u n +2 = 0 , (31)while the Dirichlet boundary condition at x = 0 im-plies: u = 0 . x = 1. Several options based on the fi-nite difference method are available. However, sincethe function u is only C at x = 1, we cannot considera second-order centered difference of the first deriva-tive. We therefore opt for the second-order one-sidedfinite difference: u (cid:48) (1) = u n − − u n − + 3 u n h + O ( h ) , that, in order to preserve the symmetric feature ofthe peridynamic model, can be rewritten using (30)and (31) in the centered form as: u (cid:48) (1) ≈ u n − − u n − + 6 u n h ≈ u n − − u n − + 6 u n − u n − + u n − h ≈ u n − − u n − + 4 u n +1 − u n +2 h . We emphasize that this is a second-order finite differ-ence of u (cid:48) (1), whether u satisfies the odd extensionor not. We propose to employ the latter formula todiscretize the Neumann boundary condition, i.e. κδ (cid:20) u n − − u n − + 4 u n +1 − u n +2 h (cid:21) = g, which can be rewritten as:2 αhu n − − αhu n − + 8 αhu n +1 − αhu n +2 = g with α = κδ / (16 h ) as before.The resulting system of linear equations associ-ated with the discretization of the extended domainmethod is summarized in Figure 3. Upon inspec-tion, it is worth noting that the use of the Dirichletboundary condition at x = 0 allows one to eliminatethe first equation of the system as the value of the de-gree of freedom u − , associated with the grid point x − in the extended domain, does not influence theother degrees of freedom. The same conclusion wouldalso hold for u n +2 if one had considered a Dirichletboundary condition at x = 1. We actually presentin the next subsection an approach to further reducethe system of equations. In fact, it is possible in the one-dimensional case toimplicitly enforce the odd extension constraints so asto reduce the system of equations and avoid solvingfor u − , u − , u n +1 , and u n +2 . Using (29) and (30),the fourth and ( n − κδ (cid:20) − u + 11 u − u − u h (cid:21) = f b ( x ) , (32) κδ (cid:20) − u n − − u n − + 11 u n − − u n h (cid:21) = f b ( x n − ) , (33)whereas using the left-sided second-order approxima-tion of u (cid:48) (1) to discretize the Neumann boundarycondition gives: κδ (cid:20) u n − − u n − + 3 u n h (cid:21) = g. The reduced system of equations is shown in Fig-ure 4. We emphasize here that the extended domainmethod and its reduced version provide exactly thesame solutions for i = 0 , . . . , n . We will therefore usethe reduced system displayed in Figure 4 for solv-ing the problems in the numerical examples. Weobserve that the matrix has in fact a similar struc-ture as the one obtained using the second-order fi-nite difference method for the discretization of thelocal model (see system of equations in Figure 7).We verify below that the equations associated withrows i = 1 , . . . , n − h , and a fortiori in δ since h = δ/ u (cid:48)(cid:48) at the points x i .Starting with (26), we see that the expression inbrackets, multiplied by the factor 1 /h , can be rewrit-ten as: − u i − − u i − + 10 u i − u i +1 − u i +2 h = − (cid:20) u i − − u i + u i +2 (2 h ) (cid:21) − (cid:20) u i − − u i + u i +1 h (cid:21) = − (cid:20) u (cid:48)(cid:48) ( x i ) + O ((2 h ) ) (cid:21) − (cid:20) u (cid:48)(cid:48) ( x i ) + O ( h ) (cid:21) = − u (cid:48)(cid:48) ( x i ) + O ( h ) . (34)8 − − · · · − α − α α − α − α − α − α α − α − α − α − α α − α − α − α − α α − α − α
00 0 0 2 αh − αh αh − αh − − u − u − u u u ... u n − u n − u n u n +1 u n +2 = f b ( x ) f b ( x )... f b ( x n − ) f b ( x n − ) g Figure 3: System of linear equations associated with the discretization of the extended domain methodwhere α = κδ h . · · ·− α α − α − α − α − α α − α − α − α − α α − α − α − α − α α − α αh − αh αh u u u ... u n − u n − u n = f b ( x ) f b ( x )... f b ( x n − ) f b ( x n − ) g Figure 4: Reduced system of linear equations associated with the discretization of the extended domainmethod where α = κδ h .In other words, the second derivative u (cid:48)(cid:48) ( x i ) is ap-proximated by the weighted average of the second-order centered differences of the second derivative us-ing, on one hand, the nearest neighbors and, on theother hand, the next-to-nearest neighbors.We now turn our attention to (32). The term inbrackets, using classical Taylor expansions, becomes: − u + 11 u − u − u h = − u (cid:48)(cid:48) ( x ) + O ( h ) , (35)where the details of the calculation have been omit-ted. This result was to be expected in view of (17).Since x /δ = 1 / / .
875 ac-tually corresponds to an approximation of the fac- tor κ/ ¯ κ ( x ) ≈ . / ≈ .
14 in order to recover − u (cid:48)(cid:48) ( x ), orconsider the material property at that point to be¯ κ ( x ) ≈ . κ according to (18), but we note thatthe approach would still remain suboptimal due tothe fact that one gets only a first-order approxima-tion. Note that the same conclusion holds in the caseof (33). We will illustrate the spurious perturbationsnear the boundaries in the numerical examples de-pending on the chosen correction method.9 Variable horizon method
Our motivation in this approach is to address some ofthe ambiguities associated with the extended domainapproach, the first one being that the problem is de-fined on a larger domain than the original domain Ωoccupied by the bar (although we have shown thanthe problem from the extended domain method canbe recast, in the one-dimensional case, as a reducedproblem defined on Ω). Our objective is therefore tointroduce a new problem restricted to the computa-tional domain Ω. Secondly, in order to avoid con-structing corrected forces near the boundaries, theset H δ ( x ) should always be a subset of domain Ω,suggesting that the horizon should tend to zero asone approaches the boundary. Finally, we supposethat one is interested in using a non-local model withconstant horizon δ in the interior of Ω sufficiently farfrom the boundaries. In order to satisfy all these re-quirements, one needs to construct a non-local modelwith variable horizon δ v ( x ) such that, ∀ x ∈ Ω:0 ≤ δ v ( x ) ≤ δ,x − δ v ( x ) ≥ ,x + δ v ( x ) ≤ . (36)The definition of the function δ v ( x ) is therefore notunique. However, the simplest continuous functionthat fulfill the above requirements is the piecewiselinear function defined as: δ v ( x ) = x, < x ≤ δ,δ, δ < x ≤ − δ, − x, − δ < x < . (37)The function δ v ( x ) is shown in Figure 5. Alter-native functions, in particular smoother functions,could also be considered as long as they take val-ues between 0 and δ v ( x ) (as defined in (37)) for all x ∈ Ω. In view of (9), the horizon being a function of x implies that the material parameter κ should alsodepends on x , i.e. κ = ¯ κ ( x ), if the non-local modelis to be compatible with the local model. The new δ v δ − δ x Figure 5: Example of a variable horizon function δ v ( x ). The circles centered at points x ∈ B δ = (0 , δ )are introduced to represent the associated domain H δ ( x ).problem consists in searching for u ( x ) such that: − (cid:90) x + δ v ( x ) x − δ v ( x ) ¯ κ ( x ) u ( y ) − u ( x ) | y − x | dy = f b ( x ) , ∀ x ∈ Ω , (38) u = 0 , at x = 0 , (39) EAu (cid:48) = g, at x = 1 . (40)In order for the model to be compatible with thelinear elasticity model everywhere in the domain witha convergence of order O ( δ ), the product ¯ κ ( x ) δ v ( x )needs to remain constant for all x ∈ Ω, so that¯ κ ( x ) δ v ( x ) = κδ , ∀ x ∈ Ω , (41)where κ satisfies (9). We note that this scaling of thematerial parameter appears to be different from thatin [3]. As before, we shall use a collocation approach to dis-cretize the problem (38)-(40) and choose the param-eter h such that δ/h = m , with m = 2. The grid inΩ is defined by the grid points: x i = ih, i = 0 , , . . . , n, (42)and u i , i = 0 , , . . . , n , denote the discrete displace-ments at the grid points x i . From (37), we observe10hat: δ v ( x i ) = δ/ h, i = 1 ,δ = 2 h, i = 2 , . . . , n − ,δ/ h, i = n − . The approximation of (38) at x and x n − by asecond-order integration scheme for the integral isgiven as:¯ κ ( x i )2 [ − u i − + 2 u i − u i +1 ] = f b ( x i ) , i = 1 , n − . However, using relation (41) and the fact that, inthose two cases, δ v ( x i ) = h , one can rewrite the aboveequation as: κδ h [ − u i − + 2 u i − u i +1 ] = f b ( x i ) , i = 1 , n − , or, using again the coefficient α = κδ / (16 h ), − αu i − + 16 αu i − αu i +1 = f b ( x i ) , i = 1 , n − . In the case of the remaining interior grid points x i , i = 2 , . . . , n −
2, the approximation of (38) is givenby (27), that is: − αu i − − αu i − + 10 αu i − αu i +1 − αu i +2 = f b ( x i ) , i = 2 , . . . , n − . The Dirichlet boundary condition is prescribed asusual as: u = 0 . Finally, it remains to discretize the Neumann bound-ary condition. A second-order approximation ofthe first derivative using a left-sided finite differencescheme is: u (cid:48) (1) = u n − − u n − + 3 u n h + O ( h ) . Therefore, the approximation of the Neumann condi-tion, using EA = κδ / α , reads:4 αhu n − − αhu n − + 12 αhu n = g. The resulting system of linear equations, corre-sponding to the discretization of the variable horizon method, is summarized in Figure 6. We note herethat the coefficient α could be factorized in front ofthe matrix. We conclude the section by emphasiz-ing that the variable horizon method provides an ap-proximation of the classical linear elasticity model oforder two in the parameter h , and a fortiori, in δ . Wewill see that this rate of convergence is consistentlyobserved in the numerical examples presented below. The objective of this section is to present several nu-merical examples in order to compare the solutions ofthe extended domain and variable horizon approachesfor implementing boundary conditions in non-localperidynamic models. For completeness, we will alsoshow the results obtained with the local linear elas-ticity model. The system of equations for the dis-cretization of the problem (1)-(3) by a second-orderfinite difference method is shown in Figure 7, wherewe have assumed the same grid on Ω than the one de-fined in Section 5.2. For simplicity, but without lossof generality, we shall consider a one-dimensional barfor which EA = 1, so that κ = 2 /δ , α = 1 / (8 h ),and β = 1 / (2 h ). Moreover, the value of the bound-ary data for the Neumann boundary condition is setin all examples to unity, i.e. g = 1, unless explicitlysaid otherwise. The first series of examples will beconcerned with δ -convergence (see e.g. [3, 16]), forwhich we shall measure the error in the peridynamicsolutions with respect to the classical local solution.In this case, we note that the discretization param-eter h implicitly goes to zero as well since h = δ/ m -convergence, for which δ is fixedand h goes to zero, in order to show that the ap-proximate solution actually converges to the corre-sponding peridynamic solution. Finally, we consideran example whose solution is given in terms of an ex-ponential function and exhibits a very steep gradientat one end point. The objective is to show that, asexpected, the rates of convergence are unaltered evenin the presence of large variations in the solution.11 · · ·− α α − α − α − α α − α − α − α − α α − α − α − α α − α αh − αh αh u u u ... u n − u n − u n = f b ( x ) f b ( x )... f b ( x n − ) f b ( x n − ) g Figure 6: System of linear equations associated with the discretization of the variable horizon method where α = κδ h . · · ·− β β − β − β β − β − β β − β
00 0 0 − β β − β βh − βh βh u u u ... u n − u n − u n = f b ( x ) f b ( x )... f b ( x n − ) f b ( x n − ) g Figure 7: System of linear equations associated with the finite difference discretization of the local linearelasticity model with β = EA h = 4 α . 12 .1 δ -convergence numerical study We consider in this section four examples for whichthe manufactured solutions are polynomial functionsof degree up to four. The objectives of this numeri-cal study are to analyze the rates of convergence in δ as the peridynamic solutions converge to the linearelasticity solution and to highlight the differences be-tween the methods. For comparison, we compute themaximum relative error: E n = max i =1 ,...,n (cid:12)(cid:12)(cid:12)(cid:12) u ( x i ) − u i u ( x i ) (cid:12)(cid:12)(cid:12)(cid:12) , with n the total number of degrees of freedom of thenumerical solution u i , i = 1 , . . . , n , for the local lin-ear elasticity model (LLEM), the extended domainmethod (EDM) without correction, and the variablehorizon method (VHM). Linear solution
The first example that we shall consider is definedwhen the body force density f b identically vanisheson Ω, i.e. f b ( x ) = 0, ∀ x ∈ Ω, in which case thesolution to Equations (1)-(3) is the linear function: u ( x ) = gx, ∀ x ∈ ¯Ω = [0 , . It is not difficult to show that the linear discrete func-tion u i = igh , i = 0 , . . . , n is the solution to the sys-tems of equations given in Figures 4, 6, and 7, for anyvalue of h with h = 1 /n . In other words, the linearfunction is exactly reproduced by the discrete solu-tions of the local linear elasticity model (LLEM), theextended domain method (EDM), and the variablehorizon method (VHM). Quadratic solution
The second example deals with the problem for whichthe body force density is given by: f b ( x ) = 1 . In that case, the exact solution to the local modelequations (1)-(3) reads: u ( x ) = x (4 − x )2 . Maximum relative error E n n δ LLEM EDM VHM4 0.5 0 0.05098 08 0.25 0 0.02443 016 0.125 0 0.01202 032 0.0625 0 0.00596 0Table 1: Maximum relative error E n for the quadraticsolution obtained with the local linear elastic-ity model (LLEM), the extended domain method(EDM), and the variable horizon method (VHM).We now solve the discrete problems with δ = 1 / /
4, 1 /
8, and 1 /
16, (i.e. h = 1 /
4, 1 /
8, 1 /
16, and 1 / h , as expected. Indeed,since the derivatives of order three and higher of theexact solution u all vanish, the discrete methods areable in this case to exactly reproduce the solution.On the other hand, we observe that the errors inthe displacement when using EDM are nonzero andseem to decrease with a rate of convergence of or-der one with respect to h . In this case, the onlysources of error arise from the fact that the secondderivative is not properly estimated at x = h and x n − = 1 − h , as explained by the result in Eq. (35).As a matter of fact, the error becomes zero if we mul-tiply the left-hand side in the second and ( n − /
7. Fi-nally, we show in Figure 8 the error e i = u ( x i ) − u i , i = 1 , . . . , n , in order to illustrate how the sources oferror committed at x and x n − pollute the solutionin the whole computational domain. Cubic solution
We repeat the previous experiment with the followingbody force density: f b ( x ) = x, . . . . . . x − . − . − . − . − . − . . E rr o r i nd i s p l a ce m e n t Example with Quadratic Solution using EDM δ =0.5 δ =0.25 δ =0.125 δ =0.0625 Figure 8: Error in displacement for the quadratic so-lution using the extended domain method (EDM) forvarious values of horizon δ .for which the exact solution to Equations (1)-(3) isgiven by: u ( x ) = x (3 − x )(3 + x )6 . The results are reported in Table 2. In the case ofLLEM and VHM, we observe that the maximum rel-ative errors decrease with a rate of convergence oforder two with respect to h , as expected. Moreover,the errors for both methods are identical. This is ex-plained by the fact that the source of error only arisesfrom the approximation of the Neumann boundarycondition. Indeed, the approximation being of ordertwo, the leading term of the truncation error involvesthe third derivative of u , which is nonzero, while theleading term of the truncation error for the approx-imation of the second derivative involves the fourthderivative of u , which is in this case zero. There-fore, since the Neumann boundary condition is dis-cretized in the same way for LLEM and VHM, thetwo methods provide the same approximate solutionof the problem, hence the same relative error.As far as EDM is concerned, we note the same rateof convergence, approximately of order one, as in theprevious problem with the quadratic solution. It isonce again obvious from Figure 9 that the issue comesfrom the approximation of the second derivative at x and x n − . Indeed, we observe a clear change in the Maximum relative error E n n δ LLEM EDM VHM4 0.5 0.01481 0.02774 0.014818 0.25 0.00379 0.01764 0.0037916 0.125 0.00096 0.00983 0.0009632 0.0625 0.00024 0.00517 0.00024Table 2: Maximum relative error E n for the cubic so-lution obtained with the local linear elasticity model(LLEM), the extended domain method (EDM), andthe variable horizon method (VHM). . . . . . . x − . − . − . − . − . − . . E rr o r i nd i s p l a ce m e n t Example with Cubic Solution using EDM δ =0.5 δ =0.25 δ =0.125 δ =0.0625 Figure 9: Error in displacement for the cubic solu-tion using the extended domain method (EDM) forvarious values of horizon δ .slope of the displacement error at x n − . The issue ishowever less visible at x as the solution is close tozero near the boundary. Quartic solution
We repeat the previous experiment and consider thecase of the quadratic body force density: f b ( x ) = x , which provides the following quartic solution toEquations (1)-(3): u ( x ) = x (16 − x )12 . E n n δ LLEM EDM VHM4 0.5 0.03226 0.00663 0.036178 0.25 0.00891 0.01267 0.0108516 0.125 0.00233 0.00889 0.0029432 0.0625 0.00060 0.00511 0.00076Table 3: Maximum relative error E n for the quar-tic solution obtained with the local linear elastic-ity model (LLEM), the extended domain method(EDM), and the variable horizon method (VHM).The maximum relative errors are shown in Table 3.We first note that the relative error for both LLEMand VHM decreases with a rate of convergence oforder two similarly to the previous example. How-ever, the errors in the solutions obtained with VHMare slightly greater than those in the solutions ob-tained by LLEM. Since the fourth derivative of u isno longer zero in Ω, the error is now partly due to theapproximation of the second derivative. In the caseof VHM, it is obtained as the average of two centereddifference stencils of order O ( h ) and O ((2 h ) ), re-spectively (see e.g. Eq. (34)) at the grid points x i , i = 2 , . . . , n −
2. It naturally follows that the solutionis slightly less accurate than the solution obtainedfrom LLEM, which only involves approximations oforder O ( h ). However, the difference between the er-rors in the two solutions correspond to the modelingerror of order δ that arises from using the non-localmodel instead of the local model.Finally, we show in Figure 10 the displacement er-ror for the EDM solution. The conclusions that wedrew in the cases of the quadratic and cubic solu-tions also hold, in addition to the fact that the be-havior of the error seems even more unpredictablehere. For comparison, we show in Figure 11 the dis-placement error in the VHM solution, which exhibitsa monotonic and smoother behavior and decreasesmuch faster as δ goes to zero. In the previous numerical examples, we considered . . . . . . x − . − . − . − . − . − . − . . E rr o r i nd i s p l a ce m e n t Example with Quartic Solution using EDM δ =0.5 δ =0.25 δ =0.125 δ =0.0625 Figure 10: Error in displacement for the quartic so-lution using the extended domain method (EDM) forvarious values of horizon δ . . . . . . . x . . . . . E rr o r i nd i s p l a ce m e n t Example with Quartic Solution using VHM δ =0.5 δ =0.25 δ =0.125 δ =0.0625 Figure 11: Error in displacement for the quartic so-lution using the variable horizon method (VHM) forvarious values of horizon δ .15he extended domain method without any correc-tion factor. In this example, we propose to studythe effect of the analytical correction factor k ( x ) /κ ,as given in (18) and (19) and shown in Figure 12for δ = 0 . k ( x ) /κ = k ( x n − ) /κ = 8 / x and node x n − only. We consider the problem with the quartic so-lution as described in the previous section and plotin Figures 13 and 14 the errors in the displacementwhen using EDM I and EDM II for various valuesof the horizon δ . We still observe some spurious ef-fects in the solutions at x n − for both methods, asin EDM, but these tend to diminish as the horizongoes to zero. We compare in Figure 15 the absoluteerror in the displacement obtained with δ = 0 . δ = 0 .
125 and δ = 0 . . . . . . . x . . . . . . κ ( x ) / κ x x x n − x n − Correction Factor for EDM I κ ( x ) / κ , x ∈ (0 ,δ ) κ ( x ) / κ , x ∈ (1 − δ, Figure 12: Plot of the correction factor k ( x ) /κ pro-vided by (18) on the left-hand side and by (19) onthe right-hand side for δ = 0 . . . . . . . x . . . . . . . E rr o r i nd i s p l a ce m e n t Example with Quartic Solution using EDM I δ = 0 . δ = 0 . δ = 0 . δ = 0 . Figure 13: Error in displacement using EDM I withcorrection factor k ( x ) /κ given by (18) and (19) forvarious values of the horizon δ .16 . . . . . . x − . − . − . − . − . − . − . − . . E rr o r i nd i s p l a ce m e n t Example with Quartic Solution using EDM II δ = 0 . δ = 0 . δ = 0 . δ = 0 . Figure 14: Error in displacement using EDM II withcorrection factor k ( x ) /κ = k ( x n − ) /κ = 8 / δ . . . . . . . x . . . . . . . A b s o l u t ee rr o r i nd i s p l a ce m e n t Accuracy Study with Quartic Solution for δ = 0 . EDMEDM IEDM II
Figure 15: Comparison of the absolute error in dis-placement for EDM without correction, EDM I, andEDM II for δ = 0 . E n n δ EDM EDM I EDM II16 0.125 0 . . . . . . . . . . . . × − Table 4: Maximum relative error E n for the quar-tic solution obtained with EDM without correction,EDM I, and EDM II for various values of the horizon δ . m h Maximum relative error E n . . . E n for the quarticsolution obtained with the variable horizon method(VHM) for δ = 1 / h = δ/m , m = 2, 4, and 8. The objective of this numerical example is to showthat the approximate solution by the variable hori-zon method converges to the exact solution of theperidynamic model. We choose for the manufacturedsolution the quartic function of the previous exam-ple, for which the corresponding loading term in theperidynamic problem is given by: f b ( x ) = x + δ v ( x )12 . We consider here two values of the horizon, δ = 1 / δ = 1 /
8, and study the convergence in h , with h = δ/m and m = 2, 4, 8. The maximum relative er-rors E n and errors in the displacement field are shownin Table 5 and Figure 16 and in Table 6 and in Fig-ure 17, for δ = 1 / δ = 1 /
8, respectively. Weclearly observe that the solutions converge to zerowith respect to the discretization parameter h at arate of about two in both cases.17 . . . . . . x . . . . . . E rr o r i nd i s p l a ce m e n t Convergence Study with Quartic Solution using VHM m = 2 m = 4 m = 8 Figure 16: Error in displacement for the quartic so-lution using the variable horizon method (VHM) for δ = 1 / m , with h = δ/m . m h Maximum relative error E n . . . × − Table 6: Maximum relative error E n for the quarticsolution obtained with the variable horizon method(VHM) for δ = 1 / h = δ/m . . . . . . . x . . . . . . E rr o r i nd i s p l a ce m e n t Convergence Study with Quartic Solution using VHM m = 2 m = 4 m = 8 Figure 17: Error in displacement for the quartic so-lution using the variable horizon method (VHM) for δ = 1 / m , with h = δ/m . . . . . . . x . . . . . u ( x ) (cid:15) = 0 . (cid:15) = 0 . Figure 18: Manufactured solution u ( x ) given in (43)for (cid:15) = 0 . (cid:15) = 0 . We consider in this example the manufactured solu-tion to the classical linear elasticity problem: u ( x ) = x − (cid:0) e − (1 − x ) /(cid:15) − e − /(cid:15) (cid:1) / (cid:0) − e − /(cid:15) (cid:1) . (43)The corresponding load f b ( x ) was calculated usingthe computer algebra system Maxima and the bound-ary data g for the Neumann boundary condition isgiven by: g = EAu (cid:48) (1) = 1 − (cid:15) (1 − e − /(cid:15) ) . We observe that the magnitude of the first derivative u (cid:48) at x = 1 increases as 1 /(cid:15) , as illustrated in Figure 18in the cases where (cid:15) = 0 . (cid:15) = 0 . x = 1, we eval-uate in this example the maximum absolute error E n : E n = max i =1 ,...,n | u ( x i ) − u i | , rather than the maximum relative error as before.These errors are reported in Tables 7 and 8 for (cid:15) = 0 . (cid:15) = 0 .
01, respectively. We observe that a largenumber of degrees of freedom, even more so when (cid:15) decreases, is necessary to attain small errors. How-ever, it is noteworthy that the solutions still converge18aximum absolute error E n n δ LLEM VHM32 / . . / . . / . . / . . / . . E n for the casewith (cid:15) = 0 . E n n δ LLEM VHM512 / . . / . . / . . / . . / . . E n for the casewith (cid:15) = 0 .
01 using the variable horizon method(VHM).with a rate of two for VHM with respect to the hori-zon δ .We show the displacement obtained with VHM forvarious values of δ for (cid:15) = 0 . (cid:15) = 0 .
01 in Fig-ures 19 and 20, respectively. It appears from thesefigures that the maximum error occurs mostly at theend point x = 1 and that a large number of degreesof freedom needs to be employed in order to recoverthe boundary layer. It implies that the source of dis-cretization error at x = 1 severely pollutes the dis-placement field inside the domain. We have presented in this paper two methods, theso-called extended domain method and variable hori-zon method, to enforce boundary conditions withinthe non-local bond-based peridynamic model. Themethods were developed based on the requirement . . . . . . x − . − . . . . . D i s p l a ce m e n t u Convergence Study with (cid:15) = δ = / δ = / δ = / u ( x ) Figure 19: Displacement for (cid:15) = 0 . δ = 1 / , 1 / ,and 1 / . . . . . . . x − . − . − . . . . D i s p l a ce m e n t u Convergence Study with (cid:15) = δ = / δ = / δ = / u ( x ) Figure 20: Displacement for (cid:15) = 0 .
01 using the vari-able horizon method (VHM) with δ = 1 / , 1 / ,and 1 / .19hat the solutions to the corresponding problems becompatible, as the horizon goes to zero, with the so-lution to the boundary-value problem derived fromthe classical linear elasticity model. The formula-tion of both methods was first derived at the con-tinuous level, independently of the choice of the dis-cretization scheme, and later discretized, as an exam-ple, using the finite difference method. The perfor-mance of both methods was assessed on a simple one-dimensional linear elastic bar, fixed at one end andsubjected to a traction at the other end, for variousbody force densities. The solutions of these manufac-tured problems were polynomial functions of degreeone to four, which allowed us to study the advantagesand disadvantages of each method.The extended domain method (EDM) without cor-rection provided the correct solution only in the caseof the problem with a linear solution, for which thebody force density is zero everywhere in the domain.In all other cases, we observed spurious artefacts inthe displacement field near the boundaries. The rea-son for this behavior was explained by the fact thatthe extended domain method fails to correctly esti-mate the second derivative at the grid points clos-est to the boundaries, as it should in order to becompatible with the local model. We have seen thatit was possible to improve the solutions obtained byEDM when using an analytical correction factor ora numerical correction factor. On the one hand, theanalytical correction factor is derived at the contin-uous level and is valid for any type of discretizationschemes. On the other hand, the numerical correc-tion factor depends on the choice of the differenti-ation scheme. We observed on a numerical examplethat the numerical correction factor actually providedbetter results than the analytical correction factor. Inview of this ambiguity, it is difficult to appreciate thesuperiority of one or the other correction approach.The variable horizon method (VHM) delivered, forall the numerical examples, consistent solutions, ifnot the same, with those obtained from the local lin-ear elasticity model. In that sense, VHM should beconsidered a better approach than EDM when ap-plying boundary conditions in non-local models. Itis also worth noting here that VHM can be viewedas a coupling method between a local model and a non-local model (see e.g. [32]), where, in that partic-ular case, the coupling interface happens to be at theboundary of the domain.VHM has shown its potential, through the prelimi-nary one-dimensional numerical examples, for enforc-ing boundary conditions in the bond-based peridy-namic model. In particular, we have theoretically andnumerically shown that VHM is a method of ordertwo in the horizon δ , with respect to the linear elastic-ity model, as well as in the discretization parameter h , with respect to the finite difference method. How-ever, additional work should be carried out in orderto confirm its performance in more general situations.For example, the method should be extended to two-and three-dimensional problems and assessed in thecase of state-based peridynamic modeling. One couldalso consider other functions δ v ( x ) than the piecewiselinear function studied here and analyze their effectson the resulting solutions. Finally, the method shouldbe assessed when other discretization methods, e.g.the Finite Element method, are used for its imple-mentation. Supplementary materials
The Python code (using numpy [33, 34], scipy [35],and matplotlib [36]) used to generate the numericalresults is available on Github and on Zenodo [37]. Acknowledgements
Serge Prudhomme is grateful for the support of thiswork by a Discovery Grant from the Natural Sci-ences and Engineering Research Council of Canada(Award https://github.com/diehlpk/paperBCBBPD eferences [1] S. A. Silling, Reformulation of elasticity theoryfor discontinuities and long-range forces, Journalof the Mechanics and Physics of Solids 48 (2000)175–209.[2] P. Diehl, S. Prudhomme, M. L´evesque, A reviewof benchmark experiments for the validation ofperidynamics models, Journal of Peridynamicsand Nonlocal Modeling 1 (1) (2019) 14–35.[3] F. Bobaru, M. Yabg, L. F. Alves, S. A. Silling,E. Askari, J. Xu, Convergence, adaptive refine-ment, and scaling in 1D peridynamics, Interna-tional Journal for Numerical Methods in Engi-neering 77 (2009) 852–877.[4] G. Sarego, Q. V. Le, F. Bobaru, M. Zaccariotto,U. Galvanetto, Linearized state-based peridy-namics for 2-d problems, International Journalfor Numerical Methods in Engineering 108 (10)(2016) 1174–1197.[5] Q. V. Le, F. Bobaru, Surface corrections for peri-dynamic models in elasticity and fracture, Com-putational Mechanics 61 (4) (2018) 499–518.[6] S. A. Silling, E. Askari, A meshfree methodbased on the peridynamic model of solid me-chanics, Computers & Structures 83 (17-18)(2005) 1526–1535.[7] Q. Du, Nonlocal calculus of variations andwell-posedness of peridynamics, in: Hand-book of peridynamic modeling, Chapman andHall/CRC, 2016, pp. 101–124.[8] B. Aksoylu, F. Celiker, O. Kilicer, Nonlocal op-erators with local boundary conditions in higherdimensions, Advances in Computational Mathe-matics 45 (1) (2019) 453–492.[9] X. Gu, E. Madenci, Q. Zhang, Revisit of non-ordinary state-based peridynamics, EngineeringFracture Mechanics 190 (2018) 31–52. [10] E. Madenci, M. Dorduncu, A. Barut, N. Phan, Astate-based peridynamic analysis in a finite ele-ment framework, Engineering Fracture Mechan-ics 195 (2018) 104–128.[11] E. Madenci, M. Dorduncu, A. Barut, N. Phan,Weak form of peridynamics for nonlocal essen-tial and natural boundary conditions, ComputerMethods in Applied Mechanics and Engineering337 (2018) 598–631.[12] M. Gunzburger, R. B. Lehoucq, A nonlocal vec-tor calculus with application to nonlocal bound-ary value problems, Multiscale Modeling & Sim-ulation 8 (5) (2010) 1581–1598.[13] B. Aksoylu, T. Mengesha, Results on nonlocalboundary value problems, Numerical FunctionalAnalysis and Optimization 31 (12) (2010) 1301–1317.[14] B. Aksoylu, M. L. Parks, Variational theory anddomain decomposition for nonlocal problems,Applied Mathematics and Computation 217 (14)(2011) 6498–6515.[15] P. Seleson, M. Gunzburger, M. L. Parks, Inter-face problems in nonlocal diffusion and sharptransitions between local and nonlocal domains,Computer Methods in Applied Mechanics andEngineering 266 (2013) 185–204.[16] P. Seleson, D. J. Littlewood, Convergence stud-ies in meshfree peridynamic simulations, Com-puters & Mathematics with Applications 71 (11)(2016) 2432–2448.[17] W. Gerstle, N. Sau, S. Silling, Peridynamic mod-eling of plain and reinforced concrete structures,in: SMiRT18: 18th lnt. Conf. Struct. Mech. Re-act. Technol., Beijing, 2005, pp. 54–68.[18] E. Madenci, E. Oterkus, Peridynamic theoryand its applications, Vol. 17, Springer, 2014.[19] S. Oterkus, E. Madenci, A. Agwai, Peridy-namic thermal diffusion, Journal of Computa-tional Physics 265 (2014) 71–96.2120] K. Zhou, Q. Du, Mathematical and numericalanalysis of linear peridynamics models with non-local boundary conditions, SIAM J. NumericalAnalysis 48 (5) (2010) 1759–1780.[21] F. Bobaru, Y. D. Ha, Adaptive refinement andmultiscale modeling in 2D peridynamics, Jour-nal for Multiscale Computational Engineering9 (6) (2011) 635–659.[22] S. Silling, D. Littlewood, P. Seleson, Variablehorizon in a peridynamic medium, Journal ofMechanics of Materials and Structures 10 (5)(2015) 591–612.[23] S. Silling, D. Littlewood, P. Seleson, Variablehorizon in a peridynamic medium, Journal ofMechanics of Materials and Structures 10 (5)(2015) 591–612.[24] P. Seleson, Q. Du, M. L. Parks, On the con-sistency between nearest-neighbor peridynamicdiscretizations and discretized classical elasticitymodels, Computer Methods in Applied Mechan-ics and Engineering 311 (2016) 698–722.[25] M. L. Parks, R. B. Lehoucq, S. J. Plimpton,S. A. Silling, Implementing peridynamics withina molecular dynamics code, Computer PhysicsCommunications 179 (11) (2008) 777–783.[26] S. A. Silling, Introduction to peridynamics, in:Handbook of peridynamic modeling, Chapmanand Hall/CRC, 2016, pp. 63–98.[27] E. Oterkus, Peridynamic theory for modelingthree-dimensional damage growth in metallicand composite structures, Ph.D. thesis, TheUniversity of Arizona (2010).[28] R. W. Macek, S. A. Silling, Peridynamics viafinite element analysis, Finite Elements in Anal-ysis and Design 43 (2007) 1169–1178.[29] J. Mitchell, S. Silling, D. Littlewood, A position-aware linear solid constitutive model for peridy-namics, Journal of Mechanics of Materials andStructures 10 (5) (2015) 539–557. [30] P. Seleson, Improved one-point quadrature al-gorithms for two-dimensional peridynamic mod-els based on analytical calculations, ComputerMethods in Applied Mechanics and Engineering282 (2014) 184–217.[31] N. Trask, H. You, Y. Yu, M. L. Parks, Anasymptotically compatible meshfree quadraturerule for nonlocal problems with applications toperidynamics, Computer Methods in AppliedMechanics and Engineering 343 (2019) 151–165.[32] M. Zaccariotto, T. Mudric, D. Tomasi, A. Sho-jaei, U. Galvanetto, Coupling of FEM mesheswith peridynamic grids, Computer Methods inApplied Mechanics and Engineering 330 (2018)471–497.[33] T. E. Oliphant, A guide to NumPy, Vol. 1, Trel-gol Publishing USA, 2006.[34] S. van der Walt, S. C. Colbert, G. Varoquaux,The numpy array: A structure for efficient nu-merical computation, Computing in Science En-gineering 13 (2) (2011) 22–30.[35] P. Virtanen, R. Gommers, T. E. Oliphant,M. Haberland, T. Reddy, D. Cournapeau,E. Burovski, P. Peterson, W. Weckesser,J. Bright, S. J. van der Walt, M. Brett, J. Wil-son, K. Jarrod Millman, N. Mayorov, A. R. J.Nelson, E. Jones, R. Kern, E. Larson, C. Carey,˙I. Polat, Y. Feng, E. W. Moore, J. Vand er-Plas, D. Laxalde, J. Perktold, R. Cimrman,I. Henriksen, E. A. Quintero, C. R. Harris,A. M. Archibald, A. H. Ribeiro, F. Pedregosa,P. van Mulbregt, and SciPy 1.0 Contributors,SciPy 1.0: Fundamental Algorithms for Scien-tific Computing in Python, Nature Methods 17(2020) 261–272.[36] J. D. Hunter, Matplotlib: A 2d graphics environ-ment, Computing in Science & Engineering 9 (3)(2007) 90–95. doi:10.1109/MCSE.2007.55 .[37] P. Diehl, S. Prudhomme, Supplementary mate-rial: On the treatment of boundary conditionsfor bond-based peridynamics models (Jul.22020). doi:10.5281/zenodo.3942681 .URL https://doi.org/10.5281/zenodo.3942681https://doi.org/10.5281/zenodo.3942681