Dirichlet absorbing boundary conditions for classical and peridynamic diffusion-type models

Diffusion-type problems in (nearly) unbounded domains play important roles in various fields of fluid dynamics, biology, and materials science. The aim of this paper is to construct accurate absorbing boundary conditions (ABCs) suitable for classical (local) as well as nonlocal peridynamic (PD) diffusion models. The main focus of the present study is on the PD diffusion formulation. The majority of the PD diffusion models proposed so far are applied to bounded domains only. In this study, we propose an effective way to handle unbounded domains both with PD and classical diffusion models. For the former, we employ a meshfree discretization, whereas for the latter the finite element method (FEM) is employed. The proposed ABCs are time-dependent and Dirichlet-type, making the approach easy to implement in the available models. The performance of the approach, in terms of accuracy and stability, is illustrated by numerical examples in 1D, 2D, and 3D.


Introduction
Peridynamics (PD) is a recent nonlocal theory that has received a widespread attention in computational mechanics. It has been widely exploited to solve various problems in mechanics and physics. The theory was originally introduced by Silling [56] and Silling et al. [58] to handle material failure in solid structures, which is not an easy task for the classical continuum mechanics (CCM) theory. In fact, the original formulation of PD introduces an equation of motion in solid mechanics based on integro-differential equations rather than on partial differential equations (PDEs) which are undefined at discontinuities. Therefore, PD models are effective tools in modelling material discontinuities, and they consider fracture and damage as natural parts of the solution. To this respect, many PD models have so far been proposed and applied to a broad range of problems in solids including fracture and crack propagation; to mention a few, see [7,16,40,52,57,59,66] and for a comprehensive review one may refer to [11,17].
In nonlocal continuum theories, the value of a physical quantity at a given material point is influenced by quantities within a finite neighborhood around that point; we refer to such influences as nonlocal effects. PD models make use of a characteristic length scale called horizon that determines the region of the nonlocal interactions. Each material point establishes direct interactions with other points within its neighborhood. The interactions are governed by a response function that includes all the material constitutive information. As the horizon shrinks and asymptotically approaches zero, the interactions become local and the formulation reduces to that of the classical local theory when suitable regularity assumptions hold [19].
The focus of the present study is on diffusion-type problems. In fact, the solution of diffusion problems is of great importance for engineering and physics applications concerning heat and mass transfer. The mathematical description of many phenomena from areas such as fluid dynamics, chemistry, biology, information, environmental and materials sciences is governed by diffusion-type equations. In this regard, numerous numerical methods (based on the classical local diffusion) have so far been employed, for example the finite element method (FEM), the finite difference method (FDM), the boundary element method (BEM), and meshfree methods (see, e.g., [4,6,25,31]). At the macroscale, most diffusion processes can be described well by local models based on Fourier's law (heat conduction) as well as Fick's law (mass transport). However, at this scale, nonlocal effects may play a crucial role; for instance, in the case of heat conduction with steep temperature gradients, stochastic jump processes and anomalous diffusion are observed in heterogeneous environments [21,35,36,73]. Therefore, nonlocal diffusion equations as well as related numerical methods have received a considerable attention in the literature [29,35,42]. Recently, miniaturization of devices has provided more cases where the application of nonlocal models rather than local ones is more appropriate [61,67].
Unfortunately, many of the nonlocal models are inapplicable for problems in which discontinuities (whether strong or weak) in the system emerge, interact, and evolve. Thermal cracking [5,33,65], hydraulic fracturing [43], and pitting corrosion [28] are just a few examples, where part of the solution is governed by diffusion, and, at the same time, they are subjected to the emergence of spontaneous discontinuities.
To address also such problems, PD is a promising approach. The first PD diffusion model was developed by Gerstle et al. [24] for electromigration that accounts for heat conduction in a 1D body. Bobaru and Duangpanya proposed a 1D PD heat conduction model in [9], and an extension of their work to 2D bodies with evolving discontinuities is presented in [10]. In these cases, the interaction between material points is prescribed by a thermal response function, and the spatial derivatives of temperature, involved in classical diffusion equations, are replaced by spatial integrals of the thermal response function. Zhao et al. [72] introduced a PD model for transient advection-diffusion problems. The majority of these models utilize the bond-based version of PD. For a mathematical investigation into the generic forms of PD diffusion models one can refer to the studies conducted in [18,60]; moreover, the first attempts to derive the state-based formulation for thermal diffusion can be found in [42,46]. There is still ongoing research on PD diffusion models [14,64].
More recently, the application of PD to corrosion problems has received a considerable attention [15,28,29,44]. Corrosion problems mainly involve a diffusion process in an electrolyte/solid system coupled with a phase-change mechanism; this results in a moving interface problem. Treatment of moving interfaces by local models is challenging since the accurate position of the interface, where the flux continuity conditions (Stefan conditions) are satisfied, must be traced over time [15,22]. The situation gets worse for corrosion problems as the interfaces evolve and increase in number in time. PD models are effective tools for treating these problems since they simulate the evolution of the corrosion damage and material disintegration in a mathematically consistent manner. Moreover, the Stefan conditions are naturally satisfied by the PD solution [13]. PD modelling of such problems only requires modifying the nonlocal interactions between material points near the interfaces [13], and the autonomous evolution of moving interfaces is a natural part of the PD models.
All the aforementioned studies focus on the solution of bounded-domain problems. In contrast, the solution of diffusion-type problems in an infinite domain via both nonlocal and local models is the main concern of the present study. The application of unbounded-domain problems can be found abundantly in different fields, including fluid dynamics, heat conduction, biology, and finance [27,32]. Modelling diffusion in unbounded domains finds application in many practical problems regarding corrosion [12]. The main challenge in modelling diffusion-type problems in unbounded domains lies in the infinity of the medium; however, the standard numerical methods, such as the FEM, are developed for bounded-problem domains.
Artificial boundary methods (ABMs) are the most common way of coping with unbounded domains numerically [53,71]. The idea in ABMs is first to truncate the unbounded domain at an artificial boundary (or absorbing boundary) at a certain distance from the region of interest. Then, by constructing appropriate absorbing boundary conditions (ABCs), applied to the artificial boundary, one can convert the original problem into a bounded one. The exterior domain and the region of interest are referred to as far field and near field, respectively. One has to place the artificial boundary so that the near field includes all the sources. The numerical accuracy and stability strongly depend on the design of the ABCs. The reduced problem domain can reproduce the original solution (within the interior region) provided that accurate and reliable boundary conditions are imposed on the artificial boundary. In the last three decades, the derivation of ABCs, for local wave-and diffusion-type problems, has been the subject of numerous studies, yielding, e.g., the high-order ABCs [26,69] and the perfectly matched layer (PML) method [3,8]. However, the solution of nonlocal problems in unbounded domains has received little attention, and it is still an open subject of research [71]. In fact, construction of nonlocal ABCs is challenging since the interactions are nonlocal; in addition, the nonlocal operators are generally associated with volume constrained boundary conditions [2,18,46]. Some recent works in the literature regarding new techniques for enabling PD models to face unbounded domains can be found in [20,51,63,68] for 1D and 2D wave and in [71] for 1D diffusion problems.
Herein, we propose a new strategy for constructing ABCs, explicitly in the time domain for both the classical local and the PD diffusion formulation. First, we construct the ABCs for the classical local diffusion equation. Then, the approach is extended and applied to a PD diffusion model. The present study is inspired by two recent works conducted in [38,53] about a new type of ABCs for scalar and elastic wave propagation problems (hyperbolic PDEs) within the classical continuum framework. Here, similar to many other techniques, we truncate the unbounded domain at an artificial boundary and restrict the computational domain. The solution of the interior domain is computed by a standard FEM solver in the case of local diffusion and by a meshfree PD solver in the case of nonlocal diffusion; moreover, for marching in time, we make use of an explicit scheme. The ABCs are derived from a simple collocation approach that is local in time. The collocation employs a least-squares scheme over a subdomain centred at nodes to which the ABCs must be prescribed. The subdomain is referred to as cloud in the terminology of meshless methods. The collocation approach is simple and shares similarities with that of the meshless finite point method (FPM) [37,41,49,50,54]. We approximate the far field solution by a series of exponential basis functions (EBFs) which satisfy the governing equation of diffusion. EBFs are obtained as time-dependent fundamental solutions (modes) based on which a semi-analytical solution of the far field can be achieved. For this purpose, EBFs are adjusted so that they can play the role of vanishing modes capable of satisfying the corresponding conditions at infinity. The modes align the flux towards the exterior domain and prevent any unphysical reflection from the absorbing boundary. The unknown coefficients of the series can be found in terms of the nodal collocation values. In this way, the approximate solution provides Dirichlet ABCs for the corresponding nodes. In the case of the FEM, the corresponding nodes are located on the absorbing boundary, while in a PD model the corresponding nodes are located within a layer adjacent to the absorbing boundary. This will be explained in detail later.
The proposed approach offers several appealing features. First, the method does not interfere with the discretization scheme of the near field, and it is simple to implement. Second, as the method is local in time, it is free from Fourier and Laplace transform procedures; this feature is especially important for PD as, in this case, the calculation of the kernel function from the Laplace transform is complicated (due to the nonlocality) [71]. Third, the proposed ABCs are Dirichlet-type, so they can be simply prescribed at the corresponding nodes. Fourth, the method neither requires any spatial differentiation of the field variable nor any auxiliary variables involved in many high-order ABCs methods [69]. In the nonlocal diffusion case, the modes are dependent on the horizon. We shall show that they converge to the local modes in the limit when the horizon approaches zero. The extension of the presented approach to 3D problems is straightforward. To the best of our knowledge, our contribution addresses for the first time a PD model solved in an unbounded domain in 3D. In the section of numerical examples, we shall show that the method has a proper level of accuracy and performs well in the case of long-term calculations.
A brief outline of the rest of this paper is as follows. Section 2 describes the governing equations and physics of the problem, and it gives a brief overview of the PD formulation for diffusion-type problems. In Sect. 3 the way of deriving the fundamental solutions necessary for the far field solution is explained. The solution strategy of the method is thoroughly described in Sect. 4. The performance of the method, in terms of accuracy and stability, is investigated in Sect. 5. The final conclusions made throughout the paper are summarized in Sect. 6

Problem description
For simplicity, hereinafter we describe the problem for 2D cases, although the explanations are extensible to 1D and 3D.
Let us consider a generic 2D unbounded domain ⊂ R 2 of diffusion as shown in Fig. 1a. The medium is assumed to be isotropic, and it surrounds some physical objects (including baffles and sources). The governing diffusion equation, according to the classical local theory, for any point in the unbounded domain is: where u(x, t) is the concentration (the main field variable) of the diffusing material at point x with coordinates x = x, y and at time t, χ 2 is referred to as the diffusion coefficient, s(x, t) is a given source function (for instance a heat source), and ∇ 2 is the well-known Laplacian operator. Meanwhile, u = ∂u/∂t is the first time derivative of u. In fact, u is the diffusive concentration that can represent different physical quantities. For instance, u can represent the temperature in heat conduction or the mass concentration in mass transport such as the metal concentration in a corrosion problem. The surface of baffles are governed by either Dirichlet or Neu- in which u * is a function that gives the concentration values prescribed at the Dirichlet boundary D , n = n x , n y is an outward unit vector normal to the surface, and q * is a function that gives the values of the flux (the amount of concentration entering through the surface per unit area during a unit time interval) at the Neumann boundary N . The domain can be excited by an initial condition as follows: where u 0 is a given function. In fact, u 0 in (3) and s in (1) are both compactly supported. Likewise, since the medium is unbounded, the following condition holds: which indicates that the concentration vanishes at infinity. The main difficulty that arises out of the solution of (1) is the infinity of the medium. One way of applying the standard numerical methods, such as the FEM, is to discretize a very large portion of the unbounded domain . However, this strategy renders the approach inefficient and computationally expensive. A remedy is to bound the domain and restrict the computations to the region of interest. In this way, the exterior domain is truncated by an artificial boundary ∞ , so-called absorbing boundary, as shown in Fig. 1b. Therefore, is divided into a bounded domain N (near field or computational domain) and an exterior domain F = / N (far field). N is the region that includes all the physical objects.
∞ is placed so that the entire support of the sources in (1) is contained in N . Accordingly, one may conclude that the homogeneous counterpart of (1) for the exterior domain is: and the medium is initially at rest. Consequently, obtaining a well-posed problem in N requires imposing the so-called absorbing boundary conditions (ABCs) on ∞ , so that the effects of the far field are captured and the bounded domain reproduces the solution of the original problem. In this sense, the concentration (the diffusive quantity in N ) must be absorbed perfectly on ∞ and have a flux direction towards the exterior domain (no unphysical flux back to the near field is permitted). In view of (4), the flux should have a vanishing behaviour so that at infinity the concentration reaches a zero state. In the present work, we shall introduce a way, inspired by the study in [53], to construct time-dependent ABCs fulfilling the aforementioned requirements. The proposed method makes use of fundamental solutions (modes) that satisfy the governing equation of the exterior domain, expressed in (5), upon which Dirichlet ABCs are constructed. The way to construct the modes, for both classical and PD diffusion, is explained in Sect. 3, and the way of applying them to the far field solution is described in Sect. 4.2. Before that, in the following subsection we review the PD diffusion formulation.

Overview of the peridynamic diffusion formulation
In this section, the nonlocal formulation of PD for diffusiontype problems is briefly explained. One may refer to [9,10] for To construct the PD model, a region H x (see Fig. 2) around each point x of the domain is considered. H x is referred to as the neighborhood of x and is usually taken as a disk in 2D (a ball in 3D) centred at x. A concentration value u(x, t) is associated to the point x at time t. The neighborhood specifies the region of points x with which the point x establishes nonlocal interactions. In this way, each point is connected to its neighboring points through pipe-like conductors, diffusion bonds, which transfer the concentration among the points (buckets). In the bond-based PD formulation, the transport in ξ := x − x, the bond connecting x and x , is independent of the transport in other bonds. In the present study, we follow the bond-based PD formulation.
The governing equation of diffusion in PD, the nonlocal counterpart of (1), is given by [9,10] where J p is the kernel of the integral operator or the concentration flow density that takes the following form [72]: where δ is the horizon (radius of the disk), p is an integer that influences the shape of the kernel, usually taken equal to 0, 1, or 2 (the reader may refer to [14] for more details), and c is known as the micro-diffusivity in PD. c can take different forms and two popular ones are constant and linear [72]. The constant micro-diffusivity is given by which is independent of the bond length. In contrast, in the linear form, the micro-diffusivity is linearly dependent on the bond length and given by The constants c 0 and c 1 can be related to the diffusion coefficient, χ 2 . A way of correlating c with χ 2 is by matching the PD flux to that of the local diffusion as studied in [9]; we shall explain an alternative approach that ends up with similar results in the subsequent section. These relations are provided for different values of p in Table 1. Moreover, other forms for the micro-diffusivity are possible when considering more general influence functions in PD [48]. It is worthwhile to mention that (6) does not contain any spatial derivatives of the field variable (see (7)). As a consequence, the PD diffusion equation is valid in the whole domain including discontinuous regions.
Similar to the explanation given for (5), the governing equation of the far field in PD takes the following form: In PD models, the boundary conditions are usually imposed in volumetric regions rather than on surfaces which is the case for the local diffusion model (see Fig. 2). This is due to the nonlocality of the kernel and the associated neighborhood to each point. The common way of imposing Dirichlet boundary conditions is by prescribing the value of the corresponding boundary points [14]. However, for the Neumann boundary condition, the rate of concentration flow, denoted byQ, entering through the boundary surface should be first calculated as follows: where N is the surface to which the flux is applied. Then, Q must be appropriately distributed over the points in the boundary region and appears as part of the source term s(x, t) in (6). In fact, the effective way of imposing Neumann boundary conditions requires a separate and comprehensive discussion; for brevity, we refer the reader to the studies in [14,42,62]. It should be remarked here that for a comprehensive investigation into the rate of convergence of the nonlocal solution to that of the classical solution, taking different forms of micro-diffusivity and kernel functions, the reader may refer to [14]. This topic is beyond the scope of the present study.
There are different ways to discretize a PD model [23]. In the present study, we pursue a common way of PD discretization given by a meshfree scheme introduced in [57]. This scheme is easy to implement and frequently used by researchers. In this meshfree approach, the domain is discretized using a uniform grid with grid spacing x as shown in Fig. 3. The domain is represented by cells centred at the nodes and hence a volume is associated to each node (length x in 1D, area ( x) 2 in 2D, and volume ( x) 3 in 3D). In Fig. 3, the node x i interacts with all its family nodes x j within its neighborhood H x i . Moreover, the horizon is expressed as δ = m x, so that the value of m determines the number of family nodes. Likewise, to compute the integral over each neighborhood, we apply a one-point quadrature with quadrature points given by the nodes and quadrature weights given by the volumes of the intersections between the neighborhood and the cells. The discertized form of the governing equation in (6) can be then written aṡ where F i is the set of family nodes of x i ; moreover, ξ i j := x j − x i , V j is the associated volume to node x j , and β plays the role of a correction function through which the actual volume of x j covered by the neighborhood of x i can be determined. In the present study, we take β as recommended in [70]. Related works employing corrections known as partial volumes, which approximate the intersections between neighbor cells and the neighborhood of a given point to use as quadrature weights, have been presented in [45,47].

Fundamental solutions (modes)
The main purpose of this section is to find fundamental solutions (modes) using exponential basis functions (EBFs) for the governing equations of diffusion in the far field. First, we attempt to find the modes for the local diffusion equation. Let us consider ψ(x, t) as a mode that satisfies the governing Eq. (5). A general form of ψ can be written as where (α x , α y , ω) ∈ R 3 are three given parameters. To take control over the direction of the diffusion flux, without loss of generality, we make the following assumption: where κ affects the magnitude of the flux, while φ determines the flux direction and, for the time being, we assume that the flux direction may be oriented towards any angle of the unit circle. We shall discuss the proper selection of the parameters later. Let us substitute ψ from (13) for u in (5); then we have To reach a non-trivial solution, the expression in square brackets must vanish (since exp(α · x + ωt) = 0), i.e., which reduces to: This equation is referred to as the characteristic equation.
Replacing (17) into (13) leads to modes capable of constructing semi-analytical solutions for the local diffusion as follows: A similar procedure is performed for PD. We again consider a generic mode as in (13) to construct fundamental solutions for the nonlocal diffusion equation. Let us assume p = 1 in (7) and consider a constant micro-diffusivity. Then, insertion of ψ into (10) leads tȯ and in view of (13) using the change of variables ξ = x − x, it can be rewritten as where H is the neighbohood around the origin. To obtain a non-trivial solution, the expression inside the parentheses must vanish. Therefore, the characteristic equation of PD (the counterpart of (17)) can be obtained in the polar coordinate system, ξ = r (cos θ, sin θ), as follows: where η = α · cos θ, sin θ . By replacing ω from (21) into (13), the PD mode can be obtained. However, to make the expression simpler and find the connection to the local modes, we rewrite the expression in (21) taking a Taylor expansion of the integrand with respect to r η and about 0, as follows: which concludes: and thus the nonlocal PD mode can be written as A glance at the obtained modes in (18) and (24) reveals that the difference is in the temporal parts. By replacing The expressions for the constants c 0 in (8) and c 1 in (9), in terms of χ 2 and δ, for different choices of p for J p in (7) in 1D, 2D, and 3D 1D 2D 3D In the limit as κδ → 0, we obtain which indicates that the nonlocal mode recovers the local one provided that the micro-diffusivity is chosen appropriately in terms of the diffusion coefficient χ 2 and horizon δ.
In Table 1 the expressions of the micro-diffusivity constants, in terms of the diffusion coefficient of the local diffusion and the horizon, are presented for various kernels in 1D, 2D, and 3D. The way of finding the 3D modes is explained in Appendix A.

Remark 1
The modes are found so that the flux direction can be adjusted easily. We shall benefit from this feature to construct semi-analytical solutions for the far field (see Sect. 4.2). From (24) it can be inferred that by taking higher-order terms (inside the parentheses), i.e., keeping more terms from the Taylor expansion in (22), the semi-analytical solution should converge to the nonlocal solution and take more nonlocal effects into account. However, since the absorbing boundary is set far from the physical objects where the nonlocal effects can play an important role, for example, at the tip of a crack, even keeping only a few terms from the Taylor expansion, e.g., making use of the local modes, can still result in a proper approximation of the far field. Now, let us consider a circular region, so-called cloud, centred at a generic point x as shown in Fig. 4a. The cloud is surrounded by an unbounded domain; moreover, a local coordinate system is centred at the point x. In Sect. 4.2, we approximate the far field solution through a series of modes with a general form as in (18). In fact, the approximate solution must be obtained in terms of the values of the nodes within the cloud of the node x. From (18), it can be inferred that a mode consists of a temporal part and a spatial part. One may rewrite it in polar coordinates as follows: To enforce the infinity condition in the direction θ * , we substitute ψ for u in (4) and set θ = θ * , to obtain At a given time t > 0, provided that χ 2 κ 2 t < ∞, the infinity condition is satisfied as long as the following condition is met: The condition (30) reduces to two different conditions depending on the sign of κ and the value of θ * towards which the infinity condition is of concern. For instance, for the first quadrant of the local coordinate system, i.e., θ * ∈ [0, π/2), the following conditions must be met so that the mode can satisfy the condition (30) at the angle θ * : which indicates two halves of the unit circle as shown in Fig. 4b. In other words, the semicircular regions illustrate the permitted range of flux direction φ (based on the sign of κ) through which the infinity condition at θ * is satisfied. The interface between the two semicircles is perpendicular to the direction of θ * , and it rotates accordingly by changing the value of θ * (see Fig. 4c). Now let us assume that the infinity condition is always taking place in the first quadrant of the local coordinate system, i.e., θ * ∈ [0, π/2). Let us consider 1 and 2 as two valid regions of φ in the unit circle corresponding to θ * = 0 and θ * = π/2, respectively. Consequently, to satisfy the condition at any angle in the first quadrant of the coordinate system, φ must be chosen from the intersection of the regions, which is: With respect to the obtained regions in (32), in Sect. 4.2 we shall explain the proper application of the modes, the suitable orientation of the local coordinate system, and the way of patching the far field solution to that of the near field in a discretized PD or FEM model.

The solution strategy
The first part of the solution strategy concerns the near field approximation; in other words, employing a PD model discretized with a meshfree method or a local governing equation discretized by the FEM in space and proceeding in time with an explicit time marching algorithm. The second part, the main focus of the present study, corresponds to the approximation of the far field and construction of the ABCs. We assume that the absorbing boundary ∞ is a circle (or a sphere in 3D) surrounding the domain of interest. The details are explained in the following subsections.

Near field
Let us discretize a generic domain with a PD grid or a FEM mesh as shown in Fig. 5. The figure illustrates only a portion of the discretized domain close to the absorbing boundary. The nodes of the near field are indicated by black circles, while the absorbing nodes are indicated by blue circles. As can be seen in the figure, in the PD case the absorbing nodes are distributed within a certain distance of the absorbing boundary ∞ ; in contrast, in the FEM case only absorbing nodes along that boundary are required. In both cases, the domain is represented by a set of nodes x i , i = 1, 2, . . . , where x i = x i , y i contains the Cartesian coordinates of the node with respect to the global coordinate system. The time is discretized into instants, as t 0 , t 1 , . . . , t n , t n+1 , . . . . In the present study, we take a constant time step t = t n+1 − t n for any n = 0, 1, . . . . Moreover, u n i := u(x i , t n ) denotes the concentration value at the node x i at time instant t n = t 0 + n t. The time marching is performed by means of an explicit forward Euler time integration scheme. Having known the nodal values of x i at time instant t n , i.e., u n i and its first time derivativeu n i , the solver can proceed to the next time instant as follows: in which t must be less than a critical value t c obtained as where min is the minimum PD grid spacing or the minimum element size in FEM in the spatial discretization, and CFL stands for the well-known Courant-Friedrichs-Lewy condition (CFL is equal to 1/2, 1/4, and 1/8, respectively, for 1D, 2D, and 3D uniform discretizations). The value ofu n+1 i is calculated directly from the discretized form of the gov-erning equation. The values of the absorbing nodes must be updated at each time step by a strategy thoroughly described in the subsequent section.

Far field
Inspired by the study in [53], an efficient updating technique is devised for the absorbing nodes. In this way, at each time instant, the far field solution (within a cloud around each absorbing node) must be first constructed using the modes introduced in Sect. 3 as the bases of the approximation. Then, the approximate solution provides the value of the corresponding absorbing node for the next time instant, and it can be conveniently imposed to the node as a Dirichlet boundary condition. Figure 6 illustrates two generic clouds with two differently oriented local coordinate systems. In the figure, x i represents an absorbing node to which a cloud denoted by ∞ i is associated, and x j represents a neighboring node of the node x i . From the figure, it can be inferred that the local coordinate system is rotated by an angle τ (with respect to the global system) so that the first quadrant of the local coordinate system falls within the exterior domain and the bisection line of the 90 • between the axes is normal to the absorbing boundary. Within each cloud, we approximate the solution by u ∞ (x, t) using a series of modes with a form obtained in (18) as follows: where a k,l are unknown coefficients of the approximation. According to (32) as well as the orientation chosen for the local coordinate system of the cloud, we choose φ k from a and we take κ l from a set of equally spaced numbers with negative values, within the following interval: A logical way of choosing κ will be discussed in Sect. 4.3. Now the challenge is to find the unknown coefficients of the approximation in (35) in terms of the nodal values of the cloud. It should be pointed out that when the solver is at the beginning of the nth time step, i.e., the time interval between t n and t n+1 , we reset the time and assume that the approximation is valid for the local time interval [0, t]. This assumption is valid since here the approximation is local in time. To this end, in view of (35), we approximate the solution at the time interval between t n and t n+1 as follows: in which the series run over n b modes and q b are the corresponding unknown coefficients; moreover, ψ and q n are: q n = q n 1 , q n 2 , . . . , q n At the beginning of the nth time step the nodal values of the cloud ∞ i , associated to node x i as shown in Fig. 6 and according to (38), can be collected in a vector as (recalling that at the beginning of each step we sett = 0): where M stands for the moment matrix of the local collocation procedure, x i j denotes the local coordinates of x j , and R is a rotation matrix defined as for transformation from the global coordinate system to the local one. By assuming that the number of nodes in the cloud are less than the number of modes n b , one can find the unknown coefficients of the approximation q n in terms of the nodal values of the cloud U n as follows: (42) in which the superscript "+" denotes the Moore-Penrose generalized inverse. By replacing q n in (38) by the expression given in (42), one can write where N is a vector that contains the shape functions of the approximation and its size is equal to the number of nodes in the cloud. Equation (43) gives an approximation of the solution of the cloud explicitly in time taking the effects of the exterior domain to the near field into account. According to (43), one can simply update the nodal value of the absorbing node at the nth time step and proceed to the time instant t n+1 as follows: (44) in which V is referred to as the updating vector of the absorbing node; it should be recalled that the right-hand side of (44) is known since the solver is using an explicit time marching algorithm. In other words, (44) provides a Dirichlet ABC for the absorbing node in terms of the nodal values of its cloud.

Remark 2
Since a constant value for t is assumed and the local time is reset at each time step, the updating vectors V of the absorbing nodes consist of constant values and thus they can be calculated once at the beginning of the simulation. This leads to a substantial reduction in the computational cost of the approach.

Remark 3
In 3D problems, ∞ is the surface of a truncating sphere in space. In this case, based on the discussion in Appendix A, we assume that the first octant of the local coordinate system falls within the exterior domain. In order to achieve that, we rotate the local coordinate system so that the unit vector 1/ √ 3 1, 1, 1 becomes normal to the absorbing boundary and points toward the exterior domain.

Remark 4
In the present study, as shown in Fig. 5a, the thickness of the absorbing layer in the PD model is taken so that the neighborhoods of points in the near field have a full circular shape in 2D (spherical shape in 3D). Taking a narrower absorbing layer contributes to a well-known source of error in PD models, referred to as skin-or surface-effect [34,55], that may deteriorate the accuracy of the near field solution (which is the main concern for the analyst). We shall show this issue in the numerical examples.
In the subsequent section, some implementation aspects of the present study are discussed.

Implementation
Since the proposed method is devised preserving the discretization scheme of the near field, it has a simple and straightforward implementation. The main steps are discretizing the computational domain using a grid of nodes (in PD) or a mesh of elements (in the FEM); specifying the absorbing nodes; constructing the neighborhoods (only in PD) as well as the clouds of the absorbing nodes; forming the final system of equations; calculating the updating vectors of the clouds; and finally marching in time with the explicit algorithm described earlier.
For both the PD and FEM cases, it is possible to construct a final system of equations, based on their discretization schemes, as follows: whereM andK are the lumped mass and stiffness matrices of the global system of equations, respectively, U T andU T are two vectors collecting the nodal values and their first time derivative in the global system, andF is a vector of known values. In the present study, the system of equations for both the bounded-domain problem and the unboundeddomain with ABCs problem is identical. However, for the unbounded-domain problem an extra procedure to update the nodal values of the absorbing nodes is required. One may refer to Algorithm 1 to get more insight into the step-by-step solution procedure of the present study. In the algorithm, d ∞ refers to the radius of the clouds centred at the absorbing nodes and t f is the final time instant of the solution.
As discussed earlier, we construct the modes taking φ and κ from the intervals expressed in (36) and (37), respectively.

Algorithm 1
The step-by-step solution process. 1: Specify the computational domain 2: Discretize the computational domain using a PD grid (and construct the corresponding neighborhoods) or a FEM mesh 3: Construct the matricesK andM and vectorF in (45) 4: Identify the absorbing nodes on ∞ or in F 5: for all x i absorbing nodes do 6: Find nodes x j such that Fig. 6) 7: Determine the rotation matrix in (41) 8: Calculate and store V for ∞ i using (44) (see Remark. 2) 9: end for 10: Initialize U 0 T from the initial conditions in (3)  In fact, κ is an influential parameter for the stability of the numerical solution. This is due to the fact that it directly influences the magnitude of the temporal part of the modes. According to (29), the temporal part should be bounded so that the satisfaction of the infinity condition can be guaranteed. Moreover, taking a large value for κ may give rise to exponentially big numbers for the moment matrix M in (40) and contribute to significant numerical errors. To confine the magnitude of the temporal part, we select κ so that the following condition is met (recalling that, at each time step,t = t is plugged into (43)): In view of (34), the following suitable range for κ is obtained: We note that the collocation in each cloud must be taken over a sufficiently large number of neighbouring nodes to ensure a proper approximation of the far field solution. Based on our experience the radius d ∞ must be large enough so that the cloud of each absorbing node in 1D, 2D, and 3D includes at least 3, 5, and 7 neighbouring nodes, respectively. It should be pointed out that d ∞ determines solely the region of the least-squares approximation and does not convey the range of nonlocal interaction as done by δ.
In this study, for the sake of computer implementation we have used the C++ programming language. To increase the performance of the developed code, the major parts have been parallelized using OpenMP directives. The library libkdtree has been employed for constructing the neighborhoods and the clouds in a high-performance way; the reader may refer to [39] for more details. In order to perform the Moore-Penrose generalized inverse in (42), we have exploited the well-known library LAPACK [1].

Numerical examples
In this section the performance of the proposed method in the solution of some problems is investigated. In the last two examples, we apply the proposed method to a PD corrosion problem. In this way, we demonstrate one of the appealing applications of the method. The way of tackling corrosion with PD corresponds to the near field solution. To this respect, we employ the PD model proposed in [29] for the near field solution. One may refer to [13] for more details about the parameters and the algorithms required to treat corrosion in PD.
To scrutinize the performance of the method, in terms of accuracy and stability, we define the energy of the solution at time t as follows: where N is the near field. In fact, provides a quantity to study the behaviour of the solution. In this sense, the energy must be a non-increasing function of time. It should be pointed out that in the first two examples, we take the solution of an extremely large bounded domain with homogeneous Neumann boundary conditions as the reference solution. The extended-domain solution can play the role of a reference solution as long as the boundary region of the extended domain is not perturbed and it remains at a zero state of concentration. In all the examples, for the PD model, we take p = 1 in the kernel of the integral operator (7) and the micro-diffusivity c is assumed to be constant as in (8).
It should be pointed out that all the examples are set up in a dimensionless system; moreover, the chosen values for the geometrical and physical parameters of the problems are not motivated by any real physical application.

Example 1
The main goal of this example is to show the performance of the proposed method in the case of a 1D unbounded problem domain (−∞, ∞). The diffusion parameter of the medium is χ 2 = 1. We consider the region (−30, 30) as the computational domain N . To validate the obtained results, the  −150, 150). The problem domain is excited by an initial concentration as follows: We begin with the solution given by a FEM model using linear 1D elements; then, we look into the solution of a PD model. For both models we take an identical discretization resolution in space and time. For the spatial discretization, we take x = 0.125 which results in 2401 and 481 nodes, respectively, for the extended domain and the truncated domain (or computational domain). The solution for a time duration t f = 250 is of concern; this duration is set so that the concentration close to the boundaries of the extended domain remains almost zero. The time step t is set to be 0.001 and thus the calculations are performed over 250,000 time steps. In the FEM model, the ABCs are constructed by placing only two clouds on the nodes positioned at x = −30 and x = 30. The local coordinate system of the cloud at x = −30 is oriented so that its axis points towards the left direction, whereas the one at x = 30 points towards the right direction (similar to that of the global coordinate system). The construction of the far field solution is carried out using n b = 20 modes. Figure 7 illustrates the concentration field obtained at different time instants for the truncated domain and the extended domain using the FEM solver. For a better comparison between the solutions, only the part of the reference solution coinciding with the truncated domain is shown in the figure. It can be inferred that the solution obtained by the proposed method conforms to the reference solution very well. Meanwhile, the solution of the truncated domain is continued further up to t f = 4,000 and, as expected, the concentration reaches a zero state. To get more insight into the performance of the proposed method, the variation of  Fig. 8. The energy is calculated within the near field in both cases. The energy is normalized with respect to the energy of the system at the first time instant. Moreover, the normalized energy difference between the solutions is defined as follows: where ref (t) is the energy of the reference solution at time t. The results are reported in Fig. 9. In the same figure, for the sake of comparison, we report the error when the truncated domain is not equipped with the ABCs, but instead homogeneous Neumann boundary conditions are imposed. The results in Fig. 8 indicate that the energy dissipates from the near field monotonically (d < 0) having an excellent match with the reference solution. As can be seen in Fig. 9, the difference is very small as long as the truncated domain is equipped with the ABCs, while this is not the case when the homogeneous Neumann boundary conditions are employed in which case the difference rapidly grows in time. It can be concluded that the proposed ABCs are performing well and no accountable spurious flux reflecting back to the near field from the boundaries can be observed. As shown in Fig. 8, in order to check the stability of the present approach, the energy is calculated up to t f = 4,000 which results in 4 × 10 6 time steps. It can be inferred that the method performs stably in time and no instability due to long-term calculations takes place.
In the second part of this example, we redone the procedure with a PD solver. The horizon is taken as δ = 4 x and to the absorbing nodes a cloud with d ∞ = δ/2 is associated. Based on the discussion in Remark 4, the absorbing layer thickness is taken as δ L = δ; the layer is wide enough so that each neighborhood in the near field has a full integration domain.  The concentration field obtained for the truncated domain and the extended domain using the PD solver at different time instants is depicted in Fig. 10. Moreover, the variation of energy versus the time step number, for both the truncated domain and the extended domain, is reported in Fig. 11. Similar to the FEM case, applying the proposed ABCs to the PD  To investigate the effects of the absorbing layer thickness δ L on the accuracy of the truncated-domain solution, we redo the simulation with the PD solver taking different layer thicknesses equal to 3 x, 2 x, and x. Figure 12 illustrates the concentration field at t = 250 taking different values for δ L . The results indicate that taking a smaller value for δ L , and consequently increasing the surface effect, results in a lower accuracy of the solution. It can be observed that the mismatch between the numerical solution and the reference solution is more pronounced near the boundaries, which is consistent with the discussion in Remark 4. In Fig. 13, the normalize energy difference between the solutions of the truncated domain, using different sizes of the absorbing layer, and the reference solution is reported. As can be seen, the method exhibits results with an appropriate level of accuracy provided that the absorbing layer is taken wide enough, i.e., δ L ≥ 3 x.

Example 2
The purpose of this example is to scrutinize the performance of the proposed method in the solution of a 2D corrosion problem in an unbounded domain. This problem concerns the corrosion of a metallic specimen (shown in Fig. 14) surrounded by an unbounded electrolyte that leads to reduction of the solid material due to dissolution into the liquid medium. In this example, the intact metal (the "HZG" 1 letters with initial concentration value u 0 = 2) is the initial solid phase and the surrounding electrolyte (with initial concentration value u 0 = 0) is the initial liquid phase of the problem domain. During the corrosion process, due to the so-called dissolution flux, the solid phase reduces and, consequently, the dissolved material diffuses out into the electrolyte. The whole process can be modelled by a PD diffusion solver. It should be pointed out that the solution of such a problem by the standard classical models is cumbersome since the interfacial dissolution flux (between the phases) cannot be defined easily in the framework of local diffusion [28]. This roots from the fact that the concentration field is not smooth at the corrosion front; therefore, capturing the existing jumps (strong discontinuities) is not an easy task for the local models. Furthermore, in this example the solid phase consists of sharp corners (weak discontinuities) which make the application of the local models more complicated.
For the PD model, it is required to associate two microdiffusivity parameters c liq and c int . The latter is referred to as the micro-dissolvability [30] and corresponds to the bonds with one end in the solid phase and one end in the liquid phase (interfacial bonds). The former corresponds to the bonds whose endings are in the liquid phase. It is assumed that no flux between two points in the solid phase occurs. It is also necessary to define a threshold concentration u sat that stands for the minimum concentration of the solid phase during the dissolution. Once the dissolution process (diffusion through the interfacial bonds) drops the concentration of a solid point below u sat , the point is regarded as a liquid point afterwards. This phase-change results in autonomous propagation of the dissolution front.
In this example it is assumed that δ = 0.04, u sat = 0.8, c liq = 7957.70, and c int = 795.77. As shown in Fig. 14, a circular domain with a radius equal to 2 is considered for the truncated domain. The domain is surrounded by an absorbing layer with a thickness equal to the horizon δ. To each absorbing node a cloud with a radius d ∞ = δ/2 is associated. In each cloud, the far field is approximated by n b = 55 modes. For the sake of verification, the solution of an extended circular domain with a radius equal to 10 is taken as the reference solution. Homogeneous Neumann boundary conditions are imposed on the boundary of the extended domain. Both domains are discretized with the same grid spacing, x = 0.02; therefore, the calculations runs over 31,417 nodes for the truncated domain (including 2500 absorbing nodes), while they run over 785,349 nodes for the extended domain. To proceed in time, in both models an incremental time step t = 3.3 × 10 −7 is applied.
In Fig. 15, the concentration level of metal captured by the truncated domain using ABCs and the extended domain at different time instants is presented. The regions with the concentration level higher than 0.8 represent the solid phase of the solution domain. The obtained results demonstrate that the solution of the truncated domain resembles that of the extended domain The metal concentration begins to diffuse out and crosses the truncating boundary with a pattern very similar to that of the reference solution. By t = 0.085 the concentration level of the near field reaches almost a zero state, which indicates that the ABCs are performing properly.
Another way to demonstrate that the applied ABCs are performing appropriately is to check the gradient field of concentration. Contour plots of the gradient field of concentration at different time instants are presented in Fig. 16. The arrows show the direction of gradient vectors while the contours indicate the magnitude. It should be pointed out that the To get a better understanding of the accuracy of the solution, the difference between the solution obtained for the truncated domain and the extended domain, on the basis of (50), is reported in Fig. 17. For comparison, in the same plot, the results of the truncated domain with homogeneous Neumann boundary conditions (when the domain is not equipped with the absorbing layer) is shown. It can be concluded that the difference between the solution of the truncated domain and the extended domain is relatively small as long as the absorbing layer is applied, while the solution of the truncated domain diverges from the reference solution when the absorbing layer is removed. To check the stability of the solution, the variation of energy obtained by the truncated domain and the extended domain is reported in Fig. 18. For the truncated domain the calculations are run up to almost 7×10 6 time steps. It can be inferred that the absorbing boundary layer is performing with a stable behaviour in time and the energy dissipates monotonically from the computational domain in a very good agreement with the reference solution.

Example 3
The goal of this example is to demonstrate that the proposed method is applicable to 3D problems. For this purpose, we simulate the dissolution process of a spherical metallic specimen (with an initial concentration value u 0 = 1) in an unbounded domain occupied by an electrolyte (with an initial concentration value u 0 = 0).
As shown in Fig. 19, the initial solid phase has a radius equal to 3.33 and the domain is truncated by a spherical domain with a radius equal to 5. The required parameters of the PD corrosion model are taken as δ = 0.2, u sat = 0.8, Fig. 18 The variation of energy in Example 2

Fig. 19
The initial concentration field in Example 3 c int = 795.77, and c liq = 7, 957.7. Similar to previous examples, the thickness of the absorbing layer is equal to the horizon δ and to each absorbing node a cloud with a radius d ∞ = δ/2 is assigned. In each cloud, the far field is approximated by n b = 605 modes. The domain is discretized with a grid spacing x = δ/2, which results in 523,155 nodes for the computational domain (including 60,374 absorbing nodes). For marching in time, an incremental time step t = 3.33 × 10 −8 is used. Since construction of a reference solution for this example requires an extremely large computational domain, the solution is validated qualitatively. Figure 20 illustrates the metal concentration level as well as the gradient field of the solution at three different time instants. Due to symmetry, we only show a portion of the solution domain in the first octant of the global coordinate system. From the concentration plots, it is clear that the solid phase is gradually reducing in time and the corrosion front is moving towards the centre of the sphere. Moreover, the results obtained for the gradient field indicates that the solid concentration is diffusing out of the domain with a flux perpendicular to the truncating surface and towards the exterior domain. It can be concluded that the absorbing layer is performing properly. For further investigation, the variation of energy in time is reported in Fig. 21. It can be inferred that the energy is dissipating monotonically in time. The calculations were done up to more than 2 × 10 6 time steps. No indication of instability was observed.

Conclusions
This paper proposes a new strategy for constructing absorbing boundary conditions (ABCs) suitable for diffusion-type problems in unbounded domains. A particular attention is paid to develop the ABCs for a peridynamic (PD) diffusion model. We develop the ABCs for both local and nonlocal diffusion at the continuum level and we devise an appropriate implementation of the method at the discrete level. The method is applicable to available meshfree PD discretizations for the nonlocal diffusion as well as the standard finite element method (FEM) for the local diffusion equation. The proposed ABCs are derived from exponential basis functions (EBFs) which play the role of vanishing modes capable of constructing semi-analytical solutions for the exterior domain. The way of adjusting the EBFs for proper satisfaction of the infinity condition is thoroughly discussed in the paper. Following that, the far field solution is patched to the near field one by prescribing accurate Dirichlet ABCs obtained from the semi-analytical solution. Since the ABCs are Dirichlet-type, the method is simple to implement and free from any approximating differential operator. The ABCs are developed explicitly in the time domain and thus the method is free from any Fourier and Laplace transform procedures. In fact, this feature is especially appealing for the PD nonlocal diffusion equation, since the calculation of the kernel function (involved in the corresponding nonlocal operators) from the Laplace transform is cumbersome. The performance of the method is studied through some examples, including a 3D corrosion problem. Our investigations demonstrate that the method produces accurate results and exhibits a stable behavior even in the case of long-term computations.