Improved wave dispersion properties in 1D and 2D bond-based peridynamic media

In this study, a novel method for improving the simulation of wave propagation in Peridynamic (PD) media is investigated. Initially, the dispersion properties of the nonlocal Bond-Based Peridynamic model are computed for 1-D and 2-D uniform grids. The optimization problem, developed through inverse analysis, is set up by comparing exact and numerical dispersion and minimizing the error. Various weighted residual techniques, i.e., point collocation, sub-domain collocation, least square approximation and the Galerkin method, are adopted and the modification of the wave dispersion is then proposed. It is found that the proposed methods are able to significantly improve the description of wave dispersion phenomena in both 1-D and 2-D PD models.


Introduction
The topic of the present work is the study of wave propagation in solid media described with a Peridynamic model. Propagation of waves is a classical problem of continuum mechanics [1,2] which is still the center of an active research effort concerning modern materials, applications and numerical methods [3][4][5][6][7][8][9][10]. The key dynamic properties of a medium supporting wave propagation are spectral characteristics. These comprise amplitude-frequency [11] and phase-frequency, or dispersion curves [12][13][14][15].
The theory of Peridynamics has been intensively developed over the last two decades. This theory was initially proposed by Silling [16] and Silling et al. [17] as an alternative nonlocal continuum theory able to deal with crack initiation and propagation. In contrast to the classical continuum mechanics which uses differential equations to describe the mechanical behavior of structures, PD applies an integrodifferential equation. The use of this type of equations is considered as an advantage for PD since it inherently enables B R. Alebrahim reza.alebrahim@unipd.it 1 modeling damage and crack propagation. Peridynamics has been found successful in modeling static and quasi-static problems and was recently also applied to simulate wave propagation, shock-waves and impact phenomena [18]. The initial investigation in this area was carried out by Silling [16] who found that the propagation of linear elastic waves with long wavelengths when simulated with PD or conventional theories is identical. In another study wave attenuation in a viscoelastic Peridynamics media is introduced by Silling [19]. The effect of non-locality in the suppression of traveling waves in a dissipative material is addressed. Detailed analysis of wave dispersion in nonlocal models in 1-D and 2-D cases was presented in [20]. It was clearly shown that Bond-Based (BB) Peridynamics displays better dynamic properties than the State-Based (SB) version, and both show dynamically softening behavior (i.e., predict lower wave speeds). A critical conclusion from [20] related to BB-and SB-PD is that due to its intrinsic integral properties, Peridynamics cannot yield a scheme more than second-order accurate. The comparison of wave dispersion in discrete and continuum Peridynamics models was addressed by Mutnuri and Gopalakrishnan [21]. They have introduced a method to modify wave dispersion in discretized Bond-Based PD model at low frequencies.
Wave propagation in coupled discrete models was investigated by some researchers in most cases proposing different approaches to implement the coupling [22][23][24][25][26][27][28]. Furthermore, Giannakeas et al. [29] studied wave dispersion in a cou-pled FEM-PD model, testing various strategies for coupling finite element method and Peridynamics. In another study Wildman and Gazonas [30] improved the wave dispersion behavior in Peridynamic through combining with a finite difference method. They have concluded that their proposed method can mitigate wave dispersion and it is able to reduce the oscillation of the stress intensity factor at the crack tip. Modeling of Lamb wave propagation in bimaterial plates was proposed by Alebrahim [31] who concluded that although PD produces accurate solution for the evaluation of symmetric Lamb waves, it is not able to model antisymmetric Lamb waves accurately.
In subsequent analyses, [13], the dispersion in SB-PD models was considered. One-, two-and three-dimensional dispersion properties of continuous and discretized PD equations were investigated. The influence of various model parameters, e.g., the horizon size, the number of nodes within the horizon and the shape of the influence function are analyzed in the context of dispersion properties. It was found that the numerical dispersion can be substantially improved by a proper choice of the influence function, in particular if it coincides with the kernel function of the nonlocal solid proposed in [32]. Such an approach was presented in [33] for computing weight coefficients correcting the numerical dispersion error. Weights distributions appear similar to multiple-point stencils in Finite Differences (FDs) and display negative coefficients that may lead to instabilities for bond-breaking schemes. In [34] dynamics of wave propagation in the split Hopkinson pressure bar with 1-D PD is analyzed. Similarly to other works, 1-D PD dispersion is investigated and analogously to many other authors, in [34] the original PD equation is first Fourier-transformed to the spectral (i.e., wavenumber-frequency) domain, and then discretized to yield the numerical dispersion formula. However, in fact, the time-space domain equation should be discretized first, followed by the spectral transform in the discrete domain. Note that the two approaches would yield the same results only under certain assumptions (see, e.g., [35] for a discussion on the time discretization influence). A proper approach can be found in [33], with a slight correction that the discrete version of the wavefield must be used to recover the dispersion relation. In [36] the mathematical equivalence of the SB-PD (and, therefore, BB-PD) with other mesh-free methods was proposed. In particular, differences between various mesh-free approaches and Peridynamics were shown in the context of influence function value that can be derived from alternative forms of the moving least squares technique. The influence of various micromodulus functions, also including negative weights, on wave propagation was provided in [37]. It was found that, for certain micromodulus functions, material instabilities due to initial perturbations with discontinuities may arise. The latter seems particularly important when wave propagation originating at crack-like discontinuity, e.g., due to a sudden energy release, is considered. According to the above-mentioned literature review, the modification of waves in 2-D PD has been implemented in a few studies. This study is mainly focused on providing a consistent and versatile framework for developing improved PD schemes. Various residual techniques are employed and the scaling coefficient corresponding to each PD bond is then computed. A set of scaling coefficients for each residual method is finally obtained to modify both pressure and shear wave dispersions in 2D PD. Employing the calculated scaling coefficients and based on the Fourier curve fitting method, a new influence function is introduced. Finally, using the influence function in the PD formulation, crack propagation in different samples is studied. This paper is organized in the following manner. Section 2 gives an overview of the peridynamic formulation. Section 3 provides numerical modeling of wave dispersion in 1-D and 2-D PD. Sections 4 and 5 introduce the modification of optimal 1-D and 2-D PD models, respectively. Section 6 investigates the minimization procedure. The proposed methods are validated against the analytical solutions and dynamic transient simulation in Sect. 7. Finally, the conclusion of the study is brought in Sect. 8.

The peridynamic formulation
In the nonlocal continuum theory of Peridynamics, the equation of motion of a material point in the 2-D BB-PD model in the presence of the body force is defined as follows (1) in Eq. (1), H x 0 represents the neighborhood of the material point X 0 . u 0 and u are vectors of displacement of the material points X 0 and X, respectively. ρ is density and b is body force density. The following definitions, according to Fig. 1, (1) is the pairwise force function and it is defined as Equation (2) can be rewritten in a linearized form as [16] f where s = |ξ |−|ξ | |ξ | in Eq. (2) is the stretch of the bond.ĉ(|ξ |) is a continuous micromodulus function and it is defined aŝ c(|ξ |) = cα(|ξ |), where for plane strain model c = 48E 5πδ 3 and Fig. 1 Schematic of a bond between node X 0 and X in a Peridynamic model before, ξ , and after deformation,ξ α(|ξ |) is an influence function to be computed (unknown). In a homogeneous linear elastic material under isotropic strain [17] x = y = = s, in plane strain conditions, the equivalence of the strain energy density in classical continuum mechanics and in 2-D BB-PD provides the following equation after canceling out the strain squared in BB-PD, s 2 , and classical continuum mechanics, 2 , and doing some simplifications, Eq. (4) can be rewritten in the following form with δ being the horizon and E is the Young's modulus. Since the influence function α(|ξ |) is unknown, we solve Eq. (1) numerically to have optimum wave dispersion phenomenon in BB-PD. In this way we are able to find the discretized form of the influence function.

Numerical dispersion in 1-D and 2-D peridynamics
In this section, we outline theory for generalized nonlocal models for wave propagation. Particular versions of these general models are then shown to yield the Bond-Based Peridynamics formulations for a particular set of model parameters. Iteration formulas for the model are considered in the displacement form and are next used to derive spectral characteristics of Peridynamics models for 1-D and 2-D cases. This section is concluded with the analysis of the non-

One-dimensional case: bond-based peridynamics
For analyzing dispersion in 1-D linear Peridynamics, the equation of motion of a material point in Eq. (1) in the absence of body force is reduced to where f =ĉ(ξ )s is the pairwise force function andĉ(ξ ) = cα(ξ ) where c = 2E δ 2 A and α(ξ ) is the influence function. A is the area of cross section. Equation (6) is written in the discretized form as follows where N , according to Fig. 2, is the number of nodes within the horizon δ (i.e., δ = N , with being the node spacing). Noting that in this study the influence function α(|ξ |) in the discretized form α(|ξ i |) is called stiffness scaling coefficient. After discretization in time, Eq. (7) can be used to compute displacements in subsequent time steps for a 1-D infinite domain. Wave propagation characteristics are found from Eq. (7) assuming a single, unit amplitude, wave solution for a linear, homogeneous, isotropic medium of the form u n (n , t) = e ikn e −iωt (thus continuous time-domain). This approximation is valid for small time steps t, which is typically the case under explicit time integration scheme assumptions. Further, i = √ −1 is the imaginary number and k is the wavenumber. For a regular grid of nodes, the following phase relation holds between any node in the stencil and the considered node n u n±i = u n e ±iki . (8) With (8) and after rearrangement, Eq. (7) becomes with the non-trivial solution for the quantity in brackets {·} = 0, yielding the dispersion relation for a PD discretized model

Two-dimensional case: bond-based peridynamics
A framework for obtaining dispersion relation in 2-D PD is introduced in this section. In contrast to the 1-D case, two wave modes-longitudinal and shear-need to be considered in a 2-D solid medium. Following, dispersion relations for a 2-D infinite PD discretized model are derived. Calling Eq. (1) and rewriting the 2-D BB-PD formulation in the discretized form in the absence of body force, one will obtain where u 0 = [u x0 , u y0 ] T and u i = [u xi , u yi ] T are displacement vectors of the considered (0) and a neighboring (ith) node from the reference to the deformed configuration, respectively. The same grid spacing in x and y directions is considered, x = y = .N is the total number of bonds in the horizon. Assuming a single, unit amplitude, wave solution for a linear, homogeneous and isotropic medium, the following phase relation holds between the grid points where k = [k x , k y ] (with k x = |k| cos β, k y = |k| sin β) is the wavevector, ξ i = [ξ xi , ξ yi ] T and ξ xi , ξ yi denote distances of the ith (i = 0) node from the central node (i = 0) in x and y directions, respectively. Assuming time dependence of the form u 0 = p 0 e −iωt and substituting Eq. 12 in Eqs. 11 and after arrangement, one will find ⎛ where p 0 = [p x0 , p y0 ] T is the wave polarization vector. For non-trivial solutions of Eq. 13, the determinant of the expression in brackets must vanish, yielding the characteristic equation the determinant from Eq. 14 yields the following solution for ω 2 (see A for details) where θ i in the linearized form of the 2-D BB-PD formulation [16] is the angle of ith bond in the reference configuration. cos(2θ i ) and sin(2θ i ) in Eq. (15) are equal to yi |ξ i | 2 and 2 ξ xi ξ yi |ξ i | 2 , respectively. Noting that Eq. (15) is valid for grids with translational symmetry, the even and odd properties of cos and sin functions allow to reduce the number of param-etersĀ i from the number of bonds in the four quadrants Q 1 − Q 4 to the number of bonds in the first quadrant Q 1 (see Fig. 3a, b). For spatially symmetric response, we require thē A i parameters to be equal for nodes located symmetrically with respect to x and y directions, i.e., for the four symmetric nodes at ξ j = [±ξ j x , ±ξ j y ] T . Then in Eqs. (16)- (18), j iterates over the four symmetric bonds for the ith bond (i = 0). Due to symmetry considerations, nodes located at x-or y-axis assume γ i = 2 (red nodes in Fig. 3a, b), while other nodes assume γ i = 4 (green nodes in Fig. 3a, b).

1-D BB peridynamic
A general form of the 1-D BB Peridynamic solid in the absence of the stiffness scaling coefficients for the long wavelength limit, i.e., k → 0 and kδ 1, can be derived from Eq. (10) The expansion in Eq. (19) is valid for kδ 1 and, consequently, for k 1 since δ ≥ , by definition. Consequently, Eq. (19) is valid for horizon radii much smaller than the wavelength. Noting that E/ρ = V 2 , where, V is phase speed. Equation (19) provides the dispersion relation of a BB-PD model solely as a function of the number of nodes within the horizon, N , The following observations can be made based on Eq. (20) -it is assumed that the product of wavenumber and horizon is small (kδ 1), -the dispersion depends on the number of nodes within the horizon, N , and converges to the exact solution as N → ∞, -the number of nodes within the horizon is related to the grid size by N = δ, indicating convergence for N → ∞ and 1/ → ∞, -the convergence of a BB-PD model for the long wavelength limit can be uniquely defined by a single parameter.

2-D BB peridynamic
In order to estimate the static error for 2-D BB-PD model in the absence of stiffness scaling coefficients, we evaluate Eq.
(15) at the long wavelength limit, i.e., for |k k k| → 0 and |k k k|δ 1, and considering Eqs. (16)-(18) at the long wavelength limit, we have c 2ρ c 2ρ withS ± i given as where i x = ξ i x and i y = ξ iy . Then, the long-wavelength approximation to the dispersion relation yields (note the summation over the first quadrant only) Please note that despite the requirement of |k k k| → 0, propagation angle β still influences Eq. (25) throughS ± i , i.e., the static error is angle-dependent due to the mesh-induced anisotropy in 2-D. Figure 4a, b show direction-dependent normalized wave velocity (i.e.,V s, p = ω s, p /|k| /V s, p ) profiles for N = {2, 3, 4} for the long wavelength limit. Recalling that the shear and pressure waves V s, p are related to V and ν by the scaling will depend on selection of ν. In the following we adopt ν = 0.25, which yields V s = 0.632V and V p = 1.095V . From Fig. 4a, b it can be clearly seen that for N = 2 wave propagation is strongly anisotropic, while for N ≥ 3 almost isotropic-propagation is observed. Figure 4c shows normalized wave velocities for two angles that are relevant for subsequent error analysis, i.e., β = {0, π/4}, for variable number N (up to N = 10). The values of shear and pressure wave velocities converge quickly towards the analytical solution (assuming ν = 0.25), i.e.,V k→0 s = 1 and V k→0 p = 1, respectively. Again, substantial differences are observed between N = 2 and N ≥ 3 for both propagation angles, indicating that N = 3 may be an optimal compromise between accuracy and computational time. It may be noted that when considering quasi-static crack propagation problems with Peridynamics (i.e., not wave propagation), higher number of nodes (e.g., N = 4 or N = 6) are frequently preferred [38].

1-D BB peridynamic
For the special case of limiting interactions to the closest neighbors only, i.e., N = 1, the dispersion relation yields Noting that for small values of k << 1, , over 40% error in the wave propagation speed is obtained for the BB-PD solid. Please note that no correction factors are considered in Eq. (27). The frequently employed volume correction method substantially improves the solution and reduces the error to zero. The result of Eq. (27) can be compared to the dispersion relation for a 1-D central-difference scheme, ω 2 = 2E/(ρ 2 )(1 − cos(k )), that reproduces the exact solution for the long wavelength limit (k << 1), ω ≈ √ E/ρk.

2-D BB peridynamic
Although not of practical importance, it is interesting to see the dispersion properties for closest-neighbors interactions in 2-D, namely N = 1. This compares directly to the 2-2 Finite Difference (FD) local approximation. From Eq. (25) for N = 1, the dispersion equation for small values of << 1 reduces to It is clear from Eq. (28) that the two wave modes propagate strongly anisotropically in the domain, with the two preferential directions, along the β = 0 and β = π/2 directions, at which only a single wave mode exists and propagate with an error of 60%. This heavily distorted wave propagation pattern is the result of the very limited stencil for the model (only five nodes), for which the diagonal nodes are missing which results in poor representation of the shear modes. For comparison with classical local models based on finite differences or cellular automata and the resulting dispersion the reader is referred to [39].
methods in terms of wave propagation. Therefore, in this section we propose a set of modified Peridynamic models with the main goal of improving PD dispersion properties. The proposed modifications are particularly attractive for regions where the product of the wavenumber and horizon (kδ) is not small. According to Eq. (7), the set of n (n = N ) scaling coefficients, α i , that best approximate dynamic properties of the model, can be found [40]. One possible way of finding α i , outlined in [33], was based on selecting a number of points in the (ω, k) space and enforcing correct behavior of the model at those points. The results, however, depend on the selection of frequency (or wavenumber points) and various scaling coefficients are obtained for different choices. The approach proposed here is based on constructing optimal approximation to the exact dispersion characteristic. Equation (10) is rewritten here Note that elements ofk ∝ k 2 . Equation 29 can be seen as a linear combination of basis functionsk, thus the coefficients α α α can be found by using various known techniques, e.g., Finite Elements, Weighted Residuals and others. Recalling the exact solution to dispersion of a 1-D medium, ω 2 = V 2 k 2 , Eq. (29) can be used to formulate the square-frequency residue The residue R in Eq. (31) can be used for computing α α α to achieve smallest possible R over the considered wavenumber domain. Following the standard approach for minimizing the residual we require where the integration is performed over the considered wavenumber domain and the weight w depends on the method adopted. Because of the discrete nature of the system, the integration will be performed over the irreducible Brillouin zone.

Physical restrictions on dispersion
The problem formulated in Eq. (32) consists of n = N unknowns α α α. Therefore, depending on the selected method, N equations must be constructed to uniquely solve the problem. However, Eq. (32) should satisfy two physically motivated conditions, i.e., zero static error (i.e., for ω = k = 0) and correct group velocity at ω → 0 and k → 0. Thus, additional constraints may be imposed on parameters α i to enforce correct solution at the long wavelength limit. Namely, we require that (a) at k → 0 we have ω → 0, and Recalling Eq. (19) for dispersion at the long wavelength limit, yields Clearly, it follows from Eq. (33) that the condition (a) is satisfied by definition. The group velocity condition at k → 0 then requires that with i = [1, 2, . . . , N ] T . Equation (34) constitutes an additional equation to the set defined in Eq. (32). Note, that the number of equations generated from Eq. (32), for the system to have unique solution, is now N − 1.

Modified optimal 2-D BB-PD model
In order to investigate the influence of the stiffness scaling coefficients, α i , on the 2-D BB-PD formulation, Eq. (15) is rearranged in the following form where the set of n =N scaling coefficients, α i , can be found through analogous methods as those used for the 1-D case and-analogously to Eq. (30)-|k| p,s ∝ |k| 2 . For a 2-D case, numerical dispersion depends on material properties of the medium and in-plane discretization parameters, see Eq. (35). Here, in contrast to the 1-D case, dispersion is direction dependent due to the discrete distribution of bonds in the Peridynamic grid that introduces a sort of anisotropy [38]. Noting that the scaling coefficients are assumed to follow the bonds distribution symmetry (as discussed above), we formulate the following residuum R e p , e s , ω 2 p , ω 2 where e p,s are scaling coefficients for the shear and longitudinal waves, respectively. The coefficients can be used to tune the accuracy of the solution, e.g., with e s = e p = 1, we treat both wave modes equally accurate. The residuum from Eq. (36) depends on the wavenumber magnitude and wave propagation angle. Error minimization requires that Due to symmetry considerations and unique characterization it is sufficient to use integration limits corresponding to the irreducible Brillouin zone in Eq. (37). Note also that due to symmetry, the number of unknown coefficients decreases significantly. For example, for N = 3 the number of independent coefficients according to Fig. 3b is 7. However, it will be explained later that due to symmetry constraints, the total number of independent coefficients for N = 3 reduces to 6.
Equation (37) generates n − 2 equations (n is the number of unique scaling coefficients for the bonds) whose form depend on the type of the method employed. The two missing equations follow from the long wavelength limit equations as presented next.

Physical restrictions on dispersion
Rewriting Eq. (35) at the long wavelength limit we require verifying two physical restrictions on the dispersion relations for the pressure and shear waves, i.e., we check (a) if ω s, p → 0 when |k| → 0 and (b) dω s, p /d|k| → V s, p g for |k| → 0. From (35) and using (26) one obtains Analogously to the 1-D case, for a 2-D BB-PD model the first condition, (a), is satisfied by definition. From (b) we get Equation (39) provides two constraints for α i parameters, that must be satisfied from physical viewpoint.

Minimization procedure
Residuals defined by Eqs. (31) and (36) for 1-D and 2-D, respectively, are weighted and minimized as given by Eqs. (32) and (37) in order to yield the sets of scaling coefficients α i . The weights, w, in Eqs. (32) and (37)  It needs to be emphasized that the proposed methodology for finding optimal scaling coefficients for the bonds, α, is based on the frequency-wavenumber domain analysis. The procedure relies not only on matching the spectral response of the model with the analytical solution, but also applies physical constrains in order to force proper low-frequency behavior (both value and slope). Consequently, Eqs. (34) and (39) for 1-D and 2-D, respectively, are used to reduce the number of residual equations required for unique solution. In practice, due to symmetry (point symmetry in 1-D and two symmetry planes in 2-D), we seek the total number of n unknowns in vector α, which is half of the total number of bonds for the considered node for 1-D, and 1/4th of all bonds for the considered node for 2-D. The constraint equations, Eqs. (34) and (39), reduce the number of residual equations to n − 1 for 1-D and to n − 2 for 2-D.
The integration in Eqs. (32) and (37) is carried out over the contour of the Irreducible Brillouin Zones (IBZ) characteristic (see Fig. 3c) for 1-D and 2-D cases, respectively, i.e., k ∈ [0, π/ ] for 1-D and Other approaches than those listed above are possible to compute coefficients α α α; however, in the following we focus on these four methods, as they are most widely used in engineering practice.

Point collocation
The weighted residual point collocation is essentially similar to the approach outlined in [33]. In the presented method, however, we implicitly link it to the weighted residuals techniques and supplement the residual by the physical constraints, thus reducing the number of equations. The collocation points distribution can be arbitrary; however, additional effort can be made in order to develop more rigorous schemes to improve the solution. In order to define the n − 1 equations via the point collocation method, the interval [0 π ] for 1-D and the IBZ contour for 2-D were uniformly divided into n − 2 divisions. For the 1-D case, from (31) and (32) we obtain while for 2-D, from Eqs. (36) and (37)

Sub-domain collocation
In this method, the interval [0 π ] for 1-D and the IBZ contour for 2-D, are uniformly split into n −1 and n −2 sub-domains, respectively. Residual equations are obtained through applying the integration in every sub-domain. For the 1-D case, from (31) and (32) while for 2-D, from Eqs. (36) and (37)

Least square approximation
With w = R, and making use of the physical conditions (34) and (39), for the 1-D case, from (31) and (32) we get while for 2-D, from Eqs. (36) and (37) The minimization procedure outlined in Eqs. (44) and (45) would lead to the minimum average squared-frequency errors over the wavenumber domain. The sets of Eqs. (44) and (45) define n − 1 and n − 2, respectively, (least-square) optimal scaling coefficients α i for the modified BB-PD approaches.

Galerkin method
The Galerkin method is applied to find the unknown scaling coefficients α α α corresponding to the optimized 1-D and 2-D non-local Peridynamic models. The set of unknown coefficients is obtained by projecting the residual, Eqs.

Rdk
where R is the residual from Eq. 31 andk is a n − 1 dimensional vector, as defined in Sect. 4. The n scaling coefficients α i can be found from Eqs. (34) and (46). It needs to be pointed out that the vector {Ā i } contains basis functions for all 4n bonds within the horizon. The scaling vector, α, on the other hand, contains only n − 2 unique coefficients. Consequently, for a 2-D system the following set of n − 2 equations is used to determine α where the summation is carried over the four i indices corresponding to the symmetric bonds (i.e., with the same α i ).

1-D wave dispersion
In this subsection, we investigate the influence of the optimal scaling coefficients on the modification of wave propagation in a 1-D domain. To validate the proposed method, we compare dispersion curves for the pressure wave obtained analytically in Sects. (3.1) and (4), with dispersion curves identified from a time-domain simulation of wave propagation in a bar. In order to obtain the dispersion curves with the time-domain simulation, we adopt the spatio-temporal Fourier transform procedure (see, e.g., [39]). We start by analyzing original (unmodified) BB-PD models in 1-D for N ∈ {1, 2, 3} and then we present numerical studies for optimal BB-PD models for N ∈ {1, 2, ≥ 3} along with their validation examples. Figure 5 shows results of timedomain simulations (gray-scale colormaps) compared with analytical predictions (black circles) for 1-D models. Figure  5a-c present results obtained for the original BB-PD models (i.e., with unmodified stiffness coefficients). Perfect agreement between numerical and analytical results can be clearly observed.

Modified BB-PD scheme for N = 1
Assuming N = 1, i.e., only closest neighbors are included in the horizon and therefore there is a single scaling coefficient parameter α 1 . Following Eq. (34) we have α 1 = 1/2. This result reproduces the coefficient correction in the original 1-D BB-PD after applying volume correction method. Using the dispersion relation of Eq. (33), ω 2 = V 2 k 2 2α 1 = V 2 k 2 , the result clearly reproduces the 2-2 FD scheme and perfectly agrees with the group velocity at the long wavelength limit (see Sect. 3.4). Note that this result holds regardless of the spatial spacing . Figure 5d shows results obtained from the numerical experiment and analytical results for this new model. The comparison of Fig. 5a, d shows that the cutoff frequency for the unmodified case is substantially higher than for the modified one, i.e., waves of higher frequencies are supported by the improved model. The proposed modified BB-PD model perfectly reproduces wave velocities at the long wavelength limit (note that the wave speed error for the unmodified model is over 40%), therefore displays much smaller wave speed errors at low frequencies than the unmodified case.

Modified BB-PD scheme for N = 2
Including second neighbors, i.e., N = 2, yields two scaling coefficient parameters α 1 and α 2 . From (34) we have the The second equation can be found using one of the methods outlined in Sect. 6. For instance employing the Galerkin method, one obtains α 1 = 3.72 and α 2 = −0.86. Inspection of Fig. 5b, e clearly shows the superiority of the proposed modified scheme. The cut-off frequency of the new modified model is nearly doubled when compared to the original BB-PD scheme. Also, the effective frequency band for which the new model can be applied, up to approximately 2.25rad/s, is much broader than for the reference method (up to approximately 1.25rad/s).

Modified BB-PD scheme for N ≥ 3
In contrast to the cases of N = 1 and N = 2, for which only the restriction equation [Eq. (34)] or the restriction equation and only one residual equation were used, respectively -we investigate the case of N ≥ 3 where we employ the whole set of residual equations to obtain optimal scaling coefficients. Scaling coefficients of the PD bonds for N ∈ {3, 5, 10} are presented in Table 1 for various minimization methods. It can be observed that there is a slight difference between the coefficients for the four methods investigated. Figure 5c, f show dispersion for N = 3 for the unmodified and modified models, respectively. Further deterioration of spectral properties can be observed for the original BB-PD, while substantial improvement is clearly seen for the modified model. Figure 6a collects scaling coefficients for higher values of N ∈ {3, 5, 10} obtained using the point collocation method. It can be noted that scaling coefficients assume values of monotonically decreasing magnitude and alternating signs. This last finding is consistent with the nonlocal media theory discussed in [32] and the analysis outlined in [20]. Negative scaling coefficients are critical for substantial improvement of spectral characteristics of the new optimal methods; however negative stiffness values are expected to cause instabilities when bond breaking is allowed in the model (e.g., in crack propagation simulations [13]). The error of the dispersion properties as a function of the horizon size was obtained and is plotted in Fig. 6b. The following formulation is employed to obtain error (square-average) values for the analyzed models It can be observed that there is a considerable difference between the original and modified PD for increasing radius of horizon. The error of the modified PD vanishes rapidly as the size of the horizon increases. According to Fig. 6b, the error of the proposed methods for N ≥ 8 drops drastically. The original and modified dispersion curves for N ∈ {3, 5} are shown in Fig. 7a, b. It can be seen that the accuracy of the models, in terms of their spectral response, has significantly improved. Clearly, all models display relatively accurate response at low frequencies, but substantial differences are noted at moderate and high frequencies. Comparing the N = 3 and N = 5 cases for unmodified models, it can be observed that the effective frequency band decreases despite the horizon is extended. Exactly the opposite effect can be seen for the modified models, where the number of scaling coefficients is greater for N = 5 than for N = 3; therefore that model can be better tuned to fit the analytical response. A comparison of errors for the three analyzed horizons, i.e., N ∈ {3, 5, 10}, and for the four modification approaches are shown in Fig. 7c. These errors are computed according to Eq. (49) and are normalized using the original BB-PD error for N = 3 (see Fig. 6b). It can be observed that among the four introduced methods, the point collocation and the Galerkin method produce the maximum and the minimum error, respectively, for all cases considered.

2-D wave dispersion
As mentioned earlier and according to the reference [33], it is clear that applying negative scaling coefficients in the PD formulation can produce instability in the PD solution upon damage initiation. Therefore, in this study two different solutions for the calculation of the scaling coefficients in 2-D BB-PD are introduced. In the first solution (hereafter it is called optimal solution), positive and negative scaling coefficients are obtained. Note that in this method no constraint is imposed on the solution procedure. However, in the second method, positive-only coefficients are obtained through introducing lower bound constraints in the solution procedure.

Optimal solution
As mentioned earlier, the optimal solution is employed to modify the PD formulation through introducing positivenegative scaling coefficients. Due to negative coefficients, this method is more suitable for the improvement of wave propagation. The scaling coefficients corresponding to the 2-D BB-PD for N = 2 and N = 3 are tabulated in Table 2.
The numerically and analytically computed wave dispersion for the unmodified 2-D BB-PD models with N = 2 and N = 3 is plotted in Fig. 8a, b, e-f, respectively. The dispersion of pressure and shear waves computed analytically according to Eq. (15) is denoted by markers and it is compared with dynamic transient simulation results (gray colormaps). Perfect fit between the analytical and numerical responses can be seen for both pressure and shear waves. As shown in Fig. 8a, b there is a cutoff frequency for the pressure and shear dispersion curves. By increasing the size of the horizon in Fig. 8e, f, the cutoff frequency shifts towards lower frequencies for the original BB-PD models. In Fig.  8a (but also for other results in 2-D) ghost higher order modes (HOM) are also present in the dynamic simulation. These higher order modes should be in general avoided but their detailed analysis and mitigation methods are beyond the scope of this paper. As it is brought in Fig. 8c, d, g, h for N = 2 and N = 3, respectively, dispersion curves in the 2-D PD were modified through defining scaling coefficients using Eq. (35). Figure 8c, d reflect the modified pressure and shear wave dispersion for N = 2 and in the same way Fig. 8g, h show the wave dispersion for that of N = 3. All unknown  Table 2. Number and order of the scaling coefficients for N ∈ {2, 3} are shown in Fig. 3a, b. According to Fig. 3a, b, due to symmetry, the number of coefficients in Table 2 is one-fourth of the whole bonds of a node. Moreover, because of the symmetry with respect to the diagonal axis, the number of independent coefficients for N > 3 is less than one-fourth of the whole bonds of a node (see Fig. 3 and Table 2). The wave dispersion plotted in Fig. 9a, b clearly demonstrate that the scaling coefficients have effectively modified the results of the original 2-D PD using optimal solution. It is clearly shown that both pressure and shear waves have effectively improved using a single set of scaling coefficients.

Positive-only solution
The scaling coefficients for N = 3 are calculated by introducing the positive-only constraint (see Table 2). In this method, MATLAB nonlinear least square solver function is employed to compute the positive scaling coefficients. In order to find proper scaling coefficients for the improvement of both crack and wave propagation, we have imposed a constraint on the lower bound of the coefficients. The lower bound value 0.5 is chosen in this study. More details regarding the positiveonly method and the lower bound value can be found in [41]. These scaling coefficients are computed using point collocation, sub-domain collocation, least square approximation and the Galerkin method. Employing positive-only solution, the results of the applied residual techniques are very close to one  another. As it is shown in Fig. 9c, d, using above positiveonly scaling coefficients in the PD formulation effectively mitigates pressure and shear dispersion relations. Based on the calculated positive-only scaling coefficients, one is able to use curve fitting method to find a new influence function. After employing the Fourier curve fitting, the influence function is computed as follows where a 0 = 2.3, a 1 = 2.0, a 2 = −0.39, a 3 = −0.48, a 4 = 0.28, a 5 = 0.14, a 6 = −0.15, a 7 = −0.02, a 8 = 0.05 and w = 1.1. Note that above influence function is computed in such a way that it is firmly compatible with strain energy equivalence mentioned in Eq. (5). The curve corresponding to the influence function is shown in Fig. 10. The red circles represent the scaling coefficients calculated through

Error calculation
The angular frequencies and the error corresponding to each method is calculated and compared in Fig. 11. Angular frequency of the wave over the IBZ contour using optimal and positive-only solutions for various minimization techniques are compared against analytical method in Fig. 11a, b, respectively. Using optimal solution, the modified curves in Fig. 11a follow the local result (black solid line) much more accurate than those curves in Fig. 11b which are obtained through positive-only solution. Errors corresponding to each method were calculated using Eq. (51).
(51) Figure 11c shows the error corresponding to each method. Error values are normalized based on the error related to the original BB-PD method. According to the information in Fig.  11c, the error of the optimal solution is much more lower than that of positive-only one. In the optimal and positive-only methods, point collocation and least square approximation are showing the minimum error, respectively.

Crack propagation
Crack propagation in the pre-notched plates with different loading conditions (see Figs. 12 and 13) is investigated using the modified influence function in Eq. (50). In order to simulate crack propagation, initially the critical stretch, s 0 , needs to be determined. According to the reference [17], the fracture energy G 0 in the BB-PD model can be computed as follows employing Eq. (50) in Eq. (52), the critical stretch in the BB-PD model for the modified influence function can be evaluated. Considering the Soda-lime glass material properties in the first example (see Fig. 12), E = 72 GPa, ρ = 2440 kg/m 3 , ν = 0.22 and G 0 = 135 J/m 2 , the critical stretch in the plane strain condition for the original and modified BB-PD is s 0 = 0.0013 and s 0 = 0.0014, respectively. In the pre-notched plate in Fig. 12a, the load σ = 12 MPa is applied uniformly at top and bottom surfaces. Horizon value and time step in this study are δ = 0.75 mm (N = 3) and t = 25 ns, respectively. Figure 12b, c demonstrate crack propagation at t = 46 µs using the modified and original BB-PD formulation, respectively. The variation of the crack branching position in the modified model is small in comparison with that of original one. This achievement is in good agreement with that which is addressed in [42] for constant and conical micro-modulus functions.
In the second crack modeling study, a square plate with an inclined crack 45 • under uniform stretch loading σ x = σ y = 23 MPa in x and y directions is considered (see Fig.  13a). Length of the plate and crack are 0.05 m and 0.01 m, respectively. The square plate is made of Duran 50 glass with following material properties, E = 65 GPa, ρ = 2235 kg/m 3 , ν = 0.2 and G 0 = 204 J/m 2 . Using Eqs. (52) and (50), the critical stretch s 0 in the original and modified BB-PD formulation is 0.0023 and 0.0026, respectively. Crack propagation in the square plate at t = 22 µs using modified and original PD formulation is shown in Fig. 13b, c. It can be observed that, similar to the original BB-PD model, branching in the modified method starts from the pre-crack tip. Moreover, the effect of the introduced influence function on the crack path is negligible. These achievements are in agreement with the results of reference [43]. Finally, based on the above crack modeling studies, it is observed that the modified BB-PD formulation can be adopted to accurately simulate crack propagation problems. Furthermore, this approach introduces a more accurate wave dispersion relation in comparison to that of original BB-PD.

Conclusion
Wave dispersion in 1-D and 2-D Peridynamic media was investigated. Deficiencies of the BB-PD formulation in the modeling of wave propagation were studied. Four different residual techniques, i.e., point collocation, sub-domain collo-cation, least square approximation and the Galerkin method, were proposed and scaling coefficients corresponding to all bonds were computed. The main achievements of the study are pointed out as follows (1) A set of scaling coefficients to be applied to the stiffness of the bonds to improve wave dispersion in 1-D PD was introduced and it was observed that, among four different minimization approaches, the Galerkin method produces more accurate wave dispersion properties. (2) A set of scaling coefficients to improve pressure and shear waves in 2-D PD was defined using optimal and positiveonly solutions. It was found that the optimal solution provides more accurate wave dispersion relation. However, it doesn't produce stable crack propagation. (3) Performing the integration over the irreducible Brillouin zone, different minimization approaches were implemented and it was seen that in the optimal and positive-only solutions, point collocation and least square approximation produce more accurate wave dispersion relations, respectively. (4) Using the positive-only scaling coefficients, a new influence function was introduced through Fourier curve fitting method. The proposed influence function improves the wave dispersion relation in the BB-PD formulation.

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

A Characteristic equation
The characteristic equation for a 2-D non-local model is obtained from Equation 53 can be rewritten in a simpler form as follows the characteristic equation assumes bi-quadratic form and yields ρ 2 ω 4 − ρ(C +S)ω 2 +CS −P 2 = 0.