Dirichlet-type absorbing boundary conditions for peridynamic scalar waves in two-dimensional viscous media

Construction of absorbing boundary conditions (ABCs) for nonlocal models is generally challenging, primarily due to the fact that nonlocal operators are commonly associated with volume constrained boundary conditions. Moreover, application of Fourier and Laplace transforms, which are essential for the major-ity of available methods for ABCs, to nonlocal models is complicated. In this paper, we propose a simple method to construct accurate ABCs for peridynamic scalar wave-type problems in viscous media. The proposed ABCs are constructed in the time and space domains and are of Dirichlet type. Consequently, their implementation is relatively simple, since no derivatives of the wave field are required. The proposed ABCs are derived at the continuum level, from a semi-analytical solution of the exterior domain using harmonic exponential basis functions in space and time (plane-wave modes). The numerical implementation is done using a meshfree collocation approach employed within a boundary layer adjacent to the interior domain boundary. The modes satisfy the peridynamic numerical dispersion relation, resulting in a compatible solution of

(ABCs) are the most prominent ABMs, and development of new ABMs is still an ongoing topic of research. [5][6][7][8] Essentially, the idea of ABMs is to limit the computational domain to a region of interest called the near-field, while referring to the area outside that region as the far-field. Provided that appropriate (artificial) boundary conditions are imposed on the boundary of the near-field (referred to as the truncating boundary), the solution of an unbounded problem can be obtained within this region. In the context of scalar waves, which is the application focus of the present study, this implies that any wave impinging on the truncating boundary is transmitted out of the near-field without any reflection.
Originally developed for simulations of electromagnetic waves, 9 the PML method provides a means for describing wave propagation in unbounded media through the imposition of an absorbing layer. In theory, an absorbing layer can be constructed such that any incident wave is exponentially attenuated. Moreover, the PML approach is applicable within a finite element framework for the implementation of unbounded domain problems. 10 While preliminary success with numerical evidence of effectiveness has been demonstrated with the PML method applied to the nonlocal Helmholtz equation, space fractional partial differential equations (PDEs) and peridynamics in two dimensions (2D), [11][12][13] the method itself suffers from several drawbacks. Most importantly, implementation of the PML method results in grid-dependent performance and numerical instability of discrete systems. 14,15 ABCs have been the subject of decade-long investigations since their first introduction in Reference 16 and their further development into high-order ABCs. 4,17 An important advantage of ABCs is that they are directly constructed in the time and space domains (in contrast to PMLs) and, in the case of high-order ABCs, up to any desired order. 4 In the context of wave propagation problems, most ABCs have been developed in the framework of the classical (local) theories. 18 The extension to peridynamic (PD) wave-type equations is challenging as the interactions between material points are nonlocal and the corresponding nonlocal operators are associated with volume constrained boundary conditions. [19][20][21] For this case, the calculation of the kernel function from the Laplace transform becomes complicated due to the nonlocality. 22 In addition, for high-order ABCs, auxiliary variables on the truncating boundary are required to cope with the calculation of higher-order spatial derivatives. For nonlocal models, this complicates the development of appropriate high-order ABCs, possibly rendering their implementation exceedingly difficult, since the accurate imposition of Neumann-and Dirichlet-type boundary conditions for PD models is still an ongoing subject of research. 23,24 More recently, a number of studies have been conducted on the application of ABCs to nonlocal models, and examples of those studies are as follows. In one dimension (1D), ABCs are considered for the PD scalar wave-type equation in Reference 25 and for the nonlocal Schrödinger equation in References 26,27. For the 1D PD equation, ABCs based on recursively related kernel functions are constructed in Reference 28, and transmitting boundary conditions with matching dispersion relations are developed in Reference 29. The nonlocal diffusion equation in unbounded domains along with the construction of ABCs in 1D is treated in References 22,30. Inspired by first-kind integral equation methods, the authors in Reference 31 developed ABCs for PD scalar wave-type equations as a generalized Dirichlet-to-Dirichlet type boundary condition in 2D. The stability and convergence analysis of high-order numerical approximations for nonlocal wave equations in unbounded domains using quadrature-based finite difference schemes and Dirichlet-to-Neumann-type ABCs in 2D is presented in Reference 32. For 2D PD models, accurate ABCs are developed in Reference 33, and ABCs based on an attenuating layer with an increased damping are proposed in Reference 34. Although various methods have been proposed to develop ABCs for nonlocal and PD models, recurring difficulties are apparent. Specifically, the implementation of ABCs often results in substantial demands of computational resources that complicate the application in practice. Additionally, the development of ABCs frequently involves the application of Fourier and Laplace transforms, rendering the computation of kernel functions especially complicated for PD wave-type models due to their inherent nonlocality. 22 In the present study, we propose a new and simple strategy to construct nonlocal ABCs for PD scalar wave-type equations based on prior work by the authors. Originally developed for the classical (local) wave equation and elastodynamics, 35,36 it has been demonstrated that the method we propose herein can be adapted and successfully applied to PD diffusion and corrosion equations of up to three dimensions (3D). 37,38 The present work emphasizes the need to consider the damping coefficient and construct reliable ABCs for scalar wave propagation problems in viscous media, a topic that is still scarce in related literature. In the proposed method, the near-field solution may be obtained using standard approaches, such as the meshfree PD solver or, when applied to the classical wave equation, the finite element method (FEM), and the ABCs are imposed on a layer of nodes at the truncating boundary. The truncating boundary is given by a surface boundary in the case of classical (local) equations and a volumetric boundary layer in the case of their nonlocal counterparts. The far-field solution is approximated using semi-analytical solutions based on harmonic exponential basis functions (EBFs) that satisfy the governing equation of the far-field; to get more insight into the features of EBFs, the reader may consult References 39-43. Using EBFs as time-dependent fundamental solutions (or plane-wave modes), we demonstrate how consistent dispersion relations may be obtained both at the continuum level and with the discretization approach employed in the near-field. In this way, the obtained ABCs are consistent with the discretization approach used for the near-field solution. The imposition of ABCs is performed numerically using a simple collocation procedure over an absorbing boundary layer. By selecting a minimal number of suitable and symmetrically distributed nodes, we demonstrate that numerical stability may be achieved while maintaining computational efficiency. In summary, the proposed method to construct ABCs has the following advantages: (i) it is consistent with the numerical dispersion relation of the near-field solution, (ii) it is easy to implement since the ABCs are of Dirichlet type and thus can simply be imposed as nodal values, (iii) it is free of any Fourier or Laplace transforms and is developed directly in the time and space domains (crucial for nonlocal models), and (iv) unlike high-order ABCs procedures, it does not involve higher-order derivatives or the introduction of auxiliary variables.
A brief outline of the present study is as follows. Section 2 describes the governing equations, starting from the classical scalar wave equation and proceeding towards the PD scalar wave-type equation. Section 3 derives the fundamental solutions required both for the development of the far-field solution and, to demonstrate that the PD scalar wave-type equation is a proper generalization of the classical scalar wave equation, the convergence of the PD fundamental solutions to the classical ones in the limit of vanishing nonlocality. The solution strategies for both the near-field and the far-field are discussed in Section 4. The effectiveness of the proposed ABCs in terms of numerical accuracy and stability is investigated with a number of selected examples in Section 5. Finally, conclusions and insights about the present study are summarized in Section 6.

PROBLEM DESCRIPTION
The following sections outline the governing equations that are the subject of the present investigation. In the general case, a viscous medium is assumed with velocity-dependent damping. The analysis is conducted in 2D.

Classical scalar wave equation
We start our discussion with the classical (local) scalar wave equation. The unbounded medium surrounding the near-field is assumed viscous and isotropic. Furthermore, as illustrated in Figure 1 (left), this medium contains the physical objects, like baffles and sources. Assuming an unbounded domain Ω ∈ R 2 , the classical scalar wave equation is where u(x, t) is the unknown field variable at point x = (x, y) T and at time t,u(x, t) andü(x, t) are the first and second time derivatives of u(x, t), respectively, 2 is the propagation wave speed, d is the damping coefficient, f (x, t) is a given source (or force) function, and ∇ 2 is the Laplace operator. Note that the factor 2 in the damping term is inserted by convention to simplify subsequent calculations. On the surface of the physical objects (scatterers) within Ω, appropriate Dirichlet and Neumann boundary conditions must be imposed as follows: where n = (n x , n y ) is an outward unit vector normal to the boundary surface, and f D (x, t) and f N (x, t) are functions that prescribe the values of the boundary conditions at the Dirichlet and Neumann boundaries, Γ D and Γ N , respectively. The governing equation in (1) is supplemented with the following initial conditions: where u 0 (x) andu 0 (x) are given functions that prescribe the first and second kind of initial conditions, respectively. The major challenge for solving the problem given by (1)-(3) arises from the unbounded medium. Following ABC methods, the domain Ω is divided into a bounded region Ω N (near-field) and an exterior region Ω F ∶= Ω ⧵ Ω N (far-field). The truncating boundary Γ ∞ specifies the regions, as depicted in Figure 1 (right). The near-field, Ω N , serves as the F I G U R E 1 Schematic representation of (left) a generic unbounded domain, Ω, including a scatterer (light blue region) and (right) the computational domain or near-field, Ω N , the far-field, Ω F , and the truncating boundary, Γ ∞ , for the classical scalar wave equation. The Dirichlet and Neumann boundaries are denoted, respectively, by Γ D and Γ N .
computational domain and encompasses all the physical objects: baffles and sources. This is required due to the radiation condition in scalar wave propagation, which mandates that all radiated energy must scatter to infinity, and thus, all source terms need to be included in the near-field to comply with this condition. 35 The homogeneous classical scalar wave equation (with homogeneous initial conditions) describes the far-field: Appropriate ABCs have to be imposed on the truncating boundary Γ ∞ , also referred to as absorbing boundary, also referred to as the absorbing boundary, to obtain the solution of the original problem within the near-field. This requires the transmission of any wave incident at the boundary to the far-field without any reflection back to the near-field.

Peridynamic scalar wave-type equation
This section is devoted to the (nonlocal) PD scalar wave-type equation, which may be regarded as a generalization of the classical equation outlined in the previous section. 44,45 We start by considering a region  x around the point x, which is called the neighborhood of x and is assumed to be circular in 2D with the radius called horizon: The point x interacts (nonlocally) with its neighboring points x ′ , as depicted in Figure 2, and the PD scalar wave-type equation isü where c is the micromodulus constant, the scalar-valued function (|| ⋅ ||) is referred to as the influence function, 46 and p is a parameter usually taken as integer, p ∈ {0, 1, 2, 3}. 47 The influence function is often defined as either constant or linear within the neighborhood  x . 48 By defining the bond ∶= x ′ − x, the constant influence function may be written as follows: F I G U R E 2 Schematic representation of the computational domain or near-field, Ω N , the far-field, Ω F , the neighborhood,  x , and the absorbing boundary layer, Γ * ∞ , for the PD scalar wave-type equation. The blue layer around the scatterer (light blue region) represents the nonlocal boundary where Dirichlet-or Neumann-type nonlocal boundary conditions are imposed.
which is independent of the bond length. Similarly, the linear influence function, which depends linearly on the bond length, may be expressed as follows: The micromodulus constant c, which generally depends on the horizon , must be determined with respect to the selected influence function such that (6) recovers the classical Equation (1) in the limit of vanishing nonlocality, that is, for → 0, under suitable regularity assumptions on the field variable u.
Analogously to the previous section, the homogeneous PD scalar wave-type equation (with homogeneous initial conditions) describes the far-field: In contrast to the classical equation, in the case of the PD scalar wave-type equation, an absorbing boundary layer Γ * ∞ has to be imposed rather than an absorbing boundary surface, as depicted in Figure 2.

FUNDAMENTAL SOLUTIONS
In this section, we develop fundamental solutions, hereafter referred to as modes, of the governing equation in the far-field using EBFs. For this purpose, we first consider a generic mode (x, t) given by with ∶= ( x , y ), where x , y , ∈ R and is the imaginary unit defined by 2 = −1. Without loss of generality, we can express as where controls the wave spatial fluctuation and determines the wave direction. Inserting the general mode from (10) into the far-field classical scalar wave Equation (4), we geẗ which results in A solution is obtained when the expression within square brackets vanishes, that is, The solution of the quadratic equation with respect to yields the local characteristic dispersion relation: Substituting (15) into (10) yields modes which play the role of semi-analytical solutions for the classical scalar wave equation: Similarly, inserting the generic mode (10) into the far-field PD scalar wave-type Equation (9), yields: Let the neighborhood with horizon around the origin be defined by Then, performing some algebraic manipulations, including the change of variable x ′ = x + , we obtain Once again, a solution is obtained when the expression within the square brackets vanishes, that is, In order to evaluate the integral term, we employ polar coordinates, so that = r(cos( ), sin( )), followed by a Taylor expansion of the integrand, as follows: where ∶= ⋅ (cos( ), sin( )) = cos( − ). Assume p = 2* in the denominator and a constant influence function (cf . (7)). Under these assumptions, and using the expansion in (21), we can write (20), as follows: The solution of the quadratic equation with respect to yields the following nonlocal characteristic dispersion relation: For the particular choice of micromodulus constant c = 4 2 2 , the above equation can be expressed as which in the limit as → 0 recovers the local characteristic dispersion relation (15): Inserting (24) into (10) yields the mode of the PD scalar wave-type equation for p = 2 given a constant influence function : Consider now a linear influence function (cf . (8)) while keeping the choice of p = 2. Using the same procedure as above yields the following equation (compare with (22)): which when solved for results in the following nonlocal characteristic dispersion relation: Choosing the micromodulus constant as c = 12 2 2 , the above equation can be written as which in the limit as → 0 also recovers the local characteristic dispersion relation (15). Insertion of (29) into (10) yields the mode of the PD scalar wave-type equation for p = 2 given a linear influence function : In the next section, we discuss the proper construction of semi-analytical solutions of the far-field using these modes.

Semi-analytical far-field solution
We begin by considering a node x i on the absorbing boundary (or boundary layer in the PD case) along with its cloud Ω ∞ i ; the situation is illustrated in Figure 3. Each cloud holds its own local Cartesian coordinate system (x, y), where F I G U R E 3 Schematic representation of a portion of a solution domain for the PD scalar wave-type equation near the absorbing boundary layer, Γ * ∞ , and a cloud, Ω ∞ i , centered at a node x i on the absorbing boundary layer with a rotated local coordinate system (rotated by an angle ) for the far-field approximation.
the positive x-axis direction is oriented towards the far-field. We consider the coordinate transformation of a representative vector x with respect to the local coordinate system with origin at x i via the 2D rotation matrix R into the local coordinates x: We denote the value of the field variable within the cloud Ω ∞ i by u ∞ (x, t) and approximate it via a linear combination of modes as follows: where a k,l and b k,l are unknown coefficients, n and n specify the number of modes, where the first one indicates the number of wave directions and the latter one the number of fluctuations of the wave in space, and̃l =̃( l ) is defined by the characteristic dispersion relation of the corresponding governing equation, as discussed in the previous section, that is, = d ±̃( ), (cf . (10)). For instance, in the case of the classical wave equation, (15) yields the following expression:̃l Note that the value of the damping coefficient d cannot be chosen arbitrarily high to ensure physical feasibility, and it is constrained by the characteristic dispersion relation; for instance, in (33), d 2 ⩽ 2 l 2 , ∀ l . The values of k and l are selected from two symmetric intervals: where Δ and Δ are parameters to be specified; a rationale for the selection of the respective numerical values is discussed in References 35-37. Given that the choice of these parameters may impact numerical stability, the selection should be made with care; this is addressed in Section 4.3.

F I G U R E 4
The propagation region of outcoming modes with respect to ; recall that is within the interval Δ [−1, 1].
Note that at Δ = 0 the incident angle of the respective mode is parallel to the positive x-axis. By considering a small region around Δ = 0, it becomes evident that the first summand in (32) represents outcoming waves, whereas the second one represents incoming waves. For the construction of the ABCs, the incoming waves from the near-field have to be absorbed at the truncating boundary. Hence, the coefficients b k,l in (32) must be set to zero, resulting in the following expression: Using polar coordinates with (x, y) = r(cos( ), sin( )), − ⩽ ⩽ , yields the equivalent representation The above equation describes outcoming waves from the near-field, for which the following condition holds: which may be reduced to the following relation: The condition in (38) results in a half-space centered at the point x i that forms the angle k with the local x-axis, i.e., the outward boundary normal. An increase (or decrease) in the angle k causes a rotation of the half-space around the point x i , as illustrated in Figure 4. For the construction of the ABCs, we adopt a symmetric interval in (34) to specify the angle k in order to transmit incoming waves from the near-field at different incident angles into their respective half-space. In addition, the value of Δ must not be selected arbitrarily high, such that incident waves propagate out and the energy is absorbed.

SOLUTION STRATEGY
This section provides a description of the numerical approximations employed in the present work. First, in Section 4.1, we review the procedure for determining the numerical solution in the near-field, where standard methods for both the classical scalar wave equation and PD scalar wave-type equation are described. Then, in Section 4.2, we turn to the central subject of the present work, namely the numerical construction of ABCs, which are derived from semi-analytical solutions of the far-field. Note that, as opposed to the classical equation that involves a boundary surface (see Γ ∞ in Figure 1), for the PD scalar wave-type equation we employ a boundary layer (see Γ * ∞ in Figure 2). Finally, in Section 4.3, we describe the numerical implementation and related steps concerning numerical stability.

Near-field
The near-field solution is discretized in time and space, so that the solution is approximated for time instants t n , n = 0, 1, 2, … , N, at the nodes x i , i = 1, 2, … , K, with Δt a fixed time step. In the course of the present work, we employ the explicit velocity-Verlet scheme for time integration, which, given the field variable, u n i , and its first and second time derivatives,u n i andü n i , respectively, at node x i and time t n , may be expressed as follows: By (39), the expressions for the field variable and its first time derivative at node x i and time t n+1 may be reformulated, respectively, as follows: Let us consider first the classical scalar wave Equation (1) in the near-field. Performing the spatial discretization by means of the FEM, the weak form may be expressed in terms of the following linear system of equations: where M is the mass matrix, K is the stiffness matrix, D ∶= 2dM is the damping matrix, F n is the force vector on the right-hand side at time t n , andÜ n ,U n , and U n are the nodal acceleration, velocity, and displacement vectors at time t n .
In analogy to (40), the time-stepping procedure for the nodal vector U n+1 and the first time derivativeU n+1 at time t n+1 is given by For the spatial discretization of the PD scalar wave-type equation, the most widely used approach, also referred to as the standard scheme, 50 is a meshfree method 51 that can be adopted to solve the strong form of (6). In the standard scheme, the computational domain is partitioned into a regular cubic grid in 3D or square grid in 2D of cells with grid spacing Δx whose centers represent the computational nodes; this is illustrated in 2D in Figure 5. Each node x i is associated to a volume V i = (Δx) 3 in 3D or an area (Δx) 2 in 2D of its respective cell. Figure 5 further depicts the neighborhood  x i in 2D with horizon around node x i , wherein nonlocal interactions are established between x i and its neighboring nodes x j ∈  x i , which are refereed to as its family nodes. The set of family nodes associated to x i is denoted by  i . Thereby, the integral in (6) is expressed as a summation over neighboring cells via a one-point quadrature over each of those cells. 52 Using this procedure, the discrete form of the PD scalar wave-type Equation (6) can be expressed as follows:ü where f n i is the right-hand side of (6) evaluated at node x i and time t n and the function (x j − x i ) is called the volume reduction factor, 53 which is used to correct the volume V j covered by the neighborhood of x i ; algorithms for volume correction (also referred to as partial-volume algorithms) are presented in References 52,54,55. F I G U R E 5 Schematic representation in 2D of the standard discretization scheme in PD and the neighborhood  x i associated to node x i .

Far-field
The ABCs are derived from the semi-analytical solution of the far-field, whose unknown coefficients are obtained using a simple collocation procedure inspired by the weighted least squares (WLS) approach commonly employed in the finite point method (FPM). 56,57 Since the procedure is essentially identical for both the classical scalar wave equation and PD scalar wave-type equation, no distinction between the two equations is made at this point. As a result, we obtain ABCs of Dirichlet-type, applicable to the nodes at every time instant by means of constant updating vectors. We begin by approximating the field variable and its first time derivative at the nodes within the cloud, x ∈ Ω ∞ i . In contrast to previous works, 35,37 where the time interval [t n , t n + Δt] was used, we assume the symmetric interval [t n − Δt, t n + Δt], in which the approximation is considered to be valid. If we reset the time and define the approximation in the local time coordinate, t ∈ [−Δt, Δt], we can express the approximation around the time instant t n as follows: where u n ∞ (x, t) andu n ∞ (x, t) are the field variable and its first time derivative, respectively, at node x and time t (in the local time), m anḋm are the mth mode and its first time derivative, respectively, n b is the number of modes, and p n m and q n m are constant unknown coefficients corresponding to the presumed local time interval around t n . Notice here that we assume two different sets of coefficient vectors p n and q n , such that the approximation of the field variable u n ∞ and its first derivativeu n ∞ can be performed independently of each other. This is possible as the construction of ABCs in the following involves a Least Squares approximation of the coefficient vectors, which in general is not unique. In the rightmost expressions in (44), we employ, for simplicity, vector notation, where , Next, we define the following vectors: that contain the field variable as well as its first time derivative at the local time steps t = 0 and t = −Δt. Using this notation, we determine the unknown coefficient vectors p n and q n in the following steps. For that purpose, we arrange all the cloud nodes in vectors G n and H n , according to (46), as follows: where M G and M H denote the moment matrices of the collocation procedure. Assuming that the number of nodes in the cloud, n c , is less than the number of modes, n b , we can obtain the coefficient vectors p n and q n by where M + G and M + H denote the Moore-Penrose generalized inverses of the respective moment matrices. Substituting the expressions on the right-hand sides in (48) for the coefficient vectors in (44), the following expressions are obtained: To determine the nodal values for the field variable u n+1 i and its first time derivativeu n+1 i at the central node of the cloud x i , we simply evaluate the equations in (49) at this node and at the local time t = Δt: where V G and V H are constant updating vectors corresponding to node x i . Since the nodal values that conform G n and H n on the right-hand sides of (50) are known at the time instant t n , we obtain Dirichlet-type ABCs that are imposed on the absorbing boundary nodes at each time instant. As a result, this yields a computationally efficient extrapolation in time procedure, since the application of the ABCs merely requires the nodal values because the updating vectors can be precomputed.

Numerical implementation
This section discusses important considerations regarding the implementation of ABCs in terms of both compatibility with the method used for the near-field solution and numerical stability. We reconsider the nonlocal characteristic dispersion relations in (24) and (29). These equations were derived at the continuum level and thus do not accurately represent the dynamics of discrete systems, which are characterized by discrete dispersion relations. Therefore, in order to construct modes according to the discrete dispersion relation, we employ the one-point quadrature of the PD standard scheme to discretize the spatial integral in (20) instead of resorting to the Taylor expansion in (21): where  0 denotes the set of family nodes around the origin. The solution that makes the quadratic expression with respect to vanish on the left-hand side of (51) is given by the following nonlocal numerical dispersion relation: which is consistent with the PD standard scheme. Note that this procedure may be employed in an analogous manner for other discretization schemes used for the near-field solution. We emphasize that, in contrast to the previous work 35 that focused on local equations, the present approach involves the selection of modes that match the dispersion relation of the discrete nonlocal system, which proves to be essential for numerical stability and accuracy of nonlocal wave problems. Figure 6 compares the nonlocal numerical dispersion relation (52) for a choice of = 2, for the case of negligible damping, that is, d = 0, and constant influence function (cf . (7)), with the corresponding nonlocal characteristic dispersion relation based on a Taylor expansion (24) and the local characteristic dispersion relation (15). For the nonlocal characteristic dispersion relation, we show two cases: one based on a lower-order (fourth-order) Taylor expansion and one based on a higher-order (26th-order) Taylor expansion. Similarly, for the nonlocal numerical dispersion relation, we show two cases: one with a coarser discretization (Δx = ∕2) and one with a finer discretization (Δx = ∕8). In all cases, the nonlocal (continuum and discrete) dispersion relations approach the linear, local dispersion relation in the limit as ( ) → 0, while exhibiting substantial deviations from it for large values of , as expected. The results also demonstrate that a high-order Taylor expansion or a fine discretization are needed to accurately recover the nonlocal characteristic dispersion relation; the two approaches seem to give a similar result. Figure 6 further illustrates that the discrete dispersion relation differs from a continuum-level Taylor series approximation. Thus, in order to construct compatible and accurate ABCs, it is necessary to match the dispersion relations of the modes employed in the far-field to the discrete near-field dynamics. The neighboring points as well as the mesh size and horizon employed in evaluating the nonlocal numerical dispersion relation in (52) must therefore match those used in the near-field, irrespective of the differences between the analytical and the discrete forms, in order to ensure that the dispersion relation of the modes used to construct the far-field corresponds to the dispersion relation of the (discrete) near-field.
Another important note, alluded to in Section 3.1, is the selection of the parameters Δ and Δ as well as the associated ranges in (34). The application of EBFs is comprehensively studied in Reference 40, and possible ranges are provided in References 35-37. In the present study, following (34), an interval of Δ [−1, 1] with ten equidistant subdivisions is employed. Empirically, it has been found that a value of Δ = ∕4 generally achieves satisfactory results. The influence of this parameter is presented by means of a sensitivity analysis in Reference 35. The second influential parameter Δ in (34) affects the spatial fluctuation of the modes and determines the range of the interval Δ [−1, 1]. To find an upper limit to the value of Δ , we adopt the Nyquist-Shannon sampling theorem (see Reference 35). Since a mode may only be resolved on a discrete lattice up to a certain maximum frequency, according to (34), the parameter Δ thereby determines the maximum frequency content. The sampling theorem thus provides the upper bound for Δ as follows: Similar to , we employ ten evenly spaced subdivisions on the interval Δ [−1, 1]. We note that the proposed extrapolation procedure for the application of the ABCs in (50) is similar to a fixed-point iteration scheme over large time periods. It is well-known that fixed-point iterations contract only if they are Lipschitz continuous. Similarly, although not a sufficient condition, for practical applications, we consider that the norm of the updating vectors V G and V H should not exceed two, since two time steps are used according to (46) for the extrapolation in (50), in order to ensure a stable simulation over the course of many time steps. Although the Euclidean norm may be adopted, in practice the maximum norm ||V|| max ∶= max i=1, … ,N |V i | for a vector V ∈ R N is found to be preferable. In the proceeding numerical examples, Δ is chosen such that the latter condition holds.
The next issue we address concerns numerical stability, which is closely related to the selection of collocation nodes in (47) and the solution procedure involving the Moore-Penrose inverse matrices in (48). It is found, in practice, that the selection of all family nodes in the cloud, Ω ∞ i , as depicted in the top-left illustration in Figure 7, leads to numerically unstable simulations. This problem may be explained by the fact that the extrapolation on the central node x i incorporates information from nodes that are located spatially ahead of it. Within the context of the PD scalar wave-type equation, this would suggest that the temporal prediction of the outcoming wavefront on the central node x i includes information from nodes to which the wavefront is yet to arrive. With that in mind, we introduce the extrapolation region  ∞ i (Δ ) ⊂ Ω ∞ i of the absorbing node x i as a conic subset of the cloud as follows: where Ω ∞ i denotes the cloud of node x i , depicted in the top-right illustration in Figure 7. Here, the angle Δ determines the spread of the cone and, according to (34) and the previous discussion, is also set to ∕4, which is assumed to be the maximum angle of incident waves and in practice yields reasonable results. The other important concern is the choice and number of collocation nodes employed in the extrapolation region, which is crucial for high resolutions or for cases with a large horizon relative to the grid spacing. It is a well-understood fact that meshfree methods should use as few nodes as possible for computational efficient while attaining the smallest possible error. 58 Moreover, a minimum number of nodes with adequate distances to one another is preferable, as it improves the conditioning of the moment matrices in (48) and thus ensures the numerical stability of the simulation. In related literature, a number of techniques are proposed to determine an optimal set of collocation nodes, often summarized under the umbrella term greedy sparse methods. 58,59 However, these methods are, in many cases, iterative and are likely to require excessively high computational costs in practice. Therefore, we propose a different approach that: (i) proves to be straightforward to implement, (ii) improves the numerical stability of the simulation, and (iii) reduces computational effort. Inspired by Reference 50, we introduce a fixed set of auxiliary nodes within the cloud, which are depicted in blue in the bottom-left illustration in Figure 7. Collocation is then performed over the nearest nodes to the auxiliary nodes within the extrapolation region,  ∞ i (Δ ), as illustrated in the bottom-right illustration in Figure 7. In this way, the total number of nodes used for collocation is reduced, resulting in improved computational performance as well as better conditioning of the moment matrices. Note that, regardless of the number of auxiliary nodes, the collocation nodes should ideally be arranged symmetrically in order to avoid eventual instabilities. 58 Furthermore, it is worth noting that in contrast to methods such as the Finite-Point-Method and the related Weighted-Least-Squares, where a higher number of nodes in the cloud can contribute to higher accuracy of the solution since there are fewer basis functions available to approximate the field variables, our method (which uses exponential basis functions) requires a higher number of basis functions to effectively absorb incoming waves within a range of wave numbers and orientations. However, in the present context, a higher number and density of nodes in the cloud may deteriorate the conditioning of the moment matrices and render the method unstable.
Remark 1. To enhance the numerical stability of the collocation procedure, we suggest a practical and straightforward method to decrease the density of cloud nodes by iteratively reducing the number of nodes along the cloud coordinate axes by a factor of 2, as conceptually illustrated in Figure 7. While this may not necessarily yield an optimal choice of collocation nodes as described in References 58 and 59, it has been found to effectively improve the conditioning of the moment matrices in (47). To ensure numerical stability prior to simulation, we recommend monitoring the condition numbers of the moment matrices. Based on our numerical simulations and investigations, we recommend that the condition numbers fall within the range 10 11 − 10 12 to guarantee numerical stability. Consequently, we suggest decreasing the density of collocation nodes by a factor of two until the condition numbers fall within this specified range. This approach provides a user-friendly and practical means of ensuring numerical stability prior to simulation, especially for clouds with a large number of nodes that may lead to bad condition numbers. While this may not be an optimal choice of collocation nodes, it has been found to be effective in practice.

NUMERICAL EXAMPLES
In this section, we present three numerical examples to demonstrate the accuracy and performance of the proposed ABCs. The first example illustrates the application of the ABCs for the classical (local) scalar wave Equation (1), while the other two examples demonstrate application of the ABCs for the (nonlocal) PD scalar wave-type Equation (6). In all the examples, a viscous medium with non-vanishing damping coefficient is assumed. The model parameters are specified in each example. For comparison purposes, in the figures below, reference solutions determined on relatively large computational domains are depicted within the near-field of the comparative truncated models equipped with ABCs.

Example 1
In the first example, we demonstrate the application of the proposed ABCs for the classical (local) scalar wave equation (1). We employ the FEM with linear shape functions to obtain the near-field solution. In Figure 8, simulation results are presented at different time instances. The near-field given by the square-shaped domain (−26, 26) × (−26, 26) around the origin is depicted, which is discretized with square elements of edge length Δx = 0.5. We present the reference solution in the left column, which was determined with homogeneous Dirichlet boundary conditions on a much larger domain given by (−126,126) × (−126,126), and the near-field solutions obtained when applying the proposed ABCs (middle column) and the conventional ABCs of first order (right column). We refer to the models used for the left, middle, and right columns as reference, proposed, and conventional, respectively. The interested reader will find a derivation of the conventional first-order ABCs for the classical scalar wave equation in viscous media in Appendix A.
Since the equation of consideration is local, the ABCs in the proposed model are implemented only on the outermost boundary nodes. The present discretization results in a total of 11,025 nodes and 10,816 elements in both the proposed and conventional models, whereas a total of 251,001 nodes and 250,000 elements are used in the reference model. Each model employs the explicit velocity-Verlet scheme from (42) for time integration with a time step of Δt = 5 × 10 −4 . The wave propagation speed is assumed to be 2 = 1, and the medium damping coefficient is set to d = 0.0125. We set the parameter Δ from (34) to 1, which was found empirically to yield reasonable simulation results as well as adequate norms of the updating vectors in (50). Initial conditions according to (3) are given as a pulse with u 0 (x) = exp In Figure 8, we compare the solution between the models in the near-field at four selected time instants, at which the wave reaches the boundary of the near-field as well as subsequent wave states. It is found that the proposed model is in good agreement with the reference model, even when the wave reaches the boundary, allowing the wave to be transmitted out without any significant reflections back into the computational domain. In contrast, the conventional first-order ABCs result in considerable reflections, which are clearly noticeable at t = 40 and t = 50. These reflections become even more apparent when examining the absolute value of the difference between the near-field solution, u(x, t), and the reference solution, u ref (x, t), that is, |u(x, t) − u ref (x, t)|, for the final time step t = 50 in Figure 9. Here, the solution difference for the proposed model is depicted in the left column, while the solution difference for the conventional model is presented in the right column. We note that, in the conventional model, the reflections appear to emanate substantially from the near-field corners, which is consistent with the fact that the corners in this model require a special consideration as they may lead to numerical instabilities. 14 To demonstrate that the proposed ABCs method is numerically stable, we provide Figure 10 (right). Here, the course of the normalized energy within the near-field, with respect to the initial energy of the system, is plotted for a total of where Ω N denotes the integration region. Throughout the overall time period, no unstable behavior is observed. In addition, the accuracy of the proposed model compared to the conventional model is illustrated in Figure 10 (left). The time evolution of the near-field solution at a selected node is plotted for all three of the models employed. As in Figures 8 and 9,

F I G U R E 9
The absolute difference between the solutions of (left) the proposed and reference models and (right) the conventional and reference models at different time instants in Example 1.
it is evident that the solution of the proposed model is in good agreement with the reference solution, whereas the solution of the conventional model exhibits significant reflections.

Example 2
In the second example, we demonstrate the application of the proposed ABCs for the PD scalar wave-type Equation (6), for the case of constant influence function (cf . (7)) and p = 2. We apply the standard scheme to discretize the governing equation in the near-field. In this case, the near-field is given by a circular domain with radius R = 26, which is depicted in Figure 11. We compare the solution of three models at different time instants. The computational domain is discretized F I G U R E 10 Left: Time evolution of the solution at x = (20, 10) (see Figure 8) for the different models in Example 1. Right: Comparison of the normalized energy evolution in time, where E 0 = E(0) is the initial energy of the system, between the proposed and reference models in Example 1.
by means of regular square cells of edge length Δx = 0.25, and a horizon of = 6Δx is assumed. The reference solution presented in the left column in Figure 11 was determined on a much larger domain (R = 126) with homogeneous Dirichlet-type boundary conditions. The middle and right columns in Figure 11 depict the proposed model with two different choices of dispersion relation. The model used for the middle column employs the nonlocal numerical dispersion relation (52) and is referred to as Proposed-A. The model used for the right column, in contrast, employs the local characteristic dispersion relation (33) and is referred to as Proposed-B. In both the Propose-A and Proposed-B models, the ABCs are imposed on a layer of width equal to the horizon, and we reduce the number of collocation nodes according to Figure 7 for overall enhanced conditioning and computational efficiency. The discretization employed herein results in a total of 33,949 nodes for both the Proposed-A and Proposed-B models, while the reference model uses 797,889 nodes. We employ the explicit velocity-Verlet scheme for time integration, which follows from (40) and (43), with a time step of Δt = 1 × 10 −4 . The wave propagation speed is set to 2 = 1, and the medium damping coefficient is set to d = 0.0125. As in Example 1, we choose the parameter Δ from (34) as Δ = 1. For the initial conditions, a pulse given by u 0 (x) = exp(− log(2)(||x||)) andu 0 (x) = 0 is assumed.
In Figure 11, we compare the solution between the models in the near-field at four selected time instants, demonstrating good agreement between the models' solutions. Both the Proposed-A and Proposed-B models absorb the incident wave. However, the Proposed-A model exhibits less reflections compared to the Proposed-B model, which can be noticed at times t = 40 and t = 50, once the wave departs the near-field. In Figure 12, this becomes more apparent, where the absolute value of the difference between the solution of the proposed model and the reference solution is depicted; this is reported for the Proposed-A model in the left column and for the Proposed-B model in the right column. This result is expected since the dispersion relation of the EBFs employed in the Proposed-A model, that is, the nonlocal numerical dispersion relation, is chosen to match the numerical method in the near-field. On the contrary, the local characteristic dispersion relation used in the Proposed-B model significantly differs from the nonlocal numerical dispersion relation, thereby resulting in the observed reflections. In Appendix B, the interested reader may find an additional numerical example, which further illustrates the performance of the proposed approach for highly dispersive nonlocal wave propagation with large horizons relative to the grid spacing.

Example 3
In the third example, we extend the application of the proposed ABCs from Example 2 to a more complex problem involving two notches in the near-field, as depicted in Figure 13. The first notch is given by the line segment between the points (−7, −7) and (7, −7), while the second notch is given by the line segment between the points (0, 10) and (7,10). The notches are implemented by breaking all the bonds that cross the corresponding line segments. Except for the horizon, which is assumed to be = 3Δx in this example, the model parameters are the same as in Example 2. The initial conditions are assumed the same as in Example 1. Here, we omit the comparison with the Proposed-B model (see Example 2) and simply compare the solution of the Proposed-A model (see Example 2) with the reference solution. The presence of the notches leads to reflections within the near-field, resulting in a complex wave state consisting of many superimposed comparatively large angles of incidence on the absorbing boundary. The solution accuracy relative to the reference solution would improve if the angles of the incident waves relative to the boundary normal vectors are reduced, which may be achieved by increasing the radius of the near-field thus increasing the distance between the absorbing boundary and the notches. For the case of a larger near-field with domain radius R = 42, Figure 14 shows an improvement in the accuracy of the solution of the Proposed-A model relative to the reference solution. This fact thus illustrates that there is always a trade-off between numerical accuracy and computational efficiency related to the selection of the domain size of the near-field. Remark 2. This numerical example poses a complex problem involving two notches in the near-field, where wave reflections can result in waves reaching the absorbing boundary from multiple incident angles. One approach to improving numerical accuracy in such cases is to enlarge the numerical domain, as depicted in Figure 14, but this can come at a significantly higher computational cost. Another parameter that can potentially affect the accuracy of the results is Δ from (34), which determines the range of angles of the modes used to construct the ABCs. However, increasing Δ may not necessarily lead to better outcomes, and even result in numerical instability due to wave reflections from the absorbing boundary layer to the near-field. 35 Hence, a careful balance between numerical accuracy and computational cost is necessary, especially for complex problems, like the one described here.

CONCLUSIONS
This paper introduces a new approach to tackle the construction of absorbing boundary conditions (ABCs) for the (nonlocal) peridynamic (PD) scalar wave-type equation in unbounded viscous media. The proposed ABCs are of Dirichlet type, which facilitates their application since their implementation neither depends on differentiation procedures nor uses auxiliary variables which are not straightforward to incorporate even for the case of the classical (local) scalar wave equation. The ABCs are derived from a semi-analytical solution of the exterior domain (far-field) consisting of plane-wave modes. In this way, the modes satisfy the nonlocal dispersion relation corresponding to the PD scalar wave-type equation and converge to the modes corresponding to the local equation in the limit of vanishing nonlocality. The ABCs are entirely developed in the time and space domains without the utilization of Fourier or Laplace transforms, which proves to be complex for nonlocal models. Our numerical studies show that, at the discrete level, the far-field solution is obtained consistent with that of near-field, in terms of the dispersion relation, provided that the modes are calculated using the same numerical dispersion relation of the near-field. The unknown coefficients of the far-field semi-analytical solution are calculated through a simple collocation procedure in space and time. A comprehensive discussion on efficient derivation of the unknown coefficients is presented. The performance of the proposed ABCs in terms of accuracy and stability is shown through several numerical examples in 2D. Our investigation demonstrates that the proposed ABCs yield a proper level of accuracy even in problems characterized by highly-dispersive propagating waves. Extension of the proposed approach for PD elasticity equations is a subject of future study.

APPENDIX A. CONVENTIONAL ABSORBING BOUNDARY CONDITIONS
The derivation of the conventional first-order ABCs for the classical (local) scalar wave equation in viscous media follows. 60 We assume that the local coordinate system at each point of the absorbing boundary Γ ∞ is oriented so that the exterior domain lies in the half-space x ⩾ 0, where the x-axis is normal to the absorbing boundary and points towards the exterior domain. To this end, the governing equation of the far-field reads as follows: where u xx and u yy denote second partial derivatives with respect to the spatial coordinates in the x-and y-directions, respectively. The initial conditions for the exterior domain are: u(x, 0) = 0, x ⩾ 0, y ∈ R.
In order to apply the Fourier transform to the field variable u in time, we have to extend it to negative times by setting the corresponding values for t < 0 to zero. We denote the prolonged field variable for negative and positive times asũ. As the initial conditions are zero according to (A2) for x ⩾ 0, the following homogeneous wave equation is obtained forũ: Then, we apply the Fourier transform in the time t and y-direction to obtain the ordinary differential equation in the dual coordinates (k, ) as follows: whereû denotes the Fourier transform of the prolonged field variableũ. For a fixed tuple (k, ) in the Fourier space, a general solution for (A4) is given in the following form: where j (k, ) and j (k, ) denote unknown coefficients that generally depend on the dual coordinates k and . By inserting the exponential ansatz function from (A5) in (A4), the following relation is obtained for a representative pair of coefficients and k: Defining the auxiliary variable z ∶= k , the two solutions of the quadratic Equation (A6) can be expressed as For the derivation of the first-order ABCs, we perform a first-order Taylor series expansion in (A7) for z around z = 0. In doing so, we obtain the following approximation: . (A8) Note that higher-order Taylor series expansions would likewise be conceivable; however, the subsequent equations turn out to be increasingly complex. Inserting the two approximate solutions from (A8) in (A5) and using the inverse Fourier transform, a general solution forũ reads as follows: where we distinguish the positive and negative solutions of (A8) in terms of the corresponding coefficients + and − , respectively. Since the term d is nonnegative and the solutionũ is bounded, the first coefficient must be zero, that is, + ≡ 0. Therefore, we can simplify (A9) and express it as follows: By evaluating (A10) at x = 0, we recognize that − (k, ) =û(0, k, ). Subsequent differentiation of (A10) with respect to x and moving the expression on the right-hand side to the left-hand side yields for x = 0: where we emphasized the explicit variable dependencies of the prolonged field variableũ x and its Fourier transformû in the notation. In Fourier space, (A11) reads as follows: u x (0, k, ) +û(0, k, ) + dû (0, k, ) = 0.
Multiplying (A12) by and applying the inverse Fourier transform finally yields the first-order approximation to the exact ABCs asu (0, y, t) + u x (0, y, t) + du(0, y, t) = 0, which is a partial differential equation that has to be enforced on the absorbing boundary Γ ∞ (at x = 0). Note that we omitted the previous notation of the prolonged solutionũ as we only consider the solution for t ≥ 0. Note that, for a vanishing damping coefficient d = 0, the relation in (A13) corresponds to the conventional first-order ABCs provided in Reference 60.

APPENDIX B. NONLOCAL WAVE PROPAGATION IN 1D: NUMERICAL PERFORMANCE FOR HIGHLY DISPERSIVE WAVES
The purpose of this section is to investigate the numerical performance of the proposed ABCs for highly dispersive nonlocal wave propagation on an unbounded viscous domain. As the waves involved in the following numerical example contain high-frequency components, a sufficiently fine spatial and temporal resolution is required. Thus, we perform the simulations in 1D, as the construction of the numerical reference solution is computationally costly in higher dimensions. monitoring the condition numbers and updating vector norms in order to guarantee numerical stability prior to the simulation and, if necessary, reducing the number of cloud nodes until the condition numbers and updating vector norms fall within the respective ranges provided in Section 4.3.
The results of the reference solution as well as the truncated domain equipped with the proposed ABCs are depicted in Figure B1. It can be seen that, for the chosen horizon = 0.25 and the given initial conditions, the waves draw high-frequency tails as they propagate to each end of the computational domain. 61 Furthermore, a good match between the solution of the truncated domain equipped with the proposed ABCs and the reference solution can be observed, demonstrating the effectiveness of the proposed method on applications involving highly dispersive waves and large values of m. The proposed method maintains numerical stability throughout the course of the simulation.